src/daComposant/daCore/BasicObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2021 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 """
  24     Définit les outils généraux élémentaires.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27 __all__ = []
  28
  29 import os
  30 import sys
  31 import logging
  32 import copy
  33 import time
  34 import numpy
  35 from functools import partial
  36 from daCore import Persistence, PlatformInfo, Interfaces
  37 from daCore import Templates
  38
  39 # ==============================================================================
  40 class CacheManager(object):
  41     """
  42     Classe générale de gestion d'un cache de calculs
  43     """
  44     def __init__(self,
  45                  toleranceInRedundancy = 1.e-18,
  46                  lenghtOfRedundancy    = -1,
  47                 ):
  48         """
  49         Les caractéristiques de tolérance peuvent être modifiées à la création.
  50         """
  51         self.__tolerBP   = float(toleranceInRedundancy)
  52         self.__lenghtOR  = int(lenghtOfRedundancy)
  53         self.__initlnOR  = self.__lenghtOR
  54         self.__seenNames = []
  55         self.__enabled   = True
  56         self.clearCache()
  57
  58     def clearCache(self):
  59         "Vide le cache"
  60         self.__listOPCV  = []
  61         self.__seenNames = []
  62
  63     def wasCalculatedIn(self, xValue, oName="" ):
  64         "Vérifie l'existence d'un calcul correspondant à la valeur"
  65         __alc = False
  66         __HxV = None
  67         if self.__enabled:
  68             for i in range(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
  69                 if not hasattr(xValue, 'size'):
  70                     pass
  71                 elif (str(oName) != self.__listOPCV[i][3]):
  72                     pass
  73                 elif (xValue.size != self.__listOPCV[i][0].size):
  74                     pass
  75                 elif (numpy.ravel(xValue)[0] - self.__listOPCV[i][0][0]) > (self.__tolerBP * self.__listOPCV[i][2] / self.__listOPCV[i][0].size):
  76                     pass
  77                 elif numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < (self.__tolerBP * self.__listOPCV[i][2]):
  78                     __alc  = True
  79                     __HxV = self.__listOPCV[i][1]
  80                     break
  81         return __alc, __HxV
  82
  83     def storeValueInX(self, xValue, HxValue, oName="" ):
  84         "Stocke pour un opérateur o un calcul Hx correspondant à la valeur x"
  85         if self.__lenghtOR < 0:
  86             self.__lenghtOR = 2 * min(xValue.size, 50) + 2 # 2 * xValue.size + 2
  87             self.__initlnOR = self.__lenghtOR
  88             self.__seenNames.append(str(oName))
  89         if str(oName) not in self.__seenNames: # Etend la liste si nouveau
  90             self.__lenghtOR += 2 * min(xValue.size, 50) + 2 # 2 * xValue.size + 2
  91             self.__initlnOR += self.__lenghtOR
  92             self.__seenNames.append(str(oName))
  93         while len(self.__listOPCV) > self.__lenghtOR:
  94             self.__listOPCV.pop(0)
  95         self.__listOPCV.append( (
  96             copy.copy(numpy.ravel(xValue)), # 0 Previous point
  97             copy.copy(HxValue),             # 1 Previous value
  98             numpy.linalg.norm(xValue),      # 2 Norm
  99             str(oName),                     # 3 Operator name
 100             ) )
 101
 102     def disable(self):
 103         "Inactive le cache"
 104         self.__initlnOR = self.__lenghtOR
 105         self.__lenghtOR = 0
 106         self.__enabled  = False
 107
 108     def enable(self):
 109         "Active le cache"
 110         self.__lenghtOR = self.__initlnOR
 111         self.__enabled  = True
 112
 113 # ==============================================================================
 114 class Operator(object):
 115     """
 116     Classe générale d'interface de type opérateur simple
 117     """
 118     NbCallsAsMatrix = 0
 119     NbCallsAsMethod = 0
 120     NbCallsOfCached = 0
 121     CM = CacheManager()
 122     #
 123     def __init__(self,
 124         name                 = "GenericOperator",
 125         fromMethod           = None,
 126         fromMatrix           = None,
 127         avoidingRedundancy   = True,
 128         inputAsMultiFunction = False,
 129         enableMultiProcess   = False,
 130         extraArguments       = None,
 131         ):
 132         """
 133         On construit un objet de ce type en fournissant, à l'aide de l'un des
 134         deux mots-clé, soit une fonction ou un multi-fonction python, soit une
 135         matrice.
 136         Arguments :
 137         - name : nom d'opérateur
 138         - fromMethod : argument de type fonction Python
 139         - fromMatrix : argument adapté au constructeur numpy.matrix
 140         - avoidingRedundancy : booléen évitant (ou pas) les calculs redondants
 141         - inputAsMultiFunction : booléen indiquant une fonction explicitement
 142           définie (ou pas) en multi-fonction
 143         - extraArguments : arguments supplémentaires passés à la fonction de
 144           base et ses dérivées (tuple ou dictionnaire)
 145         """
 146         self.__name      = str(name)
 147         self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
 148         self.__AvoidRC   = bool( avoidingRedundancy )
 149         self.__inputAsMF = bool( inputAsMultiFunction )
 150         self.__mpEnabled = bool( enableMultiProcess )
 151         self.__extraArgs = extraArguments
 152         if   fromMethod is not None and self.__inputAsMF:
 153             self.__Method = fromMethod # logtimer(fromMethod)
 154             self.__Matrix = None
 155             self.__Type   = "Method"
 156         elif fromMethod is not None and not self.__inputAsMF:
 157             self.__Method = partial( MultiFonction, _sFunction=fromMethod, _mpEnabled=self.__mpEnabled)
 158             self.__Matrix = None
 159             self.__Type   = "Method"
 160         elif fromMatrix is not None:
 161             self.__Method = None
 162             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
 163             self.__Type   = "Matrix"
 164         else:
 165             self.__Method = None
 166             self.__Matrix = None
 167             self.__Type   = None
 168
 169     def disableAvoidingRedundancy(self):
 170         "Inactive le cache"
 171         Operator.CM.disable()
 172
 173     def enableAvoidingRedundancy(self):
 174         "Active le cache"
 175         if self.__AvoidRC:
 176             Operator.CM.enable()
 177         else:
 178             Operator.CM.disable()
 179
 180     def isType(self):
 181         "Renvoie le type"
 182         return self.__Type
 183
 184     def appliedTo(self, xValue, HValue = None, argsAsSerie = False, returnSerieAsArrayMatrix = False):
 185         """
 186         Permet de restituer le résultat de l'application de l'opérateur à une
 187         série d'arguments xValue. Cette méthode se contente d'appliquer, chaque
 188         argument devant a priori être du bon type.
 189         Arguments :
 190         - les arguments par série sont :
 191             - xValue : argument adapté pour appliquer l'opérateur
 192             - HValue : valeur précalculée de l'opérateur en ce point
 193         - argsAsSerie : indique si les arguments sont une mono ou multi-valeur
 194         """
 195         if argsAsSerie:
 196             _xValue = xValue
 197             _HValue = HValue
 198         else:
 199             _xValue = (xValue,)
 200             if HValue is not None:
 201                 _HValue = (HValue,)
 202             else:
 203                 _HValue = HValue
 204         PlatformInfo.isIterable( _xValue, True, " in Operator.appliedTo" )
 205         #
 206         if _HValue is not None:
 207             assert len(_xValue) == len(_HValue), "Incompatible number of elements in xValue and HValue"
 208             _HxValue = []
 209             for i in range(len(_HValue)):
 210                 _HxValue.append( numpy.asmatrix( numpy.ravel( _HValue[i] ) ).T )
 211                 if self.__AvoidRC:
 212                     Operator.CM.storeValueInX(_xValue[i],_HxValue[-1],self.__name)
 213         else:
 214             _HxValue = []
 215             _xserie = []
 216             _hindex = []
 217             for i, xv in enumerate(_xValue):
 218                 if self.__AvoidRC:
 219                     __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xv,self.__name)
 220                 else:
 221                     __alreadyCalculated = False
 222                 #
 223                 if __alreadyCalculated:
 224                     self.__addOneCacheCall()
 225                     _hv = __HxV
 226                 else:
 227                     if self.__Matrix is not None:
 228                         self.__addOneMatrixCall()
 229                         _xv = numpy.matrix(numpy.ravel(xv)).T
 230                         _hv = self.__Matrix * _xv
 231                     else:
 232                         self.__addOneMethodCall()
 233                         _xserie.append( xv )
 234                         _hindex.append(  i )
 235                         _hv = None
 236                 _HxValue.append( _hv )
 237             #
 238             if len(_xserie)>0 and self.__Matrix is None:
 239                 if self.__extraArgs is None:
 240                     _hserie = self.__Method( _xserie ) # Calcul MF
 241                 else:
 242                     _hserie = self.__Method( _xserie, self.__extraArgs ) # Calcul MF
 243                 if not hasattr(_hserie, "pop"):
 244                     raise TypeError("The user input multi-function doesn't seem to return sequence results, behaving like a mono-function. It has to be checked.")
 245                 for i in _hindex:
 246                     _xv = _xserie.pop(0)
 247                     _hv = _hserie.pop(0)
 248                     _HxValue[i] = _hv
 249                     if self.__AvoidRC:
 250                         Operator.CM.storeValueInX(_xv,_hv,self.__name)
 251         #
 252         if returnSerieAsArrayMatrix:
 253             _HxValue = numpy.stack([numpy.ravel(_hv) for _hv in _HxValue], axis=1)
 254         #
 255         if argsAsSerie: return _HxValue
 256         else:           return _HxValue[-1]
 257
 258     def appliedControledFormTo(self, paires, argsAsSerie = False, returnSerieAsArrayMatrix = False):
 259         """
 260         Permet de restituer le résultat de l'application de l'opérateur à des
 261         paires (xValue, uValue). Cette méthode se contente d'appliquer, son
 262         argument devant a priori être du bon type. Si la uValue est None,
 263         on suppose que l'opérateur ne s'applique qu'à xValue.
 264         Arguments :
 265         - paires : les arguments par paire sont :
 266             - xValue : argument X adapté pour appliquer l'opérateur
 267             - uValue : argument U adapté pour appliquer l'opérateur
 268         - argsAsSerie : indique si l'argument est une mono ou multi-valeur
 269         """
 270         if argsAsSerie: _xuValue = paires
 271         else:           _xuValue = (paires,)
 272         PlatformInfo.isIterable( _xuValue, True, " in Operator.appliedControledFormTo" )
 273         #
 274         if self.__Matrix is not None:
 275             _HxValue = []
 276             for paire in _xuValue:
 277                 _xValue, _uValue = paire
 278                 _xValue = numpy.matrix(numpy.ravel(_xValue)).T
 279                 self.__addOneMatrixCall()
 280                 _HxValue.append( self.__Matrix * _xValue )
 281         else:
 282             _xuArgs = []
 283             for paire in _xuValue:
 284                 _xValue, _uValue = paire
 285                 if _uValue is not None:
 286                     _xuArgs.append( paire )
 287                 else:
 288                     _xuArgs.append( _xValue )
 289             self.__addOneMethodCall( len(_xuArgs) )
 290             if self.__extraArgs is None:
 291                 _HxValue = self.__Method( _xuArgs ) # Calcul MF
 292             else:
 293                 _HxValue = self.__Method( _xuArgs, self.__extraArgs ) # Calcul MF
 294         #
 295         if returnSerieAsArrayMatrix:
 296             _HxValue = numpy.stack([numpy.ravel(_hv) for _hv in _HxValue], axis=1)
 297         #
 298         if argsAsSerie: return _HxValue
 299         else:           return _HxValue[-1]
 300
 301     def appliedInXTo(self, paires, argsAsSerie = False, returnSerieAsArrayMatrix = False):
 302         """
 303         Permet de restituer le résultat de l'application de l'opérateur à une
 304         série d'arguments xValue, sachant que l'opérateur est valable en
 305         xNominal. Cette méthode se contente d'appliquer, son argument devant a
 306         priori être du bon type. Si l'opérateur est linéaire car c'est une
 307         matrice, alors il est valable en tout point nominal et xNominal peut
 308         être quelconque. Il n'y a qu'une seule paire par défaut, et argsAsSerie
 309         permet d'indiquer que l'argument est multi-paires.
 310         Arguments :
 311         - paires : les arguments par paire sont :
 312             - xNominal : série d'arguments permettant de donner le point où
 313               l'opérateur est construit pour être ensuite appliqué
 314             - xValue : série d'arguments adaptés pour appliquer l'opérateur
 315         - argsAsSerie : indique si l'argument est une mono ou multi-valeur
 316         """
 317         if argsAsSerie: _nxValue = paires
 318         else:           _nxValue = (paires,)
 319         PlatformInfo.isIterable( _nxValue, True, " in Operator.appliedInXTo" )
 320         #
 321         if self.__Matrix is not None:
 322             _HxValue = []
 323             for paire in _nxValue:
 324                 _xNominal, _xValue = paire
 325                 _xValue = numpy.matrix(numpy.ravel(_xValue)).T
 326                 self.__addOneMatrixCall()
 327                 _HxValue.append( self.__Matrix * _xValue )
 328         else:
 329             self.__addOneMethodCall( len(_nxValue) )
 330             if self.__extraArgs is None:
 331                 _HxValue = self.__Method( _nxValue ) # Calcul MF
 332             else:
 333                 _HxValue = self.__Method( _nxValue, self.__extraArgs ) # Calcul MF
 334         #
 335         if returnSerieAsArrayMatrix:
 336             _HxValue = numpy.stack([numpy.ravel(_hv) for _hv in _HxValue], axis=1)
 337         #
 338         if argsAsSerie: return _HxValue
 339         else:           return _HxValue[-1]
 340
 341     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue", argsAsSerie = False):
 342         """
 343         Permet de renvoyer l'opérateur sous la forme d'une matrice
 344         """
 345         if self.__Matrix is not None:
 346             self.__addOneMatrixCall()
 347             mValue = [self.__Matrix,]
 348         elif not isinstance(ValueForMethodForm,str) or ValueForMethodForm != "UnknownVoidValue": # Ne pas utiliser "None"
 349             mValue = []
 350             if argsAsSerie:
 351                 self.__addOneMethodCall( len(ValueForMethodForm) )
 352                 for _vfmf in ValueForMethodForm:
 353                     mValue.append( numpy.matrix( self.__Method(((_vfmf, None),)) ) )
 354             else:
 355                 self.__addOneMethodCall()
 356                 mValue = self.__Method(((ValueForMethodForm, None),))
 357         else:
 358             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
 359         #
 360         if argsAsSerie: return mValue
 361         else:           return mValue[-1]
 362
 363     def shape(self):
 364         """
 365         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
 366         la forme d'une matrice
 367         """
 368         if self.__Matrix is not None:
 369             return self.__Matrix.shape
 370         else:
 371             raise ValueError("Matrix form of the operator is not available, nor the shape")
 372
 373     def nbcalls(self, which=None):
 374         """
 375         Renvoie les nombres d'évaluations de l'opérateur
 376         """
 377         __nbcalls = (
 378             self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
 379             self.__NbCallsAsMatrix,
 380             self.__NbCallsAsMethod,
 381             self.__NbCallsOfCached,
 382             Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
 383             Operator.NbCallsAsMatrix,
 384             Operator.NbCallsAsMethod,
 385             Operator.NbCallsOfCached,
 386             )
 387         if which is None: return __nbcalls
 388         else:             return __nbcalls[which]
 389
 390     def __addOneMatrixCall(self):
 391         "Comptabilise un appel"
 392         self.__NbCallsAsMatrix   += 1 # Decompte local
 393         Operator.NbCallsAsMatrix += 1 # Decompte global
 394
 395     def __addOneMethodCall(self, nb = 1):
 396         "Comptabilise un appel"
 397         self.__NbCallsAsMethod   += nb # Decompte local
 398         Operator.NbCallsAsMethod += nb # Decompte global
 399
 400     def __addOneCacheCall(self):
 401         "Comptabilise un appel"
 402         self.__NbCallsOfCached   += 1 # Decompte local
 403         Operator.NbCallsOfCached += 1 # Decompte global
 404
 405 # ==============================================================================
 406 class FullOperator(object):
 407     """
 408     Classe générale d'interface de type opérateur complet
 409     (Direct, Linéaire Tangent, Adjoint)
 410     """
 411     def __init__(self,
 412                  name             = "GenericFullOperator",
 413                  asMatrix         = None,
 414                  asOneFunction    = None, # 1 Fonction
 415                  asThreeFunctions = None, # 3 Fonctions in a dictionary
 416                  asScript         = None, # 1 or 3 Fonction(s) by script
 417                  asDict           = None, # Parameters
 418                  appliedInX       = None,
 419                  extraArguments   = None,
 420                  avoidRC          = True,
 421                  inputAsMF        = False,# Fonction(s) as Multi-Functions
 422                  scheduledBy      = None,
 423                  toBeChecked      = False,
 424                  ):
 425         ""
 426         self.__name      = str(name)
 427         self.__check     = bool(toBeChecked)
 428         self.__extraArgs = extraArguments
 429         #
 430         self.__FO        = {}
 431         #
 432         __Parameters = {}
 433         if (asDict is not None) and isinstance(asDict, dict):
 434             __Parameters.update( asDict )
 435         # Priorité à EnableMultiProcessingInDerivatives=True
 436         if "EnableMultiProcessing" in __Parameters and __Parameters["EnableMultiProcessing"]:
 437             __Parameters["EnableMultiProcessingInDerivatives"] = True
 438             __Parameters["EnableMultiProcessingInEvaluation"]  = False
 439         if "EnableMultiProcessingInDerivatives"  not in __Parameters:
 440             __Parameters["EnableMultiProcessingInDerivatives"]  = False
 441         if __Parameters["EnableMultiProcessingInDerivatives"]:
 442             __Parameters["EnableMultiProcessingInEvaluation"]  = False
 443         if "EnableMultiProcessingInEvaluation"  not in __Parameters:
 444             __Parameters["EnableMultiProcessingInEvaluation"]  = False
 445         if "withIncrement" in __Parameters: # Temporaire
 446             __Parameters["DifferentialIncrement"] = __Parameters["withIncrement"]
 447         #
 448         if asScript is not None:
 449             __Matrix, __Function = None, None
 450             if asMatrix:
 451                 __Matrix = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
 452             elif asOneFunction:
 453                 __Function = { "Direct":Interfaces.ImportFromScript(asScript).getvalue( "DirectOperator" ) }
 454                 __Function.update({"useApproximatedDerivatives":True})
 455                 __Function.update(__Parameters)
 456             elif asThreeFunctions:
 457                 __Function = {
 458                     "Direct" :Interfaces.ImportFromScript(asScript).getvalue( "DirectOperator" ),
 459                     "Tangent":Interfaces.ImportFromScript(asScript).getvalue( "TangentOperator" ),
 460                     "Adjoint":Interfaces.ImportFromScript(asScript).getvalue( "AdjointOperator" ),
 461                     }
 462                 __Function.update(__Parameters)
 463         else:
 464             __Matrix = asMatrix
 465             if asOneFunction is not None:
 466                 if isinstance(asOneFunction, dict) and "Direct" in asOneFunction:
 467                     if asOneFunction["Direct"] is not None:
 468                         __Function = asOneFunction
 469                     else:
 470                         raise ValueError("The function has to be given in a dictionnary which have 1 key (\"Direct\")")
 471                 else:
 472                     __Function = { "Direct":asOneFunction }
 473                 __Function.update({"useApproximatedDerivatives":True})
 474                 __Function.update(__Parameters)
 475             elif asThreeFunctions is not None:
 476                 if isinstance(asThreeFunctions, dict) and \
 477                    ("Tangent" in asThreeFunctions) and (asThreeFunctions["Tangent"] is not None) and \
 478                    ("Adjoint" in asThreeFunctions) and (asThreeFunctions["Adjoint"] is not None) and \
 479                    (("useApproximatedDerivatives" not in asThreeFunctions) or not bool(asThreeFunctions["useApproximatedDerivatives"])):
 480                     __Function = asThreeFunctions
 481                 elif isinstance(asThreeFunctions, dict) and \
 482                    ("Direct" in asThreeFunctions) and (asThreeFunctions["Direct"] is not None):
 483                     __Function = asThreeFunctions
 484                     __Function.update({"useApproximatedDerivatives":True})
 485                 else:
 486                     raise ValueError("The functions has to be given in a dictionnary which have either 1 key (\"Direct\") or 3 keys (\"Direct\" (optionnal), \"Tangent\" and \"Adjoint\")")
 487                 if "Direct"  not in asThreeFunctions:
 488                     __Function["Direct"] = asThreeFunctions["Tangent"]
 489                 __Function.update(__Parameters)
 490             else:
 491                 __Function = None
 492         #
 493         # if sys.version_info[0] < 3 and isinstance(__Function, dict):
 494         #     for k in ("Direct", "Tangent", "Adjoint"):
 495         #         if k in __Function and hasattr(__Function[k],"__class__"):
 496         #             if type(__Function[k]) is type(self.__init__):
 497         #                 raise TypeError("can't use a class method (%s) as a function for the \"%s\" operator. Use a real function instead."%(type(__Function[k]),k))
 498         #
 499         if   appliedInX is not None and isinstance(appliedInX, dict):
 500             __appliedInX = appliedInX
 501         elif appliedInX is not None:
 502             __appliedInX = {"HXb":appliedInX}
 503         else:
 504             __appliedInX = None
 505         #
 506         if scheduledBy is not None:
 507             self.__T = scheduledBy
 508         #
 509         if isinstance(__Function, dict) and \
 510                 ("useApproximatedDerivatives" in __Function) and bool(__Function["useApproximatedDerivatives"]) and \
 511                 ("Direct" in __Function) and (__Function["Direct"] is not None):
 512             if "CenteredFiniteDifference"           not in __Function: __Function["CenteredFiniteDifference"]           = False
 513             if "DifferentialIncrement"              not in __Function: __Function["DifferentialIncrement"]              = 0.01
 514             if "withdX"                             not in __Function: __Function["withdX"]                             = None
 515             if "withAvoidingRedundancy"             not in __Function: __Function["withAvoidingRedundancy"]             = avoidRC
 516             if "withToleranceInRedundancy"          not in __Function: __Function["withToleranceInRedundancy"]          = 1.e-18
 517             if "withLenghtOfRedundancy"             not in __Function: __Function["withLenghtOfRedundancy"]             = -1
 518             if "NumberOfProcesses"                  not in __Function: __Function["NumberOfProcesses"]                  = None
 519             if "withmfEnabled"                      not in __Function: __Function["withmfEnabled"]                      = inputAsMF
 520             from daCore import NumericObjects
 521             FDA = NumericObjects.FDApproximation(
 522                 name                  = self.__name,
 523                 Function              = __Function["Direct"],
 524                 centeredDF            = __Function["CenteredFiniteDifference"],
 525                 increment             = __Function["DifferentialIncrement"],
 526                 dX                    = __Function["withdX"],
 527                 extraArguments        = self.__extraArgs,
 528                 avoidingRedundancy    = __Function["withAvoidingRedundancy"],
 529                 toleranceInRedundancy = __Function["withToleranceInRedundancy"],
 530                 lenghtOfRedundancy    = __Function["withLenghtOfRedundancy"],
 531                 mpEnabled             = __Function["EnableMultiProcessingInDerivatives"],
 532                 mpWorkers             = __Function["NumberOfProcesses"],
 533                 mfEnabled             = __Function["withmfEnabled"],
 534                 )
 535             self.__FO["Direct"]  = Operator( name = self.__name,           fromMethod = FDA.DirectOperator,  avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs, enableMultiProcess = __Parameters["EnableMultiProcessingInEvaluation"] )
 536             self.__FO["Tangent"] = Operator( name = self.__name+"Tangent", fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
 537             self.__FO["Adjoint"] = Operator( name = self.__name+"Adjoint", fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
 538         elif isinstance(__Function, dict) and \
 539                 ("Direct" in __Function) and ("Tangent" in __Function) and ("Adjoint" in __Function) and \
 540                 (__Function["Direct"] is not None) and (__Function["Tangent"] is not None) and (__Function["Adjoint"] is not None):
 541             self.__FO["Direct"]  = Operator( name = self.__name,           fromMethod = __Function["Direct"],  avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs, enableMultiProcess = __Parameters["EnableMultiProcessingInEvaluation"] )
 542             self.__FO["Tangent"] = Operator( name = self.__name+"Tangent", fromMethod = __Function["Tangent"], avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
 543             self.__FO["Adjoint"] = Operator( name = self.__name+"Adjoint", fromMethod = __Function["Adjoint"], avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
 544         elif asMatrix is not None:
 545             __matrice = numpy.matrix( __Matrix, numpy.float )
 546             self.__FO["Direct"]  = Operator( name = self.__name,           fromMatrix = __matrice,   avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, enableMultiProcess = __Parameters["EnableMultiProcessingInEvaluation"] )
 547             self.__FO["Tangent"] = Operator( name = self.__name+"Tangent", fromMatrix = __matrice,   avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF )
 548             self.__FO["Adjoint"] = Operator( name = self.__name+"Adjoint", fromMatrix = __matrice.T, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF )
 549             del __matrice
 550         else:
 551             raise ValueError("The %s object is improperly defined or undefined, it requires at minima either a matrix, a Direct operator for approximate derivatives or a Tangent/Adjoint operators pair. Please check your operator input."%self.__name)
 552         #
 553         if __appliedInX is not None:
 554             self.__FO["AppliedInX"] = {}
 555             for key in list(__appliedInX.keys()):
 556                 if type( __appliedInX[key] ) is type( numpy.matrix([]) ):
 557                     # Pour le cas où l'on a une vraie matrice
 558                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key].A1, numpy.float ).T
 559                 elif type( __appliedInX[key] ) is type( numpy.array([]) ) and len(__appliedInX[key].shape) > 1:
 560                     # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
 561                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key].reshape(len(__appliedInX[key]),), numpy.float ).T
 562                 else:
 563                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key],    numpy.float ).T
 564         else:
 565             self.__FO["AppliedInX"] = None
 566
 567     def getO(self):
 568         return self.__FO
 569
 570     def __repr__(self):
 571         "x.__repr__() <==> repr(x)"
 572         return repr(self.__FO)
 573
 574     def __str__(self):
 575         "x.__str__() <==> str(x)"
 576         return str(self.__FO)
 577
 578 # ==============================================================================
 579 class Algorithm(object):
 580     """
 581     Classe générale d'interface de type algorithme
 582
 583     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
 584     d'assimilation, en fournissant un container (dictionnaire) de variables
 585     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
 586
 587     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
 588     """
 589     def __init__(self, name):
 590         """
 591         L'initialisation présente permet de fabriquer des variables de stockage
 592         disponibles de manière générique dans les algorithmes élémentaires. Ces
 593         variables de stockage sont ensuite conservées dans un dictionnaire
 594         interne à l'objet, mais auquel on accède par la méthode "get".
 595
 596         Les variables prévues sont :
 597             - APosterioriCorrelations : matrice de corrélations de la matrice A
 598             - APosterioriCovariance : matrice de covariances a posteriori : A
 599             - APosterioriStandardDeviations : vecteur des écart-types de la matrice A
 600             - APosterioriVariances : vecteur des variances de la matrice A
 601             - Analysis : vecteur d'analyse : Xa
 602             - BMA : Background moins Analysis : Xa - Xb
 603             - CostFunctionJ  : fonction-coût globale, somme des deux parties suivantes Jb et Jo
 604             - CostFunctionJAtCurrentOptimum : fonction-coût globale à l'état optimal courant lors d'itérations
 605             - CostFunctionJb : partie ébauche ou background de la fonction-coût : Jb
 606             - CostFunctionJbAtCurrentOptimum : partie ébauche à l'état optimal courant lors d'itérations
 607             - CostFunctionJo : partie observations de la fonction-coût : Jo
 608             - CostFunctionJoAtCurrentOptimum : partie observations à l'état optimal courant lors d'itérations
 609             - CurrentIterationNumber : numéro courant d'itération dans les algorithmes itératifs, à partir de 0
 610             - CurrentOptimum : état optimal courant lors d'itérations
 611             - CurrentState : état courant lors d'itérations
 612             - GradientOfCostFunctionJ  : gradient de la fonction-coût globale
 613             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-coût
 614             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-coût
 615             - IndexOfOptimum : index de l'état optimal courant lors d'itérations
 616             - Innovation : l'innovation : d = Y - H(X)
 617             - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
 618             - JacobianMatrixAtBackground : matrice jacobienne à l'état d'ébauche
 619             - JacobianMatrixAtCurrentState : matrice jacobienne à l'état courant
 620             - JacobianMatrixAtOptimum : matrice jacobienne à l'optimum
 621             - KalmanGainAtOptimum : gain de Kalman à l'optimum
 622             - MahalanobisConsistency : indicateur de consistance des covariances
 623             - OMA : Observation moins Analyse : Y - Xa
 624             - OMB : Observation moins Background : Y - Xb
 625             - ForecastState : état prédit courant lors d'itérations
 626             - Residu : dans le cas des algorithmes de vérification
 627             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
 628             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
 629             - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
 630             - SimulatedObservationAtCurrentOptimum : l'état observé H(X) à l'état optimal courant
 631             - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
 632             - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
 633             - SimulationQuantiles : états observés H(X) pour les quantiles demandés
 634         On peut rajouter des variables à stocker dans l'initialisation de
 635         l'algorithme élémentaire qui va hériter de cette classe
 636         """
 637         logging.debug("%s Initialisation", str(name))
 638         self._m = PlatformInfo.SystemUsage()
 639         #
 640         self._name = str( name )
 641         self._parameters = {"StoreSupplementaryCalculations":[]}
 642         self.__required_parameters = {}
 643         self.__required_inputs = {
 644             "RequiredInputValues":{"mandatory":(), "optional":()},
 645             "ClassificationTags":[],
 646             }
 647         self.__variable_names_not_public = {"nextStep":False} # Duplication dans AlgorithmAndParameters
 648         self.__canonical_parameter_name = {} # Correspondance "lower"->"correct"
 649         self.__canonical_stored_name = {}    # Correspondance "lower"->"correct"
 650         #
 651         self.StoredVariables = {}
 652         self.StoredVariables["APosterioriCorrelations"]              = Persistence.OneMatrix(name = "APosterioriCorrelations")
 653         self.StoredVariables["APosterioriCovariance"]                = Persistence.OneMatrix(name = "APosterioriCovariance")
 654         self.StoredVariables["APosterioriStandardDeviations"]        = Persistence.OneVector(name = "APosterioriStandardDeviations")
 655         self.StoredVariables["APosterioriVariances"]                 = Persistence.OneVector(name = "APosterioriVariances")
 656         self.StoredVariables["Analysis"]                             = Persistence.OneVector(name = "Analysis")
 657         self.StoredVariables["BMA"]                                  = Persistence.OneVector(name = "BMA")
 658         self.StoredVariables["CostFunctionJ"]                        = Persistence.OneScalar(name = "CostFunctionJ")
 659         self.StoredVariables["CostFunctionJAtCurrentOptimum"]        = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
 660         self.StoredVariables["CostFunctionJb"]                       = Persistence.OneScalar(name = "CostFunctionJb")
 661         self.StoredVariables["CostFunctionJbAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
 662         self.StoredVariables["CostFunctionJo"]                       = Persistence.OneScalar(name = "CostFunctionJo")
 663         self.StoredVariables["CostFunctionJoAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
 664         self.StoredVariables["CurrentEnsembleState"]                 = Persistence.OneMatrix(name = "CurrentEnsembleState")
 665         self.StoredVariables["CurrentIterationNumber"]               = Persistence.OneIndex(name  = "CurrentIterationNumber")
 666         self.StoredVariables["CurrentOptimum"]                       = Persistence.OneVector(name = "CurrentOptimum")
 667         self.StoredVariables["CurrentState"]                         = Persistence.OneVector(name = "CurrentState")
 668         self.StoredVariables["ForecastState"]                        = Persistence.OneVector(name = "ForecastState")
 669         self.StoredVariables["GradientOfCostFunctionJ"]              = Persistence.OneVector(name = "GradientOfCostFunctionJ")
 670         self.StoredVariables["GradientOfCostFunctionJb"]             = Persistence.OneVector(name = "GradientOfCostFunctionJb")
 671         self.StoredVariables["GradientOfCostFunctionJo"]             = Persistence.OneVector(name = "GradientOfCostFunctionJo")
 672         self.StoredVariables["IndexOfOptimum"]                       = Persistence.OneIndex(name  = "IndexOfOptimum")
 673         self.StoredVariables["Innovation"]                           = Persistence.OneVector(name = "Innovation")
 674         self.StoredVariables["InnovationAtCurrentAnalysis"]          = Persistence.OneVector(name = "InnovationAtCurrentAnalysis")
 675         self.StoredVariables["InnovationAtCurrentState"]             = Persistence.OneVector(name = "InnovationAtCurrentState")
 676         self.StoredVariables["JacobianMatrixAtBackground"]           = Persistence.OneMatrix(name = "JacobianMatrixAtBackground")
 677         self.StoredVariables["JacobianMatrixAtCurrentState"]         = Persistence.OneMatrix(name = "JacobianMatrixAtCurrentState")
 678         self.StoredVariables["JacobianMatrixAtOptimum"]              = Persistence.OneMatrix(name = "JacobianMatrixAtOptimum")
 679         self.StoredVariables["KalmanGainAtOptimum"]                  = Persistence.OneMatrix(name = "KalmanGainAtOptimum")
 680         self.StoredVariables["MahalanobisConsistency"]               = Persistence.OneScalar(name = "MahalanobisConsistency")
 681         self.StoredVariables["OMA"]                                  = Persistence.OneVector(name = "OMA")
 682         self.StoredVariables["OMB"]                                  = Persistence.OneVector(name = "OMB")
 683         self.StoredVariables["Residu"]                               = Persistence.OneScalar(name = "Residu")
 684         self.StoredVariables["SigmaBck2"]                            = Persistence.OneScalar(name = "SigmaBck2")
 685         self.StoredVariables["SigmaObs2"]                            = Persistence.OneScalar(name = "SigmaObs2")
 686         self.StoredVariables["SimulatedObservationAtBackground"]     = Persistence.OneVector(name = "SimulatedObservationAtBackground")
 687         self.StoredVariables["SimulatedObservationAtCurrentAnalysis"]= Persistence.OneVector(name = "SimulatedObservationAtCurrentAnalysis")
 688         self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
 689         self.StoredVariables["SimulatedObservationAtCurrentState"]   = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
 690         self.StoredVariables["SimulatedObservationAtOptimum"]        = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
 691         self.StoredVariables["SimulationQuantiles"]                  = Persistence.OneMatrix(name = "SimulationQuantiles")
 692         #
 693         for k in self.StoredVariables:
 694             self.__canonical_stored_name[k.lower()] = k
 695         #
 696         for k, v in self.__variable_names_not_public.items():
 697             self.__canonical_parameter_name[k.lower()] = k
 698         self.__canonical_parameter_name["algorithm"] = "Algorithm"
 699         self.__canonical_parameter_name["storesupplementarycalculations"] = "StoreSupplementaryCalculations"
 700
 701     def _pre_run(self, Parameters, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None ):
 702         "Pré-calcul"
 703         logging.debug("%s Lancement", self._name)
 704         logging.debug("%s Taille mémoire utilisée de %.0f Mio"%(self._name, self._m.getUsedMemory("Mio")))
 705         self._getTimeState(reset=True)
 706         #
 707         # Mise a jour des paramètres internes avec le contenu de Parameters, en
 708         # reprenant les valeurs par défauts pour toutes celles non définies
 709         self.__setParameters(Parameters, reset=True)
 710         for k, v in self.__variable_names_not_public.items():
 711             if k not in self._parameters:  self.__setParameters( {k:v} )
 712         #
 713         # Corrections et compléments des vecteurs
 714         def __test_vvalue(argument, variable, argname, symbol=None):
 715             if symbol is None: symbol = variable
 716             if argument is None:
 717                 if variable in self.__required_inputs["RequiredInputValues"]["mandatory"]:
 718                     raise ValueError("%s %s vector %s is not set and has to be properly defined!"%(self._name,argname,symbol))
 719                 elif variable in self.__required_inputs["RequiredInputValues"]["optional"]:
 720                     logging.debug("%s %s vector %s is not set, but is optional."%(self._name,argname,symbol))
 721                 else:
 722                     logging.debug("%s %s vector %s is not set, but is not required."%(self._name,argname,symbol))
 723             else:
 724                 logging.debug("%s %s vector %s is set, and its size is %i."%(self._name,argname,symbol,numpy.array(argument).size))
 725             return 0
 726         __test_vvalue( Xb, "Xb", "Background or initial state" )
 727         __test_vvalue( Y,  "Y",  "Observation" )
 728         __test_vvalue( U,  "U",  "Control" )
 729         #
 730         # Corrections et compléments des covariances
 731         def __test_cvalue(argument, variable, argname, symbol=None):
 732             if symbol is None: symbol = variable
 733             if argument is None:
 734                 if variable in self.__required_inputs["RequiredInputValues"]["mandatory"]:
 735                     raise ValueError("%s %s error covariance matrix %s is not set and has to be properly defined!"%(self._name,argname,symbol))
 736                 elif variable in self.__required_inputs["RequiredInputValues"]["optional"]:
 737                     logging.debug("%s %s error covariance matrix %s is not set, but is optional."%(self._name,argname,symbol))
 738                 else:
 739                     logging.debug("%s %s error covariance matrix %s is not set, but is not required."%(self._name,argname,symbol))
 740             else:
 741                 logging.debug("%s %s error covariance matrix %s is set."%(self._name,argname,symbol))
 742             return 0
 743         __test_cvalue( B, "B", "Background" )
 744         __test_cvalue( R, "R", "Observation" )
 745         __test_cvalue( Q, "Q", "Evolution" )
 746         #
 747         # Corrections et compléments des opérateurs
 748         def __test_ovalue(argument, variable, argname, symbol=None):
 749             if symbol is None: symbol = variable
 750             if argument is None or (isinstance(argument,dict) and len(argument)==0):
 751                 if variable in self.__required_inputs["RequiredInputValues"]["mandatory"]:
 752                     raise ValueError("%s %s operator %s is not set and has to be properly defined!"%(self._name,argname,symbol))
 753                 elif variable in self.__required_inputs["RequiredInputValues"]["optional"]:
 754                     logging.debug("%s %s operator %s is not set, but is optional."%(self._name,argname,symbol))
 755                 else:
 756                     logging.debug("%s %s operator %s is not set, but is not required."%(self._name,argname,symbol))
 757             else:
 758                 logging.debug("%s %s operator %s is set."%(self._name,argname,symbol))
 759             return 0
 760         __test_ovalue( HO, "HO", "Observation", "H" )
 761         __test_ovalue( EM, "EM", "Evolution", "M" )
 762         __test_ovalue( CM, "CM", "Control Model", "C" )
 763         #
 764         # Corrections et compléments des bornes
 765         if ("Bounds" in self._parameters) and isinstance(self._parameters["Bounds"], (list, tuple)) and (len(self._parameters["Bounds"]) > 0):
 766             logging.debug("%s Prise en compte des bornes effectuee"%(self._name,))
 767         else:
 768             self._parameters["Bounds"] = None
 769         #
 770         # Corrections et compléments de l'initialisation en X
 771         if  "InitializationPoint" in self._parameters:
 772             if Xb is not None:
 773                 if self._parameters["InitializationPoint"] is not None and hasattr(self._parameters["InitializationPoint"],'size'):
 774                     if self._parameters["InitializationPoint"].size != numpy.ravel(Xb).size:
 775                         raise ValueError("Incompatible size %i of forced initial point that have to replace the background of size %i" \
 776                             %(self._parameters["InitializationPoint"].size,numpy.ravel(Xb).size))
 777                     # Obtenu par typecast : numpy.ravel(self._parameters["InitializationPoint"])
 778                 else:
 779                     self._parameters["InitializationPoint"] = numpy.ravel(Xb)
 780             else:
 781                 if self._parameters["InitializationPoint"] is None:
 782                     raise ValueError("Forced initial point can not be set without any given Background or required value")
 783         #
 784         # Correction pour pallier a un bug de TNC sur le retour du Minimum
 785         if "Minimizer" in self._parameters and self._parameters["Minimizer"] == "TNC":
 786             self.setParameterValue("StoreInternalVariables",True)
 787         #
 788         # Verbosité et logging
 789         if logging.getLogger().level < logging.WARNING:
 790             self._parameters["optiprint"], self._parameters["optdisp"] = 1, 1
 791             if PlatformInfo.has_scipy:
 792                 import scipy.optimize
 793                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_ALL
 794             else:
 795                 self._parameters["optmessages"] = 15
 796         else:
 797             self._parameters["optiprint"], self._parameters["optdisp"] = -1, 0
 798             if PlatformInfo.has_scipy:
 799                 import scipy.optimize
 800                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_NONE
 801             else:
 802                 self._parameters["optmessages"] = 15
 803         #
 804         return 0
 805
 806     def _post_run(self,_oH=None):
 807         "Post-calcul"
 808         if ("StoreSupplementaryCalculations" in self._parameters) and \
 809             "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
 810             for _A in self.StoredVariables["APosterioriCovariance"]:
 811                 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
 812                     self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
 813                 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
 814                     self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
 815                 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
 816                     _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
 817                     _C = numpy.dot(_EI, numpy.dot(_A, _EI))
 818                     self.StoredVariables["APosterioriCorrelations"].store( _C )
 819         if _oH is not None and "Direct" in _oH and "Tangent" in _oH and "Adjoint" in _oH:
 820             logging.debug("%s Nombre d'évaluation(s) de l'opérateur d'observation direct/tangent/adjoint.: %i/%i/%i", self._name, _oH["Direct"].nbcalls(0),_oH["Tangent"].nbcalls(0),_oH["Adjoint"].nbcalls(0))
 821             logging.debug("%s Nombre d'appels au cache d'opérateur d'observation direct/tangent/adjoint..: %i/%i/%i", self._name, _oH["Direct"].nbcalls(3),_oH["Tangent"].nbcalls(3),_oH["Adjoint"].nbcalls(3))
 822         logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
 823         logging.debug("%s Durées d'utilisation CPU de %.1fs et elapsed de %.1fs", self._name, self._getTimeState()[0], self._getTimeState()[1])
 824         logging.debug("%s Terminé", self._name)
 825         return 0
 826
 827     def _toStore(self, key):
 828         "True if in StoreSupplementaryCalculations, else False"
 829         return key in self._parameters["StoreSupplementaryCalculations"]
 830
 831     def get(self, key=None):
 832         """
 833         Renvoie l'une des variables stockées identifiée par la clé, ou le
 834         dictionnaire de l'ensemble des variables disponibles en l'absence de
 835         clé. Ce sont directement les variables sous forme objet qui sont
 836         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
 837         des classes de persistance.
 838         """
 839         if key is not None:
 840             return self.StoredVariables[self.__canonical_stored_name[key.lower()]]
 841         else:
 842             return self.StoredVariables
 843
 844     def __contains__(self, key=None):
 845         "D.__contains__(k) -> True if D has a key k, else False"
 846         if key is None or key.lower() not in self.__canonical_stored_name:
 847             return False
 848         else:
 849             return self.__canonical_stored_name[key.lower()] in self.StoredVariables
 850
 851     def keys(self):
 852         "D.keys() -> list of D's keys"
 853         if hasattr(self, "StoredVariables"):
 854             return self.StoredVariables.keys()
 855         else:
 856             return []
 857
 858     def pop(self, k, d):
 859         "D.pop(k[,d]) -> v, remove specified key and return the corresponding value"
 860         if hasattr(self, "StoredVariables") and k.lower() in self.__canonical_stored_name:
 861             return self.StoredVariables.pop(self.__canonical_stored_name[k.lower()], d)
 862         else:
 863             try:
 864                 msg = "'%s'"%k
 865             except:
 866                 raise TypeError("pop expected at least 1 arguments, got 0")
 867             "If key is not found, d is returned if given, otherwise KeyError is raised"
 868             try:
 869                 return d
 870             except:
 871                 raise KeyError(msg)
 872
 873     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
 874         """
 875         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
 876         sa forme mathématique la plus naturelle possible.
 877         """
 878         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
 879
 880     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None, listadv = None):
 881         """
 882         Permet de définir dans l'algorithme des paramètres requis et leurs
 883         caractéristiques par défaut.
 884         """
 885         if name is None:
 886             raise ValueError("A name is mandatory to define a required parameter.")
 887         #
 888         self.__required_parameters[name] = {
 889             "default"  : default,
 890             "typecast" : typecast,
 891             "minval"   : minval,
 892             "maxval"   : maxval,
 893             "listval"  : listval,
 894             "listadv"  : listadv,
 895             "message"  : message,
 896             }
 897         self.__canonical_parameter_name[name.lower()] = name
 898         logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
 899
 900     def getRequiredParameters(self, noDetails=True):
 901         """
 902         Renvoie la liste des noms de paramètres requis ou directement le
 903         dictionnaire des paramètres requis.
 904         """
 905         if noDetails:
 906             return sorted(self.__required_parameters.keys())
 907         else:
 908             return self.__required_parameters
 909
 910     def setParameterValue(self, name=None, value=None):
 911         """
 912         Renvoie la valeur d'un paramètre requis de manière contrôlée
 913         """
 914         __k = self.__canonical_parameter_name[name.lower()]
 915         default  = self.__required_parameters[__k]["default"]
 916         typecast = self.__required_parameters[__k]["typecast"]
 917         minval   = self.__required_parameters[__k]["minval"]
 918         maxval   = self.__required_parameters[__k]["maxval"]
 919         listval  = self.__required_parameters[__k]["listval"]
 920         listadv  = self.__required_parameters[__k]["listadv"]
 921         #
 922         if value is None and default is None:
 923             __val = None
 924         elif value is None and default is not None:
 925             if typecast is None: __val = default
 926             else:                __val = typecast( default )
 927         else:
 928             if typecast is None: __val = value
 929             else:
 930                 try:
 931                     __val = typecast( value )
 932                 except:
 933                     raise ValueError("The value '%s' for the parameter named '%s' can not be correctly evaluated with type '%s'."%(value, __k, typecast))
 934         #
 935         if minval is not None and (numpy.array(__val, float) < minval).any():
 936             raise ValueError("The parameter named '%s' of value '%s' can not be less than %s."%(__k, __val, minval))
 937         if maxval is not None and (numpy.array(__val, float) > maxval).any():
 938             raise ValueError("The parameter named '%s' of value '%s' can not be greater than %s."%(__k, __val, maxval))
 939         if listval is not None or listadv is not None:
 940             if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
 941                 for v in __val:
 942                     if listval is not None and v in listval: continue
 943                     elif listadv is not None and v in listadv: continue
 944                     else:
 945                         raise ValueError("The value '%s' is not allowed for the parameter named '%s', it has to be in the list %s."%(v, __k, listval))
 946             elif not (listval is not None and __val in listval) and not (listadv is not None and __val in listadv):
 947                 raise ValueError("The value '%s' is not allowed for the parameter named '%s', it has to be in the list %s."%( __val, __k,listval))
 948         #
 949         return __val
 950
 951     def requireInputArguments(self, mandatory=(), optional=()):
 952         """
 953         Permet d'imposer des arguments de calcul requis en entrée.
 954         """
 955         self.__required_inputs["RequiredInputValues"]["mandatory"] = tuple( mandatory )
 956         self.__required_inputs["RequiredInputValues"]["optional"]  = tuple( optional )
 957
 958     def getInputArguments(self):
 959         """
 960         Permet d'obtenir les listes des arguments de calcul requis en entrée.
 961         """
 962         return self.__required_inputs["RequiredInputValues"]["mandatory"], self.__required_inputs["RequiredInputValues"]["optional"]
 963
 964     def setAttributes(self, tags=()):
 965         """
 966         Permet d'adjoindre des attributs comme les tags de classification.
 967         Renvoie la liste actuelle dans tous les cas.
 968         """
 969         self.__required_inputs["ClassificationTags"].extend( tags )
 970         return self.__required_inputs["ClassificationTags"]
 971
 972     def __setParameters(self, fromDico={}, reset=False):
 973         """
 974         Permet de stocker les paramètres reçus dans le dictionnaire interne.
 975         """
 976         self._parameters.update( fromDico )
 977         __inverse_fromDico_keys = {}
 978         for k in fromDico.keys():
 979             if k.lower() in self.__canonical_parameter_name:
 980                 __inverse_fromDico_keys[self.__canonical_parameter_name[k.lower()]] = k
 981         #~ __inverse_fromDico_keys = dict([(self.__canonical_parameter_name[k.lower()],k) for k in fromDico.keys()])
 982         __canonic_fromDico_keys = __inverse_fromDico_keys.keys()
 983         for k in self.__required_parameters.keys():
 984             if k in __canonic_fromDico_keys:
 985                 self._parameters[k] = self.setParameterValue(k,fromDico[__inverse_fromDico_keys[k]])
 986             elif reset:
 987                 self._parameters[k] = self.setParameterValue(k)
 988             else:
 989                 pass
 990             logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
 991
 992     def _getTimeState(self, reset=False):
 993         """
 994         Initialise ou restitue le temps de calcul (cpu/elapsed) à la seconde
 995         """
 996         if reset:
 997             self.__initial_cpu_time      = time.process_time()
 998             self.__initial_elapsed_time  = time.perf_counter()
 999             return 0., 0.
1000         else:
1001             self.__cpu_time     = time.process_time() - self.__initial_cpu_time
1002             self.__elapsed_time = time.perf_counter() - self.__initial_elapsed_time
1003             return self.__cpu_time, self.__elapsed_time
1004
1005     def _StopOnTimeLimit(self, X=None, withReason=False):
1006         "Stop criteria on time limit: True/False [+ Reason]"
1007         c, e = self._getTimeState()
1008         if "MaximumCpuTime" in self._parameters and c > self._parameters["MaximumCpuTime"]:
1009             __SC, __SR = True, "Reached maximum CPU time (%.1fs > %.1fs)"%(c, self._parameters["MaximumCpuTime"])
1010         elif "MaximumElapsedTime" in self._parameters and e > self._parameters["MaximumElapsedTime"]:
1011             __SC, __SR = True, "Reached maximum elapsed time (%.1fs > %.1fs)"%(e, self._parameters["MaximumElapsedTime"])
1012         else:
1013             __SC, __SR = False, ""
1014         if withReason:
1015             return __SC, __SR
1016         else:
1017             return __SC
1018
1019 # ==============================================================================
1020 class AlgorithmAndParameters(object):
1021     """
1022     Classe générale d'interface d'action pour l'algorithme et ses paramètres
1023     """
1024     def __init__(self,
1025                  name               = "GenericAlgorithm",
1026                  asAlgorithm        = None,
1027                  asDict             = None,
1028                  asScript           = None,
1029                 ):
1030         """
1031         """
1032         self.__name       = str(name)
1033         self.__A          = None
1034         self.__P          = {}
1035         #
1036         self.__algorithm         = {}
1037         self.__algorithmFile     = None
1038         self.__algorithmName     = None
1039         #
1040         self.updateParameters( asDict, asScript )
1041         #
1042         if asAlgorithm is None and asScript is not None:
1043             __Algo = Interfaces.ImportFromScript(asScript).getvalue( "Algorithm" )
1044         else:
1045             __Algo = asAlgorithm
1046         #
1047         if __Algo is not None:
1048             self.__A = str(__Algo)
1049             self.__P.update( {"Algorithm":self.__A} )
1050         #
1051         self.__setAlgorithm( self.__A )
1052         #
1053         self.__variable_names_not_public = {"nextStep":False} # Duplication dans Algorithm
1054
1055     def updateParameters(self,
1056                  asDict     = None,
1057                  asScript   = None,
1058                 ):
1059         "Mise a jour des parametres"
1060         if asDict is None and asScript is not None:
1061             __Dict = Interfaces.ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
1062         else:
1063             __Dict = asDict
1064         #
1065         if __Dict is not None:
1066             self.__P.update( dict(__Dict) )
1067
1068     def executePythonScheme(self, asDictAO = None):
1069         "Permet de lancer le calcul d'assimilation"
1070         Operator.CM.clearCache()
1071         #
1072         if not isinstance(asDictAO, dict):
1073             raise ValueError("The objects for algorithm calculation have to be given together as a dictionnary, and they are not")
1074         if   hasattr(asDictAO["Background"],"getO"):        self.__Xb = asDictAO["Background"].getO()
1075         elif hasattr(asDictAO["CheckingPoint"],"getO"):     self.__Xb = asDictAO["CheckingPoint"].getO()
1076         else:                                               self.__Xb = None
1077         if hasattr(asDictAO["Observation"],"getO"):         self.__Y  = asDictAO["Observation"].getO()
1078         else:                                               self.__Y  = asDictAO["Observation"]
1079         if hasattr(asDictAO["ControlInput"],"getO"):        self.__U  = asDictAO["ControlInput"].getO()
1080         else:                                               self.__U  = asDictAO["ControlInput"]
1081         if hasattr(asDictAO["ObservationOperator"],"getO"): self.__HO = asDictAO["ObservationOperator"].getO()
1082         else:                                               self.__HO = asDictAO["ObservationOperator"]
1083         if hasattr(asDictAO["EvolutionModel"],"getO"):      self.__EM = asDictAO["EvolutionModel"].getO()
1084         else:                                               self.__EM = asDictAO["EvolutionModel"]
1085         if hasattr(asDictAO["ControlModel"],"getO"):        self.__CM = asDictAO["ControlModel"].getO()
1086         else:                                               self.__CM = asDictAO["ControlModel"]
1087         self.__B = asDictAO["BackgroundError"]
1088         self.__R = asDictAO["ObservationError"]
1089         self.__Q = asDictAO["EvolutionError"]
1090         #
1091         self.__shape_validate()
1092         #
1093         self.__algorithm.run(
1094             Xb         = self.__Xb,
1095             Y          = self.__Y,
1096             U          = self.__U,
1097             HO         = self.__HO,
1098             EM         = self.__EM,
1099             CM         = self.__CM,
1100             R          = self.__R,
1101             B          = self.__B,
1102             Q          = self.__Q,
1103             Parameters = self.__P,
1104             )
1105         return 0
1106
1107     def executeYACSScheme(self, FileName=None):
1108         "Permet de lancer le calcul d'assimilation"
1109         if FileName is None or not os.path.exists(FileName):
1110             raise ValueError("a YACS file name has to be given for YACS execution.\n")
1111         else:
1112             __file    = os.path.abspath(FileName)
1113             logging.debug("The YACS file name is \"%s\"."%__file)
1114         if not PlatformInfo.has_salome or \
1115             not PlatformInfo.has_yacs or \
1116             not PlatformInfo.has_adao:
1117             raise ImportError("\n\n"+\
1118                 "Unable to get SALOME, YACS or ADAO environnement variables.\n"+\
1119                 "Please load the right environnement before trying to use it.\n")
1120         #
1121         import pilot
1122         import SALOMERuntime
1123         import loader
1124         SALOMERuntime.RuntimeSALOME_setRuntime()
1125
1126         r = pilot.getRuntime()
1127         xmlLoader = loader.YACSLoader()
1128         xmlLoader.registerProcCataLoader()
1129         try:
1130             catalogAd = r.loadCatalog("proc", __file)
1131             r.addCatalog(catalogAd)
1132         except:
1133             pass
1134
1135         try:
1136             p = xmlLoader.load(__file)
1137         except IOError as ex:
1138             print("The YACS XML schema file can not be loaded: %s"%(ex,))
1139
1140         logger = p.getLogger("parser")
1141         if not logger.isEmpty():
1142             print("The imported YACS XML schema has errors on parsing:")
1143             print(logger.getStr())
1144
1145         if not p.isValid():
1146             print("The YACS XML schema is not valid and will not be executed:")
1147             print(p.getErrorReport())
1148
1149         info=pilot.LinkInfo(pilot.LinkInfo.ALL_DONT_STOP)
1150         p.checkConsistency(info)
1151         if info.areWarningsOrErrors():
1152             print("The YACS XML schema is not coherent and will not be executed:")
1153             print(info.getGlobalRepr())
1154
1155         e = pilot.ExecutorSwig()
1156         e.RunW(p)
1157         if p.getEffectiveState() != pilot.DONE:
1158             print(p.getErrorReport())
1159         #
1160         return 0
1161
1162     def get(self, key = None):
1163         "Vérifie l'existence d'une clé de variable ou de paramètres"
1164         if key in self.__algorithm:
1165             return self.__algorithm.get( key )
1166         elif key in self.__P:
1167             return self.__P[key]
1168         else:
1169             allvariables = self.__P
1170             for k in self.__variable_names_not_public: allvariables.pop(k, None)
1171             return allvariables
1172
1173     def pop(self, k, d):
1174         "Necessaire pour le pickling"
1175         return self.__algorithm.pop(k, d)
1176
1177     def getAlgorithmRequiredParameters(self, noDetails=True):
1178         "Renvoie la liste des paramètres requis selon l'algorithme"
1179         return self.__algorithm.getRequiredParameters(noDetails)
1180
1181     def getAlgorithmInputArguments(self):
1182         "Renvoie la liste des entrées requises selon l'algorithme"
1183         return self.__algorithm.getInputArguments()
1184
1185     def getAlgorithmAttributes(self):
1186         "Renvoie la liste des attributs selon l'algorithme"
1187         return self.__algorithm.setAttributes()
1188
1189     def setObserver(self, __V, __O, __I, __S):
1190         if self.__algorithm is None \
1191             or isinstance(self.__algorithm, dict) \
1192             or not hasattr(self.__algorithm,"StoredVariables"):
1193             raise ValueError("No observer can be build before choosing an algorithm.")
1194         if __V not in self.__algorithm:
1195             raise ValueError("An observer requires to be set on a variable named %s which does not exist."%__V)
1196         else:
1197             self.__algorithm.StoredVariables[ __V ].setDataObserver(
1198                     Scheduler      = __S,
1199                     HookFunction   = __O,
1200                     HookParameters = __I,
1201                     )
1202
1203     def removeObserver(self, __V, __O, __A = False):
1204         if self.__algorithm is None \
1205             or isinstance(self.__algorithm, dict) \
1206             or not hasattr(self.__algorithm,"StoredVariables"):
1207             raise ValueError("No observer can be removed before choosing an algorithm.")
1208         if __V not in self.__algorithm:
1209             raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%__V)
1210         else:
1211             return self.__algorithm.StoredVariables[ __V ].removeDataObserver(
1212                     HookFunction   = __O,
1213                     AllObservers   = __A,
1214                     )
1215
1216     def hasObserver(self, __V):
1217         if self.__algorithm is None \
1218             or isinstance(self.__algorithm, dict) \
1219             or not hasattr(self.__algorithm,"StoredVariables"):
1220             return False
1221         if __V not in self.__algorithm:
1222             return False
1223         return self.__algorithm.StoredVariables[ __V ].hasDataObserver()
1224
1225     def keys(self):
1226         __allvariables = list(self.__algorithm.keys()) + list(self.__P.keys())
1227         for k in self.__variable_names_not_public:
1228             if k in __allvariables: __allvariables.remove(k)
1229         return __allvariables
1230
1231     def __contains__(self, key=None):
1232         "D.__contains__(k) -> True if D has a key k, else False"
1233         return key in self.__algorithm or key in self.__P
1234
1235     def __repr__(self):
1236         "x.__repr__() <==> repr(x)"
1237         return repr(self.__A)+", "+repr(self.__P)
1238
1239     def __str__(self):
1240         "x.__str__() <==> str(x)"
1241         return str(self.__A)+", "+str(self.__P)
1242
1243     def __setAlgorithm(self, choice = None ):
1244         """
1245         Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
1246         d'assimilation. L'argument est un champ caractère se rapportant au nom
1247         d'un algorithme réalisant l'opération sur les arguments fixes.
1248         """
1249         if choice is None:
1250             raise ValueError("Error: algorithm choice has to be given")
1251         if self.__algorithmName is not None:
1252             raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
1253         daDirectory = "daAlgorithms"
1254         #
1255         # Recherche explicitement le fichier complet
1256         # ------------------------------------------
1257         module_path = None
1258         for directory in sys.path:
1259             if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
1260                 module_path = os.path.abspath(os.path.join(directory, daDirectory))
1261         if module_path is None:
1262             raise ImportError("No algorithm module named \"%s\" has been found in the search path.\n             The search path is %s"%(choice, sys.path))
1263         #
1264         # Importe le fichier complet comme un module
1265         # ------------------------------------------
1266         try:
1267             sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
1268             self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
1269             if not hasattr(self.__algorithmFile, "ElementaryAlgorithm"):
1270                 raise ImportError("this module does not define a valid elementary algorithm.")
1271             self.__algorithmName = str(choice)
1272             sys.path = sys_path_tmp ; del sys_path_tmp
1273         except ImportError as e:
1274             raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n             The import error message is: %s"%(choice,e))
1275         #
1276         # Instancie un objet du type élémentaire du fichier
1277         # -------------------------------------------------
1278         self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
1279         return 0
1280
1281     def __shape_validate(self):
1282         """
1283         Validation de la correspondance correcte des tailles des variables et
1284         des matrices s'il y en a.
1285         """
1286         if self.__Xb is None:                      __Xb_shape = (0,)
1287         elif hasattr(self.__Xb,"size"):            __Xb_shape = (self.__Xb.size,)
1288         elif hasattr(self.__Xb,"shape"):
1289             if isinstance(self.__Xb.shape, tuple): __Xb_shape = self.__Xb.shape
1290             else:                                  __Xb_shape = self.__Xb.shape()
1291         else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
1292         #
1293         if self.__Y is None:                       __Y_shape = (0,)
1294         elif hasattr(self.__Y,"size"):             __Y_shape = (self.__Y.size,)
1295         elif hasattr(self.__Y,"shape"):
1296             if isinstance(self.__Y.shape, tuple):  __Y_shape = self.__Y.shape
1297             else:                                  __Y_shape = self.__Y.shape()
1298         else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
1299         #
1300         if self.__U is None:                       __U_shape = (0,)
1301         elif hasattr(self.__U,"size"):             __U_shape = (self.__U.size,)
1302         elif hasattr(self.__U,"shape"):
1303             if isinstance(self.__U.shape, tuple):  __U_shape = self.__U.shape
1304             else:                                  __U_shape = self.__U.shape()
1305         else: raise TypeError("The control (U) has no attribute of shape: problem !")
1306         #
1307         if self.__B is None:                       __B_shape = (0,0)
1308         elif hasattr(self.__B,"shape"):
1309             if isinstance(self.__B.shape, tuple):  __B_shape = self.__B.shape
1310             else:                                  __B_shape = self.__B.shape()
1311         else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
1312         #
1313         if self.__R is None:                       __R_shape = (0,0)
1314         elif hasattr(self.__R,"shape"):
1315             if isinstance(self.__R.shape, tuple):  __R_shape = self.__R.shape
1316             else:                                  __R_shape = self.__R.shape()
1317         else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
1318         #
1319         if self.__Q is None:                       __Q_shape = (0,0)
1320         elif hasattr(self.__Q,"shape"):
1321             if isinstance(self.__Q.shape, tuple):  __Q_shape = self.__Q.shape
1322             else:                                  __Q_shape = self.__Q.shape()
1323         else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
1324         #
1325         if len(self.__HO) == 0:                              __HO_shape = (0,0)
1326         elif isinstance(self.__HO, dict):                    __HO_shape = (0,0)
1327         elif hasattr(self.__HO["Direct"],"shape"):
1328             if isinstance(self.__HO["Direct"].shape, tuple): __HO_shape = self.__HO["Direct"].shape
1329             else:                                            __HO_shape = self.__HO["Direct"].shape()
1330         else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
1331         #
1332         if len(self.__EM) == 0:                              __EM_shape = (0,0)
1333         elif isinstance(self.__EM, dict):                    __EM_shape = (0,0)
1334         elif hasattr(self.__EM["Direct"],"shape"):
1335             if isinstance(self.__EM["Direct"].shape, tuple): __EM_shape = self.__EM["Direct"].shape
1336             else:                                            __EM_shape = self.__EM["Direct"].shape()
1337         else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
1338         #
1339         if len(self.__CM) == 0:                              __CM_shape = (0,0)
1340         elif isinstance(self.__CM, dict):                    __CM_shape = (0,0)
1341         elif hasattr(self.__CM["Direct"],"shape"):
1342             if isinstance(self.__CM["Direct"].shape, tuple): __CM_shape = self.__CM["Direct"].shape
1343             else:                                            __CM_shape = self.__CM["Direct"].shape()
1344         else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
1345         #
1346         # Vérification des conditions
1347         # ---------------------------
1348         if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
1349             raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
1350         if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
1351             raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
1352         #
1353         if not( min(__B_shape) == max(__B_shape) ):
1354             raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
1355         if not( min(__R_shape) == max(__R_shape) ):
1356             raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
1357         if not( min(__Q_shape) == max(__Q_shape) ):
1358             raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
1359         if not( min(__EM_shape) == max(__EM_shape) ):
1360             raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
1361         #
1362         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and not( __HO_shape[1] == max(__Xb_shape) ):
1363             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
1364         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and not( __HO_shape[0] == max(__Y_shape) ):
1365             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
1366         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and len(self.__B) > 0 and not( __HO_shape[1] == __B_shape[0] ):
1367             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
1368         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and len(self.__R) > 0 and not( __HO_shape[0] == __R_shape[1] ):
1369             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
1370         #
1371         if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
1372             if self.__algorithmName in ["EnsembleBlue",]:
1373                 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
1374                 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
1375                 for member in asPersistentVector:
1376                     self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
1377                 __Xb_shape = min(__B_shape)
1378             else:
1379                 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
1380         #
1381         if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
1382             raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
1383         #
1384         if self.__EM is not None and len(self.__EM) > 0 and not isinstance(self.__EM, dict) and not( __EM_shape[1] == max(__Xb_shape) ):
1385             raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
1386         #
1387         if self.__CM is not None and len(self.__CM) > 0 and not isinstance(self.__CM, dict) and not( __CM_shape[1] == max(__U_shape) ):
1388             raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
1389         #
1390         if ("Bounds" in self.__P) \
1391             and (isinstance(self.__P["Bounds"], list) or isinstance(self.__P["Bounds"], tuple)) \
1392             and (len(self.__P["Bounds"]) != max(__Xb_shape)):
1393             raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
1394                 %(len(self.__P["Bounds"]),max(__Xb_shape)))
1395         #
1396         return 1
1397
1398 # ==============================================================================
1399 class RegulationAndParameters(object):
1400     """
1401     Classe générale d'interface d'action pour la régulation et ses paramètres
1402     """
1403     def __init__(self,
1404                  name               = "GenericRegulation",
1405                  asAlgorithm        = None,
1406                  asDict             = None,
1407                  asScript           = None,
1408                 ):
1409         """
1410         """
1411         self.__name       = str(name)
1412         self.__P          = {}
1413         #
1414         if asAlgorithm is None and asScript is not None:
1415             __Algo = Interfaces.ImportFromScript(asScript).getvalue( "Algorithm" )
1416         else:
1417             __Algo = asAlgorithm
1418         #
1419         if asDict is None and asScript is not None:
1420             __Dict = Interfaces.ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
1421         else:
1422             __Dict = asDict
1423         #
1424         if __Dict is not None:
1425             self.__P.update( dict(__Dict) )
1426         #
1427         if __Algo is not None:
1428             self.__P.update( {"Algorithm":str(__Algo)} )
1429
1430     def get(self, key = None):
1431         "Vérifie l'existence d'une clé de variable ou de paramètres"
1432         if key in self.__P:
1433             return self.__P[key]
1434         else:
1435             return self.__P
1436
1437 # ==============================================================================
1438 class DataObserver(object):
1439     """
1440     Classe générale d'interface de type observer
1441     """
1442     def __init__(self,
1443                  name        = "GenericObserver",
1444                  onVariable  = None,
1445                  asTemplate  = None,
1446                  asString    = None,
1447                  asScript    = None,
1448                  asObsObject = None,
1449                  withInfo    = None,
1450                  scheduledBy = None,
1451                  withAlgo    = None,
1452                 ):
1453         """
1454         """
1455         self.__name       = str(name)
1456         self.__V          = None
1457         self.__O          = None
1458         self.__I          = None
1459         #
1460         if onVariable is None:
1461             raise ValueError("setting an observer has to be done over a variable name or a list of variable names, not over None.")
1462         elif type(onVariable) in (tuple, list):
1463             self.__V = tuple(map( str, onVariable ))
1464             if withInfo is None:
1465                 self.__I = self.__V
1466             else:
1467                 self.__I = (str(withInfo),)*len(self.__V)
1468         elif isinstance(onVariable, str):
1469             self.__V = (onVariable,)
1470             if withInfo is None:
1471                 self.__I = (onVariable,)
1472             else:
1473                 self.__I = (str(withInfo),)
1474         else:
1475             raise ValueError("setting an observer has to be done over a variable name or a list of variable names.")
1476         #
1477         if asObsObject is not None:
1478             self.__O = asObsObject
1479         else:
1480             __FunctionText = str(UserScript('Observer', asTemplate, asString, asScript))
1481             __Function = Observer2Func(__FunctionText)
1482             self.__O = __Function.getfunc()
1483         #
1484         for k in range(len(self.__V)):
1485             ename = self.__V[k]
1486             einfo = self.__I[k]
1487             if ename not in withAlgo:
1488                 raise ValueError("An observer is asked to be set on a variable named %s which does not exist."%ename)
1489             else:
1490                 withAlgo.setObserver(ename, self.__O, einfo, scheduledBy)
1491
1492     def __repr__(self):
1493         "x.__repr__() <==> repr(x)"
1494         return repr(self.__V)+"\n"+repr(self.__O)
1495
1496     def __str__(self):
1497         "x.__str__() <==> str(x)"
1498         return str(self.__V)+"\n"+str(self.__O)
1499
1500 # ==============================================================================
1501 class UserScript(object):
1502     """
1503     Classe générale d'interface de type texte de script utilisateur
1504     """
1505     def __init__(self,
1506                  name       = "GenericUserScript",
1507                  asTemplate = None,
1508                  asString   = None,
1509                  asScript   = None,
1510                 ):
1511         """
1512         """
1513         self.__name       = str(name)
1514         #
1515         if asString is not None:
1516             self.__F = asString
1517         elif self.__name == "UserPostAnalysis" and (asTemplate is not None) and (asTemplate in Templates.UserPostAnalysisTemplates):
1518             self.__F = Templates.UserPostAnalysisTemplates[asTemplate]
1519         elif self.__name == "Observer" and (asTemplate is not None) and (asTemplate in Templates.ObserverTemplates):
1520             self.__F = Templates.ObserverTemplates[asTemplate]
1521         elif asScript is not None:
1522             self.__F = Interfaces.ImportFromScript(asScript).getstring()
1523         else:
1524             self.__F = ""
1525
1526     def __repr__(self):
1527         "x.__repr__() <==> repr(x)"
1528         return repr(self.__F)
1529
1530     def __str__(self):
1531         "x.__str__() <==> str(x)"
1532         return str(self.__F)
1533
1534 # ==============================================================================
1535 class ExternalParameters(object):
1536     """
1537     Classe générale d'interface de type texte de script utilisateur
1538     """
1539     def __init__(self,
1540                  name        = "GenericExternalParameters",
1541                  asDict      = None,
1542                  asScript    = None,
1543                 ):
1544         """
1545         """
1546         self.__name = str(name)
1547         self.__P    = {}
1548         #
1549         self.updateParameters( asDict, asScript )
1550
1551     def updateParameters(self,
1552                  asDict     = None,
1553                  asScript   = None,
1554                 ):
1555         "Mise a jour des parametres"
1556         if asDict is None and asScript is not None:
1557             __Dict = Interfaces.ImportFromScript(asScript).getvalue( self.__name, "ExternalParameters" )
1558         else:
1559             __Dict = asDict
1560         #
1561         if __Dict is not None:
1562             self.__P.update( dict(__Dict) )
1563
1564     def get(self, key = None):
1565         if key in self.__P:
1566             return self.__P[key]
1567         else:
1568             return list(self.__P.keys())
1569
1570     def keys(self):
1571         return list(self.__P.keys())
1572
1573     def pop(self, k, d):
1574         return self.__P.pop(k, d)
1575
1576     def items(self):
1577         return self.__P.items()
1578
1579     def __contains__(self, key=None):
1580         "D.__contains__(k) -> True if D has a key k, else False"
1581         return key in self.__P
1582
1583 # ==============================================================================
1584 class State(object):
1585     """
1586     Classe générale d'interface de type état
1587     """
1588     def __init__(self,
1589                  name               = "GenericVector",
1590                  asVector           = None,
1591                  asPersistentVector = None,
1592                  asScript           = None,
1593                  asDataFile         = None,
1594                  colNames           = None,
1595                  colMajor           = False,
1596                  scheduledBy        = None,
1597                  toBeChecked        = False,
1598                 ):
1599         """
1600         Permet de définir un vecteur :
1601         - asVector : entrée des données, comme un vecteur compatible avec le
1602           constructeur de numpy.matrix, ou "True" si entrée par script.
1603         - asPersistentVector : entrée des données, comme une série de vecteurs
1604           compatible avec le constructeur de numpy.matrix, ou comme un objet de
1605           type Persistence, ou "True" si entrée par script.
1606         - asScript : si un script valide est donné contenant une variable
1607           nommée "name", la variable est de type "asVector" (par défaut) ou
1608           "asPersistentVector" selon que l'une de ces variables est placée à
1609           "True".
1610         - asDataFile : si un ou plusieurs fichiers valides sont donnés
1611           contenant des valeurs en colonnes, elles-mêmes nommées "colNames"
1612           (s'il n'y a pas de nom de colonne indiquée, on cherche une colonne
1613           nommée "name"), on récupère les colonnes et on les range ligne après
1614           ligne (colMajor=False, par défaut) ou colonne après colonne
1615           (colMajor=True). La variable résultante est de type "asVector" (par
1616           défaut) ou "asPersistentVector" selon que l'une de ces variables est
1617           placée à "True".
1618         """
1619         self.__name       = str(name)
1620         self.__check      = bool(toBeChecked)
1621         #
1622         self.__V          = None
1623         self.__T          = None
1624         self.__is_vector  = False
1625         self.__is_series  = False
1626         #
1627         if asScript is not None:
1628             __Vector, __Series = None, None
1629             if asPersistentVector:
1630                 __Series = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1631             else:
1632                 __Vector = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1633         elif asDataFile is not None:
1634             __Vector, __Series = None, None
1635             if asPersistentVector:
1636                 if colNames is not None:
1637                     __Series = Interfaces.ImportFromFile(asDataFile).getvalue( colNames )[1]
1638                 else:
1639                     __Series = Interfaces.ImportFromFile(asDataFile).getvalue( [self.__name,] )[1]
1640                 if bool(colMajor) and not Interfaces.ImportFromFile(asDataFile).getformat() == "application/numpy.npz":
1641                     __Series = numpy.transpose(__Series)
1642                 elif not bool(colMajor) and Interfaces.ImportFromFile(asDataFile).getformat() == "application/numpy.npz":
1643                     __Series = numpy.transpose(__Series)
1644             else:
1645                 if colNames is not None:
1646                     __Vector = Interfaces.ImportFromFile(asDataFile).getvalue( colNames )[1]
1647                 else:
1648                     __Vector = Interfaces.ImportFromFile(asDataFile).getvalue( [self.__name,] )[1]
1649                 if bool(colMajor):
1650                     __Vector = numpy.ravel(__Vector, order = "F")
1651                 else:
1652                     __Vector = numpy.ravel(__Vector, order = "C")
1653         else:
1654             __Vector, __Series = asVector, asPersistentVector
1655         #
1656         if __Vector is not None:
1657             self.__is_vector = True
1658             self.__V         = numpy.matrix( numpy.asmatrix(__Vector).A1, numpy.float ).T
1659             self.shape       = self.__V.shape
1660             self.size        = self.__V.size
1661         elif __Series is not None:
1662             self.__is_series  = True
1663             if isinstance(__Series, (tuple, list, numpy.ndarray, numpy.matrix, str)):
1664                 self.__V = Persistence.OneVector(self.__name, basetype=numpy.matrix)
1665                 if isinstance(__Series, str): __Series = eval(__Series)
1666                 for member in __Series:
1667                     self.__V.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
1668             else:
1669                 self.__V = __Series
1670             if isinstance(self.__V.shape, (tuple, list)):
1671                 self.shape       = self.__V.shape
1672             else:
1673                 self.shape       = self.__V.shape()
1674             if len(self.shape) == 1:
1675                 self.shape       = (self.shape[0],1)
1676             self.size        = self.shape[0] * self.shape[1]
1677         else:
1678             raise ValueError("The %s object is improperly defined or undefined, it requires at minima either a vector, a list/tuple of vectors or a persistent object. Please check your vector input."%self.__name)
1679         #
1680         if scheduledBy is not None:
1681             self.__T = scheduledBy
1682
1683     def getO(self, withScheduler=False):
1684         if withScheduler:
1685             return self.__V, self.__T
1686         elif self.__T is None:
1687             return self.__V
1688         else:
1689             return self.__V
1690
1691     def isvector(self):
1692         "Vérification du type interne"
1693         return self.__is_vector
1694
1695     def isseries(self):
1696         "Vérification du type interne"
1697         return self.__is_series
1698
1699     def __repr__(self):
1700         "x.__repr__() <==> repr(x)"
1701         return repr(self.__V)
1702
1703     def __str__(self):
1704         "x.__str__() <==> str(x)"
1705         return str(self.__V)
1706
1707 # ==============================================================================
1708 class Covariance(object):
1709     """
1710     Classe générale d'interface de type covariance
1711     """
1712     def __init__(self,
1713                  name          = "GenericCovariance",
1714                  asCovariance  = None,
1715                  asEyeByScalar = None,
1716                  asEyeByVector = None,
1717                  asCovObject   = None,
1718                  asScript      = None,
1719                  toBeChecked   = False,
1720                 ):
1721         """
1722         Permet de définir une covariance :
1723         - asCovariance : entrée des données, comme une matrice compatible avec
1724           le constructeur de numpy.matrix
1725         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
1726           multiplicatif d'une matrice de corrélation identité, aucune matrice
1727           n'étant donc explicitement à donner
1728         - asEyeByVector : entrée des données comme un seul vecteur de variance,
1729           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
1730           n'étant donc explicitement à donner
1731         - asCovObject : entrée des données comme un objet python, qui a les
1732           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
1733           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
1734           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
1735         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
1736           pleine doit être vérifié
1737         """
1738         self.__name       = str(name)
1739         self.__check      = bool(toBeChecked)
1740         #
1741         self.__C          = None
1742         self.__is_scalar  = False
1743         self.__is_vector  = False
1744         self.__is_matrix  = False
1745         self.__is_object  = False
1746         #
1747         if asScript is not None:
1748             __Matrix, __Scalar, __Vector, __Object = None, None, None, None
1749             if asEyeByScalar:
1750                 __Scalar = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1751             elif asEyeByVector:
1752                 __Vector = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1753             elif asCovObject:
1754                 __Object = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1755             else:
1756                 __Matrix = Interfaces.ImportFromScript(asScript).getvalue( self.__name )
1757         else:
1758             __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject
1759         #
1760         if __Scalar is not None:
1761             if numpy.matrix(__Scalar).size != 1:
1762                 raise ValueError('  The diagonal multiplier given to define a sparse matrix is not a unique scalar value.\n  Its actual measured size is %i. Please check your scalar input.'%numpy.matrix(__Scalar).size)
1763             self.__is_scalar = True
1764             self.__C         = numpy.abs( float(__Scalar) )
1765             self.shape       = (0,0)
1766             self.size        = 0
1767         elif __Vector is not None:
1768             self.__is_vector = True
1769             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(__Vector, float ) ) ) )
1770             self.shape       = (self.__C.size,self.__C.size)
1771             self.size        = self.__C.size**2
1772         elif __Matrix is not None:
1773             self.__is_matrix = True
1774             self.__C         = numpy.matrix( __Matrix, float )
1775             self.shape       = self.__C.shape
1776             self.size        = self.__C.size
1777         elif __Object is not None:
1778             self.__is_object = True
1779             self.__C         = __Object
1780             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__matmul__","__mul__","__rmatmul__","__rmul__"):
1781                 if not hasattr(self.__C,at):
1782                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
1783             if hasattr(self.__C,"shape"):
1784                 self.shape       = self.__C.shape
1785             else:
1786                 self.shape       = (0,0)
1787             if hasattr(self.__C,"size"):
1788                 self.size        = self.__C.size
1789             else:
1790                 self.size        = 0
1791         else:
1792             pass
1793             # raise ValueError("The %s covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix."%self.__name)
1794         #
1795         self.__validate()
1796
1797     def __validate(self):
1798         "Validation"
1799         if self.__C is None:
1800             raise UnboundLocalError("%s covariance matrix value has not been set!"%(self.__name,))
1801         if self.ismatrix() and min(self.shape) != max(self.shape):
1802             raise ValueError("The given matrix for %s is not a square one, its shape is %s. Please check your matrix input."%(self.__name,self.shape))
1803         if self.isobject() and min(self.shape) != max(self.shape):
1804             raise ValueError("The matrix given for \"%s\" is not a square one, its shape is %s. Please check your object input."%(self.__name,self.shape))
1805         if self.isscalar() and self.__C <= 0:
1806             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
1807         if self.isvector() and (self.__C <= 0).any():
1808             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
1809         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
1810             try:
1811                 L = numpy.linalg.cholesky( self.__C )
1812             except:
1813                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
1814         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
1815             try:
1816                 L = self.__C.cholesky()
1817             except:
1818                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
1819
1820     def isscalar(self):
1821         "Vérification du type interne"
1822         return self.__is_scalar
1823
1824     def isvector(self):
1825         "Vérification du type interne"
1826         return self.__is_vector
1827
1828     def ismatrix(self):
1829         "Vérification du type interne"
1830         return self.__is_matrix
1831
1832     def isobject(self):
1833         "Vérification du type interne"
1834         return self.__is_object
1835
1836     def getI(self):
1837         "Inversion"
1838         if   self.ismatrix():
1839             return Covariance(self.__name+"I", asCovariance  = numpy.linalg.inv(self.__C) )
1840         elif self.isvector():
1841             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
1842         elif self.isscalar():
1843             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
1844         elif self.isobject() and hasattr(self.__C,"getI"):
1845             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
1846         else:
1847             return None # Indispensable
1848
1849     def getT(self):
1850         "Transposition"
1851         if   self.ismatrix():
1852             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
1853         elif self.isvector():
1854             return Covariance(self.__name+"T", asEyeByVector = self.__C )
1855         elif self.isscalar():
1856             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
1857         elif self.isobject() and hasattr(self.__C,"getT"):
1858             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
1859         else:
1860             raise AttributeError("the %s covariance matrix has no getT attribute."%(self.__name,))
1861
1862     def cholesky(self):
1863         "Décomposition de Cholesky"
1864         if   self.ismatrix():
1865             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
1866         elif self.isvector():
1867             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
1868         elif self.isscalar():
1869             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
1870         elif self.isobject() and hasattr(self.__C,"cholesky"):
1871             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
1872         else:
1873             raise AttributeError("the %s covariance matrix has no cholesky attribute."%(self.__name,))
1874
1875     def choleskyI(self):
1876         "Inversion de la décomposition de Cholesky"
1877         if   self.ismatrix():
1878             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.inv(numpy.linalg.cholesky(self.__C)) )
1879         elif self.isvector():
1880             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
1881         elif self.isscalar():
1882             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
1883         elif self.isobject() and hasattr(self.__C,"choleskyI"):
1884             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
1885         else:
1886             raise AttributeError("the %s covariance matrix has no choleskyI attribute."%(self.__name,))
1887
1888     def sqrtm(self):
1889         "Racine carrée matricielle"
1890         if   self.ismatrix():
1891             import scipy
1892             return Covariance(self.__name+"C", asCovariance  = numpy.real(scipy.linalg.sqrtm(self.__C)) )
1893         elif self.isvector():
1894             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
1895         elif self.isscalar():
1896             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
1897         elif self.isobject() and hasattr(self.__C,"sqrtm"):
1898             return Covariance(self.__name+"C", asCovObject   = self.__C.sqrtm() )
1899         else:
1900             raise AttributeError("the %s covariance matrix has no sqrtm attribute."%(self.__name,))
1901
1902     def sqrtmI(self):
1903         "Inversion de la racine carrée matricielle"
1904         if   self.ismatrix():
1905             import scipy
1906             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.inv(numpy.real(scipy.linalg.sqrtm(self.__C))) )
1907         elif self.isvector():
1908             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
1909         elif self.isscalar():
1910             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
1911         elif self.isobject() and hasattr(self.__C,"sqrtmI"):
1912             return Covariance(self.__name+"H", asCovObject   = self.__C.sqrtmI() )
1913         else:
1914             raise AttributeError("the %s covariance matrix has no sqrtmI attribute."%(self.__name,))
1915
1916     def diag(self, msize=None):
1917         "Diagonale de la matrice"
1918         if   self.ismatrix():
1919             return numpy.diag(self.__C)
1920         elif self.isvector():
1921             return self.__C
1922         elif self.isscalar():
1923             if msize is None:
1924                 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
1925             else:
1926                 return self.__C * numpy.ones(int(msize))
1927         elif self.isobject() and hasattr(self.__C,"diag"):
1928             return self.__C.diag()
1929         else:
1930             raise AttributeError("the %s covariance matrix has no diag attribute."%(self.__name,))
1931
1932     def asfullmatrix(self, msize=None):
1933         "Matrice pleine"
1934         if   self.ismatrix():
1935             return numpy.asarray(self.__C)
1936         elif self.isvector():
1937             return numpy.asarray( numpy.diag(self.__C), float )
1938         elif self.isscalar():
1939             if msize is None:
1940                 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
1941             else:
1942                 return numpy.asarray( self.__C * numpy.eye(int(msize)), float )
1943         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
1944             return self.__C.asfullmatrix()
1945         else:
1946             raise AttributeError("the %s covariance matrix has no asfullmatrix attribute."%(self.__name,))
1947
1948     def trace(self, msize=None):
1949         "Trace de la matrice"
1950         if   self.ismatrix():
1951             return numpy.trace(self.__C)
1952         elif self.isvector():
1953             return float(numpy.sum(self.__C))
1954         elif self.isscalar():
1955             if msize is None:
1956                 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
1957             else:
1958                 return self.__C * int(msize)
1959         elif self.isobject():
1960             return self.__C.trace()
1961         else:
1962             raise AttributeError("the %s covariance matrix has no trace attribute."%(self.__name,))
1963
1964     def getO(self):
1965         return self
1966
1967     def __repr__(self):
1968         "x.__repr__() <==> repr(x)"
1969         return repr(self.__C)
1970
1971     def __str__(self):
1972         "x.__str__() <==> str(x)"
1973         return str(self.__C)
1974
1975     def __add__(self, other):
1976         "x.__add__(y) <==> x+y"
1977         if   self.ismatrix() or self.isobject():
1978             return self.__C + numpy.asmatrix(other)
1979         elif self.isvector() or self.isscalar():
1980             _A = numpy.asarray(other)
1981             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
1982             return numpy.asmatrix(_A)
1983
1984     def __radd__(self, other):
1985         "x.__radd__(y) <==> y+x"
1986         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
1987
1988     def __sub__(self, other):
1989         "x.__sub__(y) <==> x-y"
1990         if   self.ismatrix() or self.isobject():
1991             return self.__C - numpy.asmatrix(other)
1992         elif self.isvector() or self.isscalar():
1993             _A = numpy.asarray(other)
1994             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
1995             return numpy.asmatrix(_A)
1996
1997     def __rsub__(self, other):
1998         "x.__rsub__(y) <==> y-x"
1999         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
2000
2001     def __neg__(self):
2002         "x.__neg__() <==> -x"
2003         return - self.__C
2004
2005     def __matmul__(self, other):
2006         "x.__mul__(y) <==> x@y"
2007         if   self.ismatrix() and isinstance(other, (int, float)):
2008             return numpy.asarray(self.__C) * other
2009         elif self.ismatrix() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
2010             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2011                 return numpy.ravel(self.__C @ numpy.ravel(other))
2012             elif numpy.asarray(other).shape[0] == self.shape[1]: # Matrice
2013                 return numpy.asarray(self.__C) @ numpy.asarray(other)
2014             else:
2015                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asarray(other).shape,self.__name))
2016         elif self.isvector() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
2017             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2018                 return numpy.ravel(self.__C) * numpy.ravel(other)
2019             elif numpy.asarray(other).shape[0] == self.shape[1]: # Matrice
2020                 return numpy.ravel(self.__C).reshape((-1,1)) * numpy.asarray(other)
2021             else:
2022                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
2023         elif self.isscalar() and isinstance(other,numpy.matrix):
2024             return numpy.asarray(self.__C * other)
2025         elif self.isscalar() and isinstance(other, (list, numpy.ndarray, tuple)):
2026             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
2027                 return self.__C * numpy.ravel(other)
2028             else:
2029                 return self.__C * numpy.asarray(other)
2030         elif self.isobject():
2031             return self.__C.__matmul__(other)
2032         else:
2033             raise NotImplementedError("%s covariance matrix __matmul__ method not available for %s type!"%(self.__name,type(other)))
2034
2035     def __mul__(self, other):
2036         "x.__mul__(y) <==> x*y"
2037         if   self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
2038             return self.__C * other
2039         elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
2040             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2041                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
2042             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
2043                 return self.__C * numpy.asmatrix(other)
2044             else:
2045                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
2046         elif self.isvector() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
2047             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2048                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
2049             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
2050                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
2051             else:
2052                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
2053         elif self.isscalar() and isinstance(other,numpy.matrix):
2054             return self.__C * other
2055         elif self.isscalar() and isinstance(other, (list, numpy.ndarray, tuple)):
2056             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
2057                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
2058             else:
2059                 return self.__C * numpy.asmatrix(other)
2060         elif self.isobject():
2061             return self.__C.__mul__(other)
2062         else:
2063             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
2064
2065     def __rmatmul__(self, other):
2066         "x.__rmul__(y) <==> y@x"
2067         if self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
2068             return other * self.__C
2069         elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
2070             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2071                 return numpy.asmatrix(numpy.ravel(other)) * self.__C
2072             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
2073                 return numpy.asmatrix(other) * self.__C
2074             else:
2075                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.asmatrix(other).shape,self.shape,self.__name))
2076         elif self.isvector() and isinstance(other,numpy.matrix):
2077             if numpy.ravel(other).size == self.shape[0]: # Vecteur
2078                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
2079             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
2080                 return numpy.asmatrix(numpy.array(other) * self.__C)
2081             else:
2082                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.ravel(other).shape,self.shape,self.__name))
2083         elif self.isscalar() and isinstance(other,numpy.matrix):
2084             return other * self.__C
2085         elif self.isobject():
2086             return self.__C.__rmatmul__(other)
2087         else:
2088             raise NotImplementedError("%s covariance matrix __rmatmul__ method not available for %s type!"%(self.__name,type(other)))
2089
2090     def __rmul__(self, other):
2091         "x.__rmul__(y) <==> y*x"
2092         if self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
2093             return other * self.__C
2094         elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
2095             if numpy.ravel(other).size == self.shape[1]: # Vecteur
2096                 return numpy.asmatrix(numpy.ravel(other)) * self.__C
2097             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
2098                 return numpy.asmatrix(other) * self.__C
2099             else:
2100                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.asmatrix(other).shape,self.shape,self.__name))
2101         elif self.isvector() and isinstance(other,numpy.matrix):
2102             if numpy.ravel(other).size == self.shape[0]: # Vecteur
2103                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
2104             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
2105                 return numpy.asmatrix(numpy.array(other) * self.__C)
2106             else:
2107                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.ravel(other).shape,self.shape,self.__name))
2108         elif self.isscalar() and isinstance(other,numpy.matrix):
2109             return other * self.__C
2110         elif self.isobject():
2111             return self.__C.__rmul__(other)
2112         else:
2113             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
2114
2115     def __len__(self):
2116         "x.__len__() <==> len(x)"
2117         return self.shape[0]
2118
2119 # ==============================================================================
2120 class Observer2Func(object):
2121     """
2122     Creation d'une fonction d'observateur a partir de son texte
2123     """
2124     def __init__(self, corps=""):
2125         self.__corps = corps
2126     def func(self,var,info):
2127         "Fonction d'observation"
2128         exec(self.__corps)
2129     def getfunc(self):
2130         "Restitution du pointeur de fonction dans l'objet"
2131         return self.func
2132
2133 # ==============================================================================
2134 class CaseLogger(object):
2135     """
2136     Conservation des commandes de creation d'un cas
2137     """
2138     def __init__(self, __name="", __objname="case", __addViewers=None, __addLoaders=None):
2139         self.__name     = str(__name)
2140         self.__objname  = str(__objname)
2141         self.__logSerie = []
2142         self.__switchoff = False
2143         self.__viewers = {
2144             "TUI" :Interfaces._TUIViewer,
2145             "SCD" :Interfaces._SCDViewer,
2146             "YACS":Interfaces._YACSViewer,
2147             }
2148         self.__loaders = {
2149             "TUI" :Interfaces._TUIViewer,
2150             "COM" :Interfaces._COMViewer,
2151             }
2152         if __addViewers is not None:
2153             self.__viewers.update(dict(__addViewers))
2154         if __addLoaders is not None:
2155             self.__loaders.update(dict(__addLoaders))
2156
2157     def register(self, __command=None, __keys=None, __local=None, __pre=None, __switchoff=False):
2158         "Enregistrement d'une commande individuelle"
2159         if __command is not None and __keys is not None and __local is not None and not self.__switchoff:
2160             if "self" in __keys: __keys.remove("self")
2161             self.__logSerie.append( (str(__command), __keys, __local, __pre, __switchoff) )
2162             if __switchoff:
2163                 self.__switchoff = True
2164         if not __switchoff:
2165             self.__switchoff = False
2166
2167     def dump(self, __filename=None, __format="TUI", __upa=""):
2168         "Restitution normalisée des commandes"
2169         if __format in self.__viewers:
2170             __formater = self.__viewers[__format](self.__name, self.__objname, self.__logSerie)
2171         else:
2172             raise ValueError("Dumping as \"%s\" is not available"%__format)
2173         return __formater.dump(__filename, __upa)
2174
2175     def load(self, __filename=None, __content=None, __object=None, __format="TUI"):
2176         "Chargement normalisé des commandes"
2177         if __format in self.__loaders:
2178             __formater = self.__loaders[__format]()
2179         else:
2180             raise ValueError("Loading as \"%s\" is not available"%__format)
2181         return __formater.load(__filename, __content, __object)
2182
2183 # ==============================================================================
2184 def MultiFonction(
2185         __xserie,
2186         _extraArguments = None,
2187         _sFunction      = lambda x: x,
2188         _mpEnabled      = False,
2189         _mpWorkers      = None,
2190         ):
2191     """
2192     Pour une liste ordonnée de vecteurs en entrée, renvoie en sortie la liste
2193     correspondante de valeurs de la fonction en argument
2194     """
2195     # Vérifications et définitions initiales
2196     # logging.debug("MULTF Internal multifonction calculations begin with function %s"%(_sFunction.__name__,))
2197     if not PlatformInfo.isIterable( __xserie ):
2198         raise TypeError("MultiFonction not iterable unkown input type: %s"%(type(__xserie),))
2199     if _mpEnabled:
2200         if (_mpWorkers is None) or (_mpWorkers is not None and _mpWorkers < 1):
2201             __mpWorkers = None
2202         else:
2203             __mpWorkers = int(_mpWorkers)
2204         try:
2205             import multiprocessing
2206             __mpEnabled = True
2207         except ImportError:
2208             __mpEnabled = False
2209     else:
2210         __mpEnabled = False
2211         __mpWorkers = None
2212     #
2213     # Calculs effectifs
2214     if __mpEnabled:
2215         _jobs = __xserie
2216         # logging.debug("MULTF Internal multiprocessing calculations begin : evaluation of %i point(s)"%(len(_jobs),))
2217         import multiprocessing
2218         with multiprocessing.Pool(__mpWorkers) as pool:
2219             __multiHX = pool.map( _sFunction, _jobs )
2220             pool.close()
2221             pool.join()
2222         # logging.debug("MULTF Internal multiprocessing calculation end")
2223     else:
2224         # logging.debug("MULTF Internal monoprocessing calculation begin")
2225         __multiHX = []
2226         if _extraArguments is None:
2227             for __xvalue in __xserie:
2228                 __multiHX.append( _sFunction( __xvalue ) )
2229         elif _extraArguments is not None and isinstance(_extraArguments, (list, tuple, map)):
2230             for __xvalue in __xserie:
2231                 __multiHX.append( _sFunction( __xvalue, *_extraArguments ) )
2232         elif _extraArguments is not None and isinstance(_extraArguments, dict):
2233             for __xvalue in __xserie:
2234                 __multiHX.append( _sFunction( __xvalue, **_extraArguments ) )
2235         else:
2236             raise TypeError("MultiFonction extra arguments unkown input type: %s"%(type(_extraArguments),))
2237         # logging.debug("MULTF Internal monoprocessing calculation end")
2238     #
2239     # logging.debug("MULTF Internal multifonction calculations end")
2240     return __multiHX
2241
2242 # ==============================================================================
2243 def CostFunction3D(_x,
2244                    _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
2245                    _HmX = None,  # Simulation déjà faite de Hm(x)
2246                    _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
2247                    _BI  = None,
2248                    _RI  = None,
2249                    _Xb  = None,
2250                    _Y   = None,
2251                    _SIV = False, # A résorber pour la 8.0
2252                    _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
2253                    _nPS = 0,     # nbPreviousSteps
2254                    _QM  = "DA",  # QualityMeasure
2255                    _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
2256                    _fRt = False, # Restitue ou pas la sortie étendue
2257                    _sSc = True,  # Stocke ou pas les SSC
2258                   ):
2259     """
2260     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
2261     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
2262     DFO, QuantileRegression
2263     """
2264     if not _sSc:
2265         _SIV = False
2266         _SSC = {}
2267     else:
2268         for k in ["CostFunctionJ",
2269                   "CostFunctionJb",
2270                   "CostFunctionJo",
2271                   "CurrentOptimum",
2272                   "CurrentState",
2273                   "IndexOfOptimum",
2274                   "SimulatedObservationAtCurrentOptimum",
2275                   "SimulatedObservationAtCurrentState",
2276                  ]:
2277             if k not in _SSV:
2278                 _SSV[k] = []
2279             if hasattr(_SSV[k],"store"):
2280                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
2281     #
2282     _X  = numpy.asmatrix(numpy.ravel( _x )).T
2283     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
2284         _SSV["CurrentState"].append( _X )
2285     #
2286     if _HmX is not None:
2287         _HX = _HmX
2288     else:
2289         if _Hm is None:
2290             raise ValueError("COSTFUNCTION3D Operator has to be defined.")
2291         if _arg is None:
2292             _HX = _Hm( _X )
2293         else:
2294             _HX = _Hm( _X, *_arg )
2295     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
2296     #
2297     if "SimulatedObservationAtCurrentState" in _SSC or \
2298        "SimulatedObservationAtCurrentOptimum" in _SSC:
2299         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
2300     #
2301     if numpy.any(numpy.isnan(_HX)):
2302         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
2303     else:
2304         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
2305         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
2306             if _BI is None or _RI is None:
2307                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
2308             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
2309             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
2310             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
2311         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
2312             if _RI is None:
2313                 raise ValueError("Observation error covariance matrix has to be properly defined!")
2314             Jb  = 0.
2315             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
2316         elif _QM in ["LeastSquares", "LS", "L2"]:
2317             Jb  = 0.
2318             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
2319         elif _QM in ["AbsoluteValue", "L1"]:
2320             Jb  = 0.
2321             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
2322         elif _QM in ["MaximumError", "ME"]:
2323             Jb  = 0.
2324             Jo  = numpy.max( numpy.abs(_Y - _HX) )
2325         elif _QM in ["QR", "Null"]:
2326             Jb  = 0.
2327             Jo  = 0.
2328         else:
2329             raise ValueError("Unknown asked quality measure!")
2330         #
2331         J   = float( Jb ) + float( Jo )
2332     #
2333     if _sSc:
2334         _SSV["CostFunctionJb"].append( Jb )
2335         _SSV["CostFunctionJo"].append( Jo )
2336         _SSV["CostFunctionJ" ].append( J )
2337     #
2338     if "IndexOfOptimum" in _SSC or \
2339        "CurrentOptimum" in _SSC or \
2340        "SimulatedObservationAtCurrentOptimum" in _SSC:
2341         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
2342     if "IndexOfOptimum" in _SSC:
2343         _SSV["IndexOfOptimum"].append( IndexMin )
2344     if "CurrentOptimum" in _SSC:
2345         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
2346     if "SimulatedObservationAtCurrentOptimum" in _SSC:
2347         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
2348     #
2349     if _fRt:
2350         return _SSV
2351     else:
2352         if _QM in ["QR"]: # Pour le QuantileRegression
2353             return _HX
2354         else:
2355             return J
2356
2357 # ==============================================================================
2358 if __name__ == "__main__":
2359     print('\n AUTODIAGNOSTIC\n')