src/daComposant/daCore/NumericObjects.py

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