src/daComposant/daCore/NumericObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2023 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 __doc__ = """
  24     Définit les objets numériques génériques.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, copy, types, sys, logging, numpy, itertools
  29 from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm
  30 from daCore.PlatformInfo import PlatformInfo
  31 mpr = PlatformInfo().MachinePrecision()
  32 mfp = PlatformInfo().MaximumPrecision()
  33 # logging.getLogger().setLevel(logging.DEBUG)
  34
  35 # ==============================================================================
  36 def ExecuteFunction( triplet ):
  37     assert len(triplet) == 3, "Incorrect number of arguments"
  38     X, xArgs, funcrepr = triplet
  39     __X = numpy.ravel( X ).reshape((-1,1))
  40     __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
  41     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  42     __fonction = getattr(__module,funcrepr["__userFunction__name"])
  43     sys.path = __sys_path_tmp ; del __sys_path_tmp
  44     if isinstance(xArgs, dict):
  45         __HX  = __fonction( __X, **xArgs )
  46     else:
  47         __HX  = __fonction( __X )
  48     return numpy.ravel( __HX )
  49
  50 # ==============================================================================
  51 class FDApproximation(object):
  52     """
  53     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  54     création d'un objet, en fournissant une fonction "Function", on obtient un
  55     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  56     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  57     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  58     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  59     centrées si le booléen "centeredDF" est vrai.
  60     """
  61     def __init__(self,
  62             name                  = "FDApproximation",
  63             Function              = None,
  64             centeredDF            = False,
  65             increment             = 0.01,
  66             dX                    = None,
  67             extraArguments        = None,
  68             reducingMemoryUse     = False,
  69             avoidingRedundancy    = True,
  70             toleranceInRedundancy = 1.e-18,
  71             lengthOfRedundancy    = -1,
  72             mpEnabled             = False,
  73             mpWorkers             = None,
  74             mfEnabled             = False,
  75             ):
  76         self.__name = str(name)
  77         self.__extraArgs = extraArguments
  78         #
  79         if mpEnabled:
  80             try:
  81                 import multiprocessing
  82                 self.__mpEnabled = True
  83             except ImportError:
  84                 self.__mpEnabled = False
  85         else:
  86             self.__mpEnabled = False
  87         self.__mpWorkers = mpWorkers
  88         if self.__mpWorkers is not None and self.__mpWorkers < 1:
  89             self.__mpWorkers = None
  90         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
  91         #
  92         self.__mfEnabled = bool(mfEnabled)
  93         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
  94         #
  95         self.__rmEnabled = bool(reducingMemoryUse)
  96         logging.debug("FDA Calculs avec réduction mémoire : %s"%(self.__rmEnabled,))
  97         #
  98         if avoidingRedundancy:
  99             self.__avoidRC = True
 100             self.__tolerBP = float(toleranceInRedundancy)
 101             self.__lengthRJ = int(lengthOfRedundancy)
 102             self.__listJPCP = [] # Jacobian Previous Calculated Points
 103             self.__listJPCI = [] # Jacobian Previous Calculated Increment
 104             self.__listJPCR = [] # Jacobian Previous Calculated Results
 105             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
 106             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
 107         else:
 108             self.__avoidRC = False
 109         logging.debug("FDA Calculs avec réduction des doublons : %s"%self.__avoidRC)
 110         if self.__avoidRC:
 111             logging.debug("FDA Tolérance de détermination des doublons : %.2e"%self.__tolerBP)
 112         #
 113         if self.__mpEnabled:
 114             if isinstance(Function,types.FunctionType):
 115                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 116                 self.__userFunction__name = Function.__name__
 117                 try:
 118                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 119                 except Exception:
 120                     mod = os.path.abspath(Function.__globals__['__file__'])
 121                 if not os.path.isfile(mod):
 122                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 123                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 124                 self.__userFunction__path = os.path.dirname(mod)
 125                 del mod
 126                 self.__userOperator = Operator(
 127                     name                 = self.__name,
 128                     fromMethod           = Function,
 129                     avoidingRedundancy   = self.__avoidRC,
 130                     inputAsMultiFunction = self.__mfEnabled,
 131                     extraArguments       = self.__extraArgs )
 132                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 133             elif isinstance(Function,types.MethodType):
 134                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 135                 self.__userFunction__name = Function.__name__
 136                 try:
 137                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 138                 except Exception:
 139                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 140                 if not os.path.isfile(mod):
 141                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 142                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 143                 self.__userFunction__path = os.path.dirname(mod)
 144                 del mod
 145                 self.__userOperator = Operator(
 146                     name                 = self.__name,
 147                     fromMethod           = Function,
 148                     avoidingRedundancy   = self.__avoidRC,
 149                     inputAsMultiFunction = self.__mfEnabled,
 150                     extraArguments       = self.__extraArgs )
 151                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 152             else:
 153                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 154         else:
 155             self.__userOperator = Operator(
 156                 name                 = self.__name,
 157                 fromMethod           = Function,
 158                 avoidingRedundancy   = self.__avoidRC,
 159                 inputAsMultiFunction = self.__mfEnabled,
 160                 extraArguments       = self.__extraArgs )
 161             self.__userFunction = self.__userOperator.appliedTo
 162         #
 163         self.__centeredDF = bool(centeredDF)
 164         if abs(float(increment)) > 1.e-15:
 165             self.__increment  = float(increment)
 166         else:
 167             self.__increment  = 0.01
 168         if dX is None:
 169             self.__dX     = None
 170         else:
 171             self.__dX     = numpy.ravel( dX )
 172
 173     # ---------------------------------------------------------
 174     def __doublon__(self, __e, __l, __n, __v=None):
 175         __ac, __iac = False, -1
 176         for i in range(len(__l)-1,-1,-1):
 177             if numpy.linalg.norm(__e - __l[i]) < self.__tolerBP * __n[i]:
 178                 __ac, __iac = True, i
 179                 if __v is not None: logging.debug("FDA Cas%s déjà calculé, récupération du doublon %i"%(__v,__iac))
 180                 break
 181         return __ac, __iac
 182
 183     # ---------------------------------------------------------
 184     def __listdotwith__(self, __LMatrix, __dotWith = None, __dotTWith = None):
 185         "Produit incrémental d'une matrice liste de colonnes avec un vecteur"
 186         if not isinstance(__LMatrix, (list,tuple)):
 187             raise TypeError("Columnwise list matrix has not the proper type: %s"%type(__LMatrix))
 188         if __dotWith is not None:
 189             __Idwx = numpy.ravel( __dotWith )
 190             assert len(__LMatrix) == __Idwx.size, "Incorrect size of elements"
 191             __Produit = numpy.zeros(__LMatrix[0].size)
 192             for i, col in enumerate(__LMatrix):
 193                 __Produit += float(__Idwx[i]) * col
 194             return __Produit
 195         elif __dotTWith is not None:
 196             _Idwy = numpy.ravel( __dotTWith ).T
 197             assert __LMatrix[0].size == _Idwy.size, "Incorrect size of elements"
 198             __Produit = numpy.zeros(len(__LMatrix))
 199             for i, col in enumerate(__LMatrix):
 200                 __Produit[i] = float( _Idwy @ col)
 201             return __Produit
 202         else:
 203             __Produit = None
 204         return __Produit
 205
 206     # ---------------------------------------------------------
 207     def DirectOperator(self, X, **extraArgs ):
 208         """
 209         Calcul du direct à l'aide de la fonction fournie.
 210
 211         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 212         ne doivent pas être données ici à la fonction utilisateur.
 213         """
 214         logging.debug("FDA Calcul DirectOperator (explicite)")
 215         if self.__mfEnabled:
 216             _HX = self.__userFunction( X, argsAsSerie = True )
 217         else:
 218             _HX = numpy.ravel(self.__userFunction( numpy.ravel(X) ))
 219         #
 220         return _HX
 221
 222     # ---------------------------------------------------------
 223     def TangentMatrix(self, X, dotWith = None, dotTWith = None ):
 224         """
 225         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 226         c'est-à-dire le gradient de H en X. On utilise des différences finies
 227         directionnelles autour du point X. X est un numpy.ndarray.
 228
 229         Différences finies centrées (approximation d'ordre 2):
 230         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 231            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 232            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 233            H( X_moins_dXi )
 234         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 235            le pas 2*dXi
 236         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 237
 238         Différences finies non centrées (approximation d'ordre 1):
 239         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 240            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 241            HX_plus_dXi = H( X_plus_dXi )
 242         2/ On calcule la valeur centrale HX = H(X)
 243         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 244            le pas dXi
 245         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 246
 247         """
 248         logging.debug("FDA Début du calcul de la Jacobienne")
 249         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 250         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 251         #
 252         if X is None or len(X)==0:
 253             raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
 254         #
 255         _X = numpy.ravel( X )
 256         #
 257         if self.__dX is None:
 258             _dX  = self.__increment * _X
 259         else:
 260             _dX = numpy.ravel( self.__dX )
 261         assert len(_X) == len(_dX), "Inconsistent dX increment length with respect to the X one"
 262         assert _X.size == _dX.size, "Inconsistent dX increment size with respect to the X one"
 263         #
 264         if (_dX == 0.).any():
 265             moyenne = _dX.mean()
 266             if moyenne == 0.:
 267                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 268             else:
 269                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 270         #
 271         __alreadyCalculated  = False
 272         if self.__avoidRC:
 273             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 274             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 275             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 276                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 277                 logging.debug("FDA Cas J déjà calculé, récupération du doublon %i"%__i)
 278         #
 279         if __alreadyCalculated:
 280             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 281             _Jacobienne = self.__listJPCR[__i]
 282             logging.debug("FDA Fin du calcul de la Jacobienne")
 283             if dotWith is not None:
 284                 return numpy.dot(_Jacobienne,   numpy.ravel( dotWith ))
 285             elif dotTWith is not None:
 286                 return numpy.dot(_Jacobienne.T, numpy.ravel( dotTWith ))
 287         else:
 288             logging.debug("FDA   Calcul Jacobienne (explicite)")
 289             if self.__centeredDF:
 290                 #
 291                 if self.__mpEnabled and not self.__mfEnabled:
 292                     funcrepr = {
 293                         "__userFunction__path" : self.__userFunction__path,
 294                         "__userFunction__modl" : self.__userFunction__modl,
 295                         "__userFunction__name" : self.__userFunction__name,
 296                     }
 297                     _jobs = []
 298                     for i in range( len(_dX) ):
 299                         _dXi            = _dX[i]
 300                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 301                         _X_plus_dXi[i]  = _X[i] + _dXi
 302                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 303                         _X_moins_dXi[i] = _X[i] - _dXi
 304                         #
 305                         _jobs.append( (_X_plus_dXi,  self.__extraArgs, funcrepr) )
 306                         _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
 307                     #
 308                     import multiprocessing
 309                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 310                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 311                     self.__pool.close()
 312                     self.__pool.join()
 313                     #
 314                     _Jacobienne  = []
 315                     for i in range( len(_dX) ):
 316                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 317                     #
 318                 elif self.__mfEnabled:
 319                     _xserie = []
 320                     for i in range( len(_dX) ):
 321                         _dXi            = _dX[i]
 322                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 323                         _X_plus_dXi[i]  = _X[i] + _dXi
 324                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 325                         _X_moins_dXi[i] = _X[i] - _dXi
 326                         #
 327                         _xserie.append( _X_plus_dXi )
 328                         _xserie.append( _X_moins_dXi )
 329                     #
 330                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 331                     #
 332                     _Jacobienne  = []
 333                     for i in range( len(_dX) ):
 334                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 335                     #
 336                 else:
 337                     _Jacobienne  = []
 338                     for i in range( _dX.size ):
 339                         _dXi            = _dX[i]
 340                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 341                         _X_plus_dXi[i]  = _X[i] + _dXi
 342                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 343                         _X_moins_dXi[i] = _X[i] - _dXi
 344                         #
 345                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 346                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 347                         #
 348                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 349                 #
 350             else:
 351                 #
 352                 if self.__mpEnabled and not self.__mfEnabled:
 353                     funcrepr = {
 354                         "__userFunction__path" : self.__userFunction__path,
 355                         "__userFunction__modl" : self.__userFunction__modl,
 356                         "__userFunction__name" : self.__userFunction__name,
 357                     }
 358                     _jobs = []
 359                     _jobs.append( (_X, self.__extraArgs, funcrepr) )
 360                     for i in range( len(_dX) ):
 361                         _X_plus_dXi    = numpy.array( _X, dtype=float )
 362                         _X_plus_dXi[i] = _X[i] + _dX[i]
 363                         #
 364                         _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
 365                     #
 366                     import multiprocessing
 367                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 368                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 369                     self.__pool.close()
 370                     self.__pool.join()
 371                     #
 372                     _HX = _HX_plus_dX.pop(0)
 373                     #
 374                     _Jacobienne = []
 375                     for i in range( len(_dX) ):
 376                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 377                     #
 378                 elif self.__mfEnabled:
 379                     _xserie = []
 380                     _xserie.append( _X )
 381                     for i in range( len(_dX) ):
 382                         _X_plus_dXi    = numpy.array( _X, dtype=float )
 383                         _X_plus_dXi[i] = _X[i] + _dX[i]
 384                         #
 385                         _xserie.append( _X_plus_dXi )
 386                     #
 387                     _HX_plus_dX = self.DirectOperator( _xserie )
 388                     #
 389                     _HX = _HX_plus_dX.pop(0)
 390                     #
 391                     _Jacobienne = []
 392                     for i in range( len(_dX) ):
 393                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 394                     #
 395                 else:
 396                     _Jacobienne  = []
 397                     _HX = self.DirectOperator( _X )
 398                     for i in range( _dX.size ):
 399                         _dXi            = _dX[i]
 400                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 401                         _X_plus_dXi[i]  = _X[i] + _dXi
 402                         #
 403                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 404                         #
 405                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 406             #
 407             if (dotWith is not None) or (dotTWith is not None):
 408                 __Produit = self.__listdotwith__(_Jacobienne, dotWith, dotTWith)
 409             else:
 410                 __Produit = None
 411             if __Produit is None or self.__avoidRC:
 412                 _Jacobienne = numpy.transpose( numpy.vstack( _Jacobienne ) )
 413                 if self.__avoidRC:
 414                     if self.__lengthRJ < 0: self.__lengthRJ = 2 * _X.size
 415                     while len(self.__listJPCP) > self.__lengthRJ:
 416                         self.__listJPCP.pop(0)
 417                         self.__listJPCI.pop(0)
 418                         self.__listJPCR.pop(0)
 419                         self.__listJPPN.pop(0)
 420                         self.__listJPIN.pop(0)
 421                     self.__listJPCP.append( copy.copy(_X) )
 422                     self.__listJPCI.append( copy.copy(_dX) )
 423                     self.__listJPCR.append( copy.copy(_Jacobienne) )
 424                     self.__listJPPN.append( numpy.linalg.norm(_X) )
 425                     self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 426             logging.debug("FDA Fin du calcul de la Jacobienne")
 427             if __Produit is not None:
 428                 return __Produit
 429         #
 430         return _Jacobienne
 431
 432     # ---------------------------------------------------------
 433     def TangentOperator(self, paire, **extraArgs ):
 434         """
 435         Calcul du tangent à l'aide de la Jacobienne.
 436
 437         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 438         ne doivent pas être données ici à la fonction utilisateur.
 439         """
 440         if self.__mfEnabled:
 441             assert len(paire) == 1, "Incorrect length of arguments"
 442             _paire = paire[0]
 443             assert len(_paire) == 2, "Incorrect number of arguments"
 444         else:
 445             assert len(paire) == 2, "Incorrect number of arguments"
 446             _paire = paire
 447         X, dX = _paire
 448         if dX is None or len(dX) == 0:
 449             #
 450             # Calcul de la forme matricielle si le second argument est None
 451             # -------------------------------------------------------------
 452             _Jacobienne = self.TangentMatrix( X )
 453             if self.__mfEnabled: return [_Jacobienne,]
 454             else:                return _Jacobienne
 455         else:
 456             #
 457             # Calcul de la valeur linéarisée de H en X appliqué à dX
 458             # ------------------------------------------------------
 459             _HtX = self.TangentMatrix( X, dotWith = dX )
 460             if self.__mfEnabled: return [_HtX,]
 461             else:                return _HtX
 462
 463     # ---------------------------------------------------------
 464     def AdjointOperator(self, paire, **extraArgs ):
 465         """
 466         Calcul de l'adjoint à l'aide de la Jacobienne.
 467
 468         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 469         ne doivent pas être données ici à la fonction utilisateur.
 470         """
 471         if self.__mfEnabled:
 472             assert len(paire) == 1, "Incorrect length of arguments"
 473             _paire = paire[0]
 474             assert len(_paire) == 2, "Incorrect number of arguments"
 475         else:
 476             assert len(paire) == 2, "Incorrect number of arguments"
 477             _paire = paire
 478         X, Y = _paire
 479         if Y is None or len(Y) == 0:
 480             #
 481             # Calcul de la forme matricielle si le second argument est None
 482             # -------------------------------------------------------------
 483             _JacobienneT = self.TangentMatrix( X ).T
 484             if self.__mfEnabled: return [_JacobienneT,]
 485             else:                return _JacobienneT
 486         else:
 487             #
 488             # Calcul de la valeur de l'adjoint en X appliqué à Y
 489             # --------------------------------------------------
 490             _HaY = self.TangentMatrix( X, dotTWith = Y )
 491             if self.__mfEnabled: return [_HaY,]
 492             else:                return _HaY
 493
 494 # ==============================================================================
 495 def SetInitialDirection( __Direction = [], __Amplitude = 1., __Position = None ):
 496     "Établit ou élabore une direction avec une amplitude"
 497     #
 498     if len(__Direction) == 0 and __Position is None:
 499         raise ValueError("If initial direction is void, current position has to be given")
 500     if abs(float(__Amplitude)) < mpr:
 501         raise ValueError("Amplitude of perturbation can not be zero")
 502     #
 503     if len(__Direction) > 0:
 504         __dX0 = numpy.asarray(__Direction)
 505     else:
 506         __dX0 = []
 507         __X0 = numpy.ravel(numpy.asarray(__Position))
 508         __mX0 = numpy.mean( __X0, dtype=mfp )
 509         if abs(__mX0) < 2*mpr: __mX0 = 1. # Évite le problème de position nulle
 510         for v in __X0:
 511             if abs(v) > 1.e-8:
 512                 __dX0.append( numpy.random.normal(0.,abs(v)) )
 513             else:
 514                 __dX0.append( numpy.random.normal(0.,__mX0) )
 515     #
 516     __dX0 = numpy.asarray(__dX0,float) # Évite le problème d'array de taille 1
 517     __dX0 = numpy.ravel( __dX0 )       # Redresse les vecteurs
 518     __dX0 = float(__Amplitude) * __dX0
 519     #
 520     return __dX0
 521
 522 # ==============================================================================
 523 def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
 524     "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
 525     #
 526     __bgCenter = numpy.ravel(__bgCenter)[:,None]
 527     if __nbMembers < 1:
 528         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
 529     #
 530     if __bgCovariance is None:
 531         _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 532     else:
 533         _Z = numpy.random.multivariate_normal(numpy.zeros(__bgCenter.size), __bgCovariance, size=__nbMembers).T
 534         _Perturbations = numpy.tile( __bgCenter, __nbMembers) + _Z
 535     #
 536     return _Perturbations
 537
 538 # ==============================================================================
 539 def EnsembleOfBackgroundPerturbations( __bgCenter, __bgCovariance, __nbMembers, __withSVD = True):
 540     "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
 541     def __CenteredRandomAnomalies(Zr, N):
 542         """
 543         Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
 544         notes manuscrites de MB et conforme au code de PS avec eps = -1
 545         """
 546         eps = -1
 547         Q = numpy.identity(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
 548         Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
 549         R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
 550         Q = numpy.dot(Q,R)
 551         Zr = numpy.dot(Q,Zr)
 552         return Zr.T
 553     #
 554     __bgCenter = numpy.ravel(__bgCenter).reshape((-1,1))
 555     if __nbMembers < 1:
 556         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
 557     if __bgCovariance is None:
 558         _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 559     else:
 560         if __withSVD:
 561             _U, _s, _V = numpy.linalg.svd(__bgCovariance, full_matrices=False)
 562             _nbctl = __bgCenter.size
 563             if __nbMembers > _nbctl:
 564                 _Z = numpy.concatenate((numpy.dot(
 565                     numpy.diag(numpy.sqrt(_s[:_nbctl])), _V[:_nbctl]),
 566                     numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1-_nbctl)), axis = 0)
 567             else:
 568                 _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:__nbMembers-1])), _V[:__nbMembers-1])
 569             _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
 570             _Perturbations = __bgCenter + _Zca
 571         else:
 572             if max(abs(__bgCovariance.flatten())) > 0:
 573                 _nbctl = __bgCenter.size
 574                 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1)
 575                 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
 576                 _Perturbations = __bgCenter + _Zca
 577             else:
 578                 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 579     #
 580     return _Perturbations
 581
 582 # ==============================================================================
 583 def EnsembleMean( __Ensemble ):
 584     "Renvoie la moyenne empirique d'un ensemble"
 585     return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
 586
 587 # ==============================================================================
 588 def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.):
 589     "Renvoie les anomalies centrées à partir d'un ensemble"
 590     if __OptMean is None:
 591         __Em = EnsembleMean( __Ensemble )
 592     else:
 593         __Em = numpy.ravel( __OptMean ).reshape((-1,1))
 594     #
 595     return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
 596
 597 # ==============================================================================
 598 def EnsembleErrorCovariance( __Ensemble, __Quick = False ):
 599     "Renvoie l'estimation empirique de la covariance d'ensemble"
 600     if __Quick:
 601         # Covariance rapide mais rarement définie positive
 602         __Covariance = numpy.cov( __Ensemble )
 603     else:
 604         # Résultat souvent identique à numpy.cov, mais plus robuste
 605         __n, __m = numpy.asarray( __Ensemble ).shape
 606         __Anomalies = EnsembleOfAnomalies( __Ensemble )
 607         # Estimation empirique
 608         __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m-1)
 609         # Assure la symétrie
 610         __Covariance = ( __Covariance + __Covariance.T ) * 0.5
 611         # Assure la positivité
 612         __epsilon    = mpr*numpy.trace( __Covariance )
 613         __Covariance = __Covariance + __epsilon * numpy.identity(__n)
 614     #
 615     return __Covariance
 616
 617 # ==============================================================================
 618 def EnsemblePerturbationWithGivenCovariance(
 619         __Ensemble,
 620         __Covariance,
 621         __Seed = None,
 622         ):
 623     "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
 624     if hasattr(__Covariance,"assparsematrix"):
 625         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
 626             # Traitement d'une covariance nulle ou presque
 627             return __Ensemble
 628         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
 629             # Traitement d'une covariance nulle ou presque
 630             return __Ensemble
 631     else:
 632         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
 633             # Traitement d'une covariance nulle ou presque
 634             return __Ensemble
 635         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
 636             # Traitement d'une covariance nulle ou presque
 637             return __Ensemble
 638     #
 639     __n, __m = __Ensemble.shape
 640     if __Seed is not None: numpy.random.seed(__Seed)
 641     #
 642     if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
 643         # Traitement d'une covariance multiple de l'identité
 644         __zero = 0.
 645         __std  = numpy.sqrt(__Covariance.assparsematrix())
 646         __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
 647     #
 648     elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
 649         # Traitement d'une covariance diagonale avec variances non identiques
 650         __zero = numpy.zeros(__n)
 651         __std  = numpy.sqrt(__Covariance.assparsematrix())
 652         __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
 653     #
 654     elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
 655         # Traitement d'une covariance pleine
 656         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
 657     #
 658     elif isinstance(__Covariance, numpy.ndarray):
 659         # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
 660         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
 661     #
 662     else:
 663         raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
 664     #
 665     return __Ensemble
 666
 667 # ==============================================================================
 668 def CovarianceInflation(
 669         __InputCovOrEns,
 670         __InflationType   = None,
 671         __InflationFactor = None,
 672         __BackgroundCov   = None,
 673         ):
 674     """
 675     Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
 676
 677     Synthèse : Hunt 2007, section 2.3.5
 678     """
 679     if __InflationFactor is None:
 680         return __InputCovOrEns
 681     else:
 682         __InflationFactor = float(__InflationFactor)
 683     #
 684     __InputCovOrEns = numpy.asarray(__InputCovOrEns)
 685     if __InputCovOrEns.size == 0: return __InputCovOrEns
 686     #
 687     if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
 688         if __InflationFactor < 1.:
 689             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 690         if __InflationFactor < 1.+mpr: # No inflation = 1
 691             return __InputCovOrEns
 692         __OutputCovOrEns = __InflationFactor**2 * __InputCovOrEns
 693     #
 694     elif __InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
 695         if __InflationFactor < 1.:
 696             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 697         if __InflationFactor < 1.+mpr: # No inflation = 1
 698             return __InputCovOrEns
 699         __InputCovOrEnsMean = __InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
 700         __OutputCovOrEns = __InputCovOrEnsMean[:,numpy.newaxis] \
 701             + __InflationFactor * (__InputCovOrEns - __InputCovOrEnsMean[:,numpy.newaxis])
 702     #
 703     elif __InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
 704         if __InflationFactor < 0.:
 705             raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
 706         if __InflationFactor < mpr: # No inflation = 0
 707             return __InputCovOrEns
 708         __n, __m = __InputCovOrEns.shape
 709         if __n != __m:
 710             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 711         __tr = __InputCovOrEns.trace()/__n
 712         if __InflationFactor > __tr:
 713             raise ValueError("Inflation factor for additive inflation has to be small over %.0e."%__tr)
 714         __OutputCovOrEns = (1. - __InflationFactor)*__InputCovOrEns + __InflationFactor * numpy.identity(__n)
 715     #
 716     elif __InflationType == "HybridOnBackgroundCovariance":
 717         if __InflationFactor < 0.:
 718             raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
 719         if __InflationFactor < mpr: # No inflation = 0
 720             return __InputCovOrEns
 721         __n, __m = __InputCovOrEns.shape
 722         if __n != __m:
 723             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 724         if __BackgroundCov is None:
 725             raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
 726         if __InputCovOrEns.shape != __BackgroundCov.shape:
 727             raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
 728         __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov
 729     #
 730     elif __InflationType == "Relaxation":
 731         raise NotImplementedError("Relaxation inflation type not implemented")
 732     #
 733     else:
 734         raise ValueError("Error in inflation type, '%s' is not a valid keyword."%__InflationType)
 735     #
 736     return __OutputCovOrEns
 737
 738 # ==============================================================================
 739 def HessienneEstimation(__selfA, __nb, __HaM, __HtM, __BI, __RI):
 740     "Estimation de la Hessienne"
 741     #
 742     __HessienneI = []
 743     for i in range(int(__nb)):
 744         __ee    = numpy.zeros((__nb,1))
 745         __ee[i] = 1.
 746         __HtEE  = numpy.dot(__HtM,__ee).reshape((-1,1))
 747         __HessienneI.append( numpy.ravel( __BI * __ee + __HaM * (__RI * __HtEE) ) )
 748     #
 749     __A = numpy.linalg.inv(numpy.array( __HessienneI ))
 750     __A = (__A + __A.T) * 0.5 # Symétrie
 751     __A = __A + mpr*numpy.trace( __A ) * numpy.identity(__nb) # Positivité
 752     #
 753     if min(__A.shape) != max(__A.shape):
 754         raise ValueError(
 755             "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
 756             " is of shape %s, despites it has to be a"%(str(__A.shape),)+\
 757             " squared matrix. There is an error in the observation operator,"+\
 758             " please check it.")
 759     if (numpy.diag(__A) < 0).any():
 760         raise ValueError(
 761             "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
 762             " has at least one negative value on its diagonal. There is an"+\
 763             " error in the observation operator, please check it.")
 764     if logging.getLogger().level < logging.WARNING: # La vérification n'a lieu qu'en debug
 765         try:
 766             numpy.linalg.cholesky( __A )
 767         except Exception:
 768             raise ValueError(
 769                 "The %s a posteriori covariance matrix A"%(__selfA._name,)+\
 770                 " is not symmetric positive-definite. Please check your a"+\
 771                 " priori covariances and your observation operator.")
 772     #
 773     return __A
 774
 775 # ==============================================================================
 776 def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
 777     "Estimation des quantiles a posteriori à partir de A>0 (selfA est modifié)"
 778     nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
 779     #
 780     # Traitement des bornes
 781     if "StateBoundsForQuantiles" in selfA._parameters:
 782         LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
 783     elif "Bounds" in selfA._parameters:
 784         LBounds = selfA._parameters["Bounds"]  # Défaut raisonnable
 785     else:
 786         LBounds = None
 787     if LBounds is not None:
 788         LBounds = ForceNumericBounds( LBounds )
 789     __Xa = numpy.ravel(Xa)
 790     #
 791     # Échantillonnage des états
 792     YfQ  = None
 793     EXr  = None
 794     for i in range(nbsamples):
 795         if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None and HXa is not None:
 796             dXr = (numpy.random.multivariate_normal(__Xa,A) - __Xa).reshape((-1,1))
 797             if LBounds is not None: # "EstimateProjection" par défaut
 798                 dXr = numpy.max(numpy.hstack((dXr,LBounds[:,0].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
 799                 dXr = numpy.min(numpy.hstack((dXr,LBounds[:,1].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
 800             dYr = HtM @ dXr
 801             Yr = HXa.reshape((-1,1)) + dYr
 802             if selfA._toStore("SampledStateForQuantiles"): Xr = __Xa + numpy.ravel(dXr)
 803         elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
 804             Xr = numpy.random.multivariate_normal(__Xa,A)
 805             if LBounds is not None: # "EstimateProjection" par défaut
 806                 Xr = numpy.max(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,0].reshape((-1,1)))),axis=1)
 807                 Xr = numpy.min(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,1].reshape((-1,1)))),axis=1)
 808             Yr = numpy.asarray(Hm( Xr ))
 809         else:
 810             raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
 811         #
 812         if YfQ is None:
 813             YfQ = Yr.reshape((-1,1))
 814             if selfA._toStore("SampledStateForQuantiles"): EXr = Xr.reshape((-1,1))
 815         else:
 816             YfQ = numpy.hstack((YfQ,Yr.reshape((-1,1))))
 817             if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.hstack((EXr,Xr.reshape((-1,1))))
 818     #
 819     # Extraction des quantiles
 820     YfQ.sort(axis=-1)
 821     YQ = None
 822     for quantile in selfA._parameters["Quantiles"]:
 823         if not (0. <= float(quantile) <= 1.): continue
 824         indice = int(nbsamples * float(quantile) - 1./nbsamples)
 825         if YQ is None: YQ = YfQ[:,indice].reshape((-1,1))
 826         else:          YQ = numpy.hstack((YQ,YfQ[:,indice].reshape((-1,1))))
 827     if YQ is not None: # Liste non vide de quantiles
 828         selfA.StoredVariables["SimulationQuantiles"].store( YQ )
 829     if selfA._toStore("SampledStateForQuantiles"):
 830         selfA.StoredVariables["SampledStateForQuantiles"].store( EXr )
 831     #
 832     return 0
 833
 834 # ==============================================================================
 835 def ForceNumericBounds( __Bounds, __infNumbers = True ):
 836     "Force les bornes à être des valeurs numériques, sauf si globalement None"
 837     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
 838     if __Bounds is None: return None
 839     # Converti toutes les bornes individuelles None à +/- l'infini chiffré
 840     __Bounds = numpy.asarray( __Bounds, dtype=float )
 841     if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
 842         raise ValueError("Incorrectly shaped bounds data")
 843     if __infNumbers:
 844         __Bounds[numpy.isnan(__Bounds[:,0]),0] = -float('inf')
 845         __Bounds[numpy.isnan(__Bounds[:,1]),1] =  float('inf')
 846     else:
 847         __Bounds[numpy.isnan(__Bounds[:,0]),0] = -sys.float_info.max
 848         __Bounds[numpy.isnan(__Bounds[:,1]),1] =  sys.float_info.max
 849     return __Bounds
 850
 851 # ==============================================================================
 852 def RecentredBounds( __Bounds, __Center, __Scale = None):
 853     "Recentre les bornes autour de 0, sauf si globalement None"
 854     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
 855     if __Bounds is None: return None
 856     if __Scale is None:
 857         # Recentre les valeurs numériques de bornes
 858         return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1))
 859     else:
 860         # Recentre les valeurs numériques de bornes et change l'échelle par une matrice
 861         return __Scale @ (ForceNumericBounds( __Bounds, False ) - numpy.ravel( __Center ).reshape((-1,1)))
 862
 863 # ==============================================================================
 864 def ApplyBounds( __Vector, __Bounds, __newClip = True):
 865     "Applique des bornes numériques à un point"
 866     # Conserve une valeur par défaut s'il n'y a pas de bornes
 867     if __Bounds is None: return __Vector
 868     #
 869     if not isinstance(__Vector, numpy.ndarray): # Is an array
 870         raise ValueError("Incorrect array definition of vector data")
 871     if not isinstance(__Bounds, numpy.ndarray): # Is an array
 872         raise ValueError("Incorrect array definition of bounds data")
 873     if 2*__Vector.size != __Bounds.size: # Is a 2 column array of vector length
 874         raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size,__Vector.size))
 875     if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
 876         raise ValueError("Incorrectly shaped bounds data")
 877     #
 878     if __newClip:
 879         __Vector = __Vector.clip(
 880             __Bounds[:,0].reshape(__Vector.shape),
 881             __Bounds[:,1].reshape(__Vector.shape),
 882             )
 883     else:
 884         __Vector = numpy.max(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,0])),axis=1)
 885         __Vector = numpy.min(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,1])),axis=1)
 886         __Vector = numpy.asarray(__Vector)
 887     #
 888     return __Vector
 889
 890 # ==============================================================================
 891 def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __SuppPars):
 892     "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
 893     __Betaf = __SuppPars["HybridCovarianceEquilibrium"]
 894     #
 895     Xf = EnsembleMean( __EnXf )
 896     Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
 897     Pf = (1 - __Betaf) * __B.asfullmatrix(Xf.size) + __Betaf * Pf
 898     #
 899     selfB = PartialAlgorithm("3DVAR")
 900     selfB._parameters["Minimizer"] = "LBFGSB"
 901     selfB._parameters["MaximumNumberOfIterations"] = __SuppPars["HybridMaximumNumberOfIterations"]
 902     selfB._parameters["CostDecrementTolerance"] = __SuppPars["HybridCostDecrementTolerance"]
 903     selfB._parameters["ProjectedGradientTolerance"] = -1
 904     selfB._parameters["GradientNormTolerance"] = 1.e-05
 905     selfB._parameters["StoreInternalVariables"] = False
 906     selfB._parameters["optiprint"] = -1
 907     selfB._parameters["optdisp"] = 0
 908     selfB._parameters["Bounds"] = None
 909     selfB._parameters["InitializationPoint"] = Xf
 910     from daAlgorithms.Atoms import std3dvar
 911     std3dvar.std3dvar(selfB, Xf, __Ynpu, None, __HO, None, __R, Pf)
 912     Xa = selfB.get("Analysis")[-1].reshape((-1,1))
 913     del selfB
 914     #
 915     return Xa + EnsembleOfAnomalies( __EnXn )
 916
 917 # ==============================================================================
 918 def BuildComplexSampleList(
 919     __SampleAsnUplet,
 920     __SampleAsExplicitHyperCube,
 921     __SampleAsMinMaxStepHyperCube,
 922     __SampleAsIndependantRandomVariables,
 923     __X0,
 924     ):
 925     # ---------------------------
 926     if len(__SampleAsnUplet) > 0:
 927         sampleList = __SampleAsnUplet
 928         for i,Xx in enumerate(sampleList):
 929             if numpy.ravel(Xx).size != __X0.size:
 930                 raise ValueError("The size %i of the %ith state X in the sample and %i of the checking point Xb are different, they have to be identical."%(numpy.ravel(Xx).size,i+1,X0.size))
 931     elif len(__SampleAsExplicitHyperCube) > 0:
 932         sampleList = itertools.product(*list(__SampleAsExplicitHyperCube))
 933     elif len(__SampleAsMinMaxStepHyperCube) > 0:
 934         coordinatesList = []
 935         for i,dim in enumerate(__SampleAsMinMaxStepHyperCube):
 936             if len(dim) != 3:
 937                 raise ValueError("For dimension %i, the variable definition \"%s\" is incorrect, it should be [min,max,step]."%(i,dim))
 938             else:
 939                 coordinatesList.append(numpy.linspace(dim[0],dim[1],1+int((float(dim[1])-float(dim[0]))/float(dim[2]))))
 940         sampleList = itertools.product(*coordinatesList)
 941     elif len(__SampleAsIndependantRandomVariables) > 0:
 942         coordinatesList = []
 943         for i,dim in enumerate(__SampleAsIndependantRandomVariables):
 944             if len(dim) != 3:
 945                 raise ValueError("For dimension %i, the variable definition \"%s\" is incorrect, it should be ('distribution',(parameters),length) with distribution in ['normal'(mean,std),'lognormal'(mean,sigma),'uniform'(low,high),'weibull'(shape)]."%(i,dim))
 946             elif not( str(dim[0]) in ['normal','lognormal','uniform','weibull'] and hasattr(numpy.random,dim[0]) ):
 947                 raise ValueError("For dimension %i, the distribution name \"%s\" is not allowed, please choose in ['normal'(mean,std),'lognormal'(mean,sigma),'uniform'(low,high),'weibull'(shape)]"%(i,dim[0]))
 948             else:
 949                 distribution = getattr(numpy.random,str(dim[0]),'normal')
 950                 coordinatesList.append(distribution(*dim[1], size=max(1,int(dim[2]))))
 951         sampleList = itertools.product(*coordinatesList)
 952     else:
 953         sampleList = iter([__X0,])
 954     # ----------
 955     return sampleList
 956
 957 # ==============================================================================
 958 def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle,
 959         __CovForecast = False, __LinEvolution = False,
 960         ):
 961     """
 962     Prévision multi-pas avec une correction par pas (multi-méthodes)
 963     """
 964     #
 965     # Initialisation
 966     # --------------
 967     if selfA._parameters["EstimationOf"] == "State":
 968         if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
 969             Xn = numpy.asarray(Xb)
 970             if __CovForecast: Pn = B
 971             selfA.StoredVariables["Analysis"].store( Xn )
 972             if selfA._toStore("APosterioriCovariance"):
 973                 if hasattr(B,"asfullmatrix"):
 974                     selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
 975                 else:
 976                     selfA.StoredVariables["APosterioriCovariance"].store( B )
 977             selfA._setInternalState("seed", numpy.random.get_state())
 978         elif selfA._parameters["nextStep"]:
 979             Xn = selfA._getInternalState("Xn")
 980             if __CovForecast: Pn = selfA._getInternalState("Pn")
 981     else:
 982         Xn = numpy.asarray(Xb)
 983         if __CovForecast: Pn = B
 984     #
 985     if hasattr(Y,"stepnumber"):
 986         duration = Y.stepnumber()
 987     else:
 988         duration = 2
 989     #
 990     # Multi-steps
 991     # -----------
 992     for step in range(duration-1):
 993         selfA.StoredVariables["CurrentStepNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 994         #
 995         if hasattr(Y,"store"):
 996             Ynpu = numpy.asarray( Y[step+1] ).reshape((-1,1))
 997         else:
 998             Ynpu = numpy.asarray( Y ).reshape((-1,1))
 999         #
1000         if U is not None:
1001             if hasattr(U,"store") and len(U)>1:
1002                 Un = numpy.asarray( U[step] ).reshape((-1,1))
1003             elif hasattr(U,"store") and len(U)==1:
1004                 Un = numpy.asarray( U[0] ).reshape((-1,1))
1005             else:
1006                 Un = numpy.asarray( U ).reshape((-1,1))
1007         else:
1008             Un = None
1009         #
1010         # Predict (Time Update)
1011         # ---------------------
1012         if selfA._parameters["EstimationOf"] == "State":
1013             if __CovForecast or __LinEvolution:
1014                 Mt = EM["Tangent"].asMatrix(Xn)
1015                 Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
1016             if __CovForecast:
1017                 Ma = EM["Adjoint"].asMatrix(Xn)
1018                 Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
1019                 Pn_predicted = Q + Mt @ (Pn @ Ma)
1020             if __LinEvolution:
1021                 Xn_predicted = Mt @ Xn
1022             else:
1023                 M  = EM["Direct"].appliedControledFormTo
1024                 Xn_predicted = M( (Xn, Un) )
1025             if CM is not None and "Tangent" in CM and Un is not None: # Attention : si Cm est aussi dans M, doublon !
1026                 Cm = CM["Tangent"].asMatrix(Xn_predicted)
1027                 Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
1028                 Xn_predicted = Xn_predicted + (Cm @ Un).reshape((-1,1))
1029         elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
1030             # --- > Par principe, M = Id, Q = 0
1031             Xn_predicted = Xn
1032             if __CovForecast: Pn_predicted = Pn
1033         Xn_predicted = numpy.asarray(Xn_predicted).reshape((-1,1))
1034         if selfA._toStore("ForecastState"):
1035             selfA.StoredVariables["ForecastState"].store( Xn_predicted )
1036         if __CovForecast:
1037             if hasattr(Pn_predicted,"asfullmatrix"):
1038                 Pn_predicted = Pn_predicted.asfullmatrix(Xn.size)
1039             else:
1040                 Pn_predicted = numpy.asarray(Pn_predicted).reshape((Xn.size,Xn.size))
1041             if selfA._toStore("ForecastCovariance"):
1042                 selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
1043         #
1044         # Correct (Measurement Update)
1045         # ----------------------------
1046         if __CovForecast:
1047             oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, Pn_predicted, True)
1048         else:
1049             oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, B, True)
1050         #
1051         #--------------------------
1052         Xn = selfA._getInternalState("Xn")
1053         if __CovForecast: Pn = selfA._getInternalState("Pn")
1054     #
1055     return 0
1056
1057 # ==============================================================================
1058 if __name__ == "__main__":
1059     print('\n AUTODIAGNOSTIC\n')