src/daComposant/daCore/NumericObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2021 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Définit les objets numériques génériques.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, time, copy, types, sys, logging
  29 import math, numpy, scipy, scipy.optimize, scipy.version
  30 from daCore.BasicObjects import Operator
  31 from daCore.PlatformInfo import PlatformInfo
  32 mpr = PlatformInfo().MachinePrecision()
  33 mfp = PlatformInfo().MaximumPrecision()
  34 # logging.getLogger().setLevel(logging.DEBUG)
  35
  36 # ==============================================================================
  37 def ExecuteFunction( triplet ):
  38     assert len(triplet) == 3, "Incorrect number of arguments"
  39     X, xArgs, funcrepr = triplet
  40     __X = numpy.asmatrix(numpy.ravel( X )).T
  41     __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
  42     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  43     __fonction = getattr(__module,funcrepr["__userFunction__name"])
  44     sys.path = __sys_path_tmp ; del __sys_path_tmp
  45     if isinstance(xArgs, dict):
  46         __HX  = __fonction( __X, **xArgs )
  47     else:
  48         __HX  = __fonction( __X )
  49     return numpy.ravel( __HX )
  50
  51 # ==============================================================================
  52 class FDApproximation(object):
  53     """
  54     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  55     création d'un objet, en fournissant une fonction "Function", on obtient un
  56     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  57     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  58     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  59     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  60     centrées si le booléen "centeredDF" est vrai.
  61     """
  62     def __init__(self,
  63             name                  = "FDApproximation",
  64             Function              = None,
  65             centeredDF            = False,
  66             increment             = 0.01,
  67             dX                    = None,
  68             extraArguments        = None,
  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         if mpEnabled:
  79             try:
  80                 import multiprocessing
  81                 self.__mpEnabled = True
  82             except ImportError:
  83                 self.__mpEnabled = False
  84         else:
  85             self.__mpEnabled = False
  86         self.__mpWorkers = mpWorkers
  87         if self.__mpWorkers is not None and self.__mpWorkers < 1:
  88             self.__mpWorkers = None
  89         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
  90         #
  91         if mfEnabled:
  92             self.__mfEnabled = True
  93         else:
  94             self.__mfEnabled = False
  95         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
  96         #
  97         if avoidingRedundancy:
  98             self.__avoidRC = True
  99             self.__tolerBP = float(toleranceInRedundancy)
 100             self.__lenghtRJ = int(lenghtOfRedundancy)
 101             self.__listJPCP = [] # Jacobian Previous Calculated Points
 102             self.__listJPCI = [] # Jacobian Previous Calculated Increment
 103             self.__listJPCR = [] # Jacobian Previous Calculated Results
 104             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
 105             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
 106         else:
 107             self.__avoidRC = False
 108         #
 109         if self.__mpEnabled:
 110             if isinstance(Function,types.FunctionType):
 111                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 112                 self.__userFunction__name = Function.__name__
 113                 try:
 114                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 115                 except:
 116                     mod = os.path.abspath(Function.__globals__['__file__'])
 117                 if not os.path.isfile(mod):
 118                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 119                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 120                 self.__userFunction__path = os.path.dirname(mod)
 121                 del mod
 122                 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
 123                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 124             elif isinstance(Function,types.MethodType):
 125                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 126                 self.__userFunction__name = Function.__name__
 127                 try:
 128                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 129                 except:
 130                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 131                 if not os.path.isfile(mod):
 132                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 133                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 134                 self.__userFunction__path = os.path.dirname(mod)
 135                 del mod
 136                 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
 137                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 138             else:
 139                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 140         else:
 141             self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
 142             self.__userFunction = self.__userOperator.appliedTo
 143         #
 144         self.__centeredDF = bool(centeredDF)
 145         if abs(float(increment)) > 1.e-15:
 146             self.__increment  = float(increment)
 147         else:
 148             self.__increment  = 0.01
 149         if dX is None:
 150             self.__dX     = None
 151         else:
 152             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
 153         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
 154         if self.__avoidRC:
 155             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
 156
 157     # ---------------------------------------------------------
 158     def __doublon__(self, e, l, n, v=None):
 159         __ac, __iac = False, -1
 160         for i in range(len(l)-1,-1,-1):
 161             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
 162                 __ac, __iac = True, i
 163                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
 164                 break
 165         return __ac, __iac
 166
 167     # ---------------------------------------------------------
 168     def DirectOperator(self, X, **extraArgs ):
 169         """
 170         Calcul du direct à l'aide de la fonction fournie.
 171
 172         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 173         ne doivent pas être données ici à la fonction utilisateur.
 174         """
 175         logging.debug("FDA Calcul DirectOperator (explicite)")
 176         if self.__mfEnabled:
 177             _HX = self.__userFunction( X, argsAsSerie = True )
 178         else:
 179             _X = numpy.asmatrix(numpy.ravel( X )).T
 180             _HX = numpy.ravel(self.__userFunction( _X ))
 181         #
 182         return _HX
 183
 184     # ---------------------------------------------------------
 185     def TangentMatrix(self, X ):
 186         """
 187         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 188         c'est-à-dire le gradient de H en X. On utilise des différences finies
 189         directionnelles autour du point X. X est un numpy.matrix.
 190
 191         Différences finies centrées (approximation d'ordre 2):
 192         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 193            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 194            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 195            H( X_moins_dXi )
 196         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 197            le pas 2*dXi
 198         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 199
 200         Différences finies non centrées (approximation d'ordre 1):
 201         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 202            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 203            HX_plus_dXi = H( X_plus_dXi )
 204         2/ On calcule la valeur centrale HX = H(X)
 205         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 206            le pas dXi
 207         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 208
 209         """
 210         logging.debug("FDA Début du calcul de la Jacobienne")
 211         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 212         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 213         #
 214         if X is None or len(X)==0:
 215             raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
 216         #
 217         _X = numpy.asmatrix(numpy.ravel( X )).T
 218         #
 219         if self.__dX is None:
 220             _dX  = self.__increment * _X
 221         else:
 222             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 223         #
 224         if (_dX == 0.).any():
 225             moyenne = _dX.mean()
 226             if moyenne == 0.:
 227                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 228             else:
 229                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 230         #
 231         __alreadyCalculated  = False
 232         if self.__avoidRC:
 233             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 234             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 235             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 236                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 237                 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
 238         #
 239         if __alreadyCalculated:
 240             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 241             _Jacobienne = self.__listJPCR[__i]
 242         else:
 243             logging.debug("FDA   Calcul Jacobienne (explicite)")
 244             if self.__centeredDF:
 245                 #
 246                 if self.__mpEnabled and not self.__mfEnabled:
 247                     funcrepr = {
 248                         "__userFunction__path" : self.__userFunction__path,
 249                         "__userFunction__modl" : self.__userFunction__modl,
 250                         "__userFunction__name" : self.__userFunction__name,
 251                     }
 252                     _jobs = []
 253                     for i in range( len(_dX) ):
 254                         _dXi            = _dX[i]
 255                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 256                         _X_plus_dXi[i]  = _X[i] + _dXi
 257                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 258                         _X_moins_dXi[i] = _X[i] - _dXi
 259                         #
 260                         _jobs.append( (_X_plus_dXi,  self.__extraArgs, funcrepr) )
 261                         _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
 262                     #
 263                     import multiprocessing
 264                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 265                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 266                     self.__pool.close()
 267                     self.__pool.join()
 268                     #
 269                     _Jacobienne  = []
 270                     for i in range( len(_dX) ):
 271                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 272                     #
 273                 elif self.__mfEnabled:
 274                     _xserie = []
 275                     for i in range( len(_dX) ):
 276                         _dXi            = _dX[i]
 277                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 278                         _X_plus_dXi[i]  = _X[i] + _dXi
 279                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 280                         _X_moins_dXi[i] = _X[i] - _dXi
 281                         #
 282                         _xserie.append( _X_plus_dXi )
 283                         _xserie.append( _X_moins_dXi )
 284                     #
 285                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 286                      #
 287                     _Jacobienne  = []
 288                     for i in range( len(_dX) ):
 289                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 290                     #
 291                 else:
 292                     _Jacobienne  = []
 293                     for i in range( _dX.size ):
 294                         _dXi            = _dX[i]
 295                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 296                         _X_plus_dXi[i]  = _X[i] + _dXi
 297                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 298                         _X_moins_dXi[i] = _X[i] - _dXi
 299                         #
 300                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 301                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 302                         #
 303                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 304                 #
 305             else:
 306                 #
 307                 if self.__mpEnabled and not self.__mfEnabled:
 308                     funcrepr = {
 309                         "__userFunction__path" : self.__userFunction__path,
 310                         "__userFunction__modl" : self.__userFunction__modl,
 311                         "__userFunction__name" : self.__userFunction__name,
 312                     }
 313                     _jobs = []
 314                     _jobs.append( (_X.A1, self.__extraArgs, funcrepr) )
 315                     for i in range( len(_dX) ):
 316                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 317                         _X_plus_dXi[i] = _X[i] + _dX[i]
 318                         #
 319                         _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
 320                     #
 321                     import multiprocessing
 322                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 323                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 324                     self.__pool.close()
 325                     self.__pool.join()
 326                     #
 327                     _HX = _HX_plus_dX.pop(0)
 328                     #
 329                     _Jacobienne = []
 330                     for i in range( len(_dX) ):
 331                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 332                     #
 333                 elif self.__mfEnabled:
 334                     _xserie = []
 335                     _xserie.append( _X.A1 )
 336                     for i in range( len(_dX) ):
 337                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 338                         _X_plus_dXi[i] = _X[i] + _dX[i]
 339                         #
 340                         _xserie.append( _X_plus_dXi )
 341                     #
 342                     _HX_plus_dX = self.DirectOperator( _xserie )
 343                     #
 344                     _HX = _HX_plus_dX.pop(0)
 345                     #
 346                     _Jacobienne = []
 347                     for i in range( len(_dX) ):
 348                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 349                    #
 350                 else:
 351                     _Jacobienne  = []
 352                     _HX = self.DirectOperator( _X )
 353                     for i in range( _dX.size ):
 354                         _dXi            = _dX[i]
 355                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 356                         _X_plus_dXi[i]  = _X[i] + _dXi
 357                         #
 358                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 359                         #
 360                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 361                 #
 362             #
 363             _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
 364             if self.__avoidRC:
 365                 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
 366                 while len(self.__listJPCP) > self.__lenghtRJ:
 367                     self.__listJPCP.pop(0)
 368                     self.__listJPCI.pop(0)
 369                     self.__listJPCR.pop(0)
 370                     self.__listJPPN.pop(0)
 371                     self.__listJPIN.pop(0)
 372                 self.__listJPCP.append( copy.copy(_X) )
 373                 self.__listJPCI.append( copy.copy(_dX) )
 374                 self.__listJPCR.append( copy.copy(_Jacobienne) )
 375                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 376                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 377         #
 378         logging.debug("FDA Fin du calcul de la Jacobienne")
 379         #
 380         return _Jacobienne
 381
 382     # ---------------------------------------------------------
 383     def TangentOperator(self, paire, **extraArgs ):
 384         """
 385         Calcul du tangent à l'aide de la Jacobienne.
 386
 387         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 388         ne doivent pas être données ici à la fonction utilisateur.
 389         """
 390         if self.__mfEnabled:
 391             assert len(paire) == 1, "Incorrect lenght of arguments"
 392             _paire = paire[0]
 393             assert len(_paire) == 2, "Incorrect number of arguments"
 394         else:
 395             assert len(paire) == 2, "Incorrect number of arguments"
 396             _paire = paire
 397         X, dX = _paire
 398         _Jacobienne = self.TangentMatrix( X )
 399         if dX is None or len(dX) == 0:
 400             #
 401             # Calcul de la forme matricielle si le second argument est None
 402             # -------------------------------------------------------------
 403             if self.__mfEnabled: return [_Jacobienne,]
 404             else:                return _Jacobienne
 405         else:
 406             #
 407             # Calcul de la valeur linéarisée de H en X appliqué à dX
 408             # ------------------------------------------------------
 409             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 410             _HtX = numpy.dot(_Jacobienne, _dX)
 411             if self.__mfEnabled: return [_HtX.A1,]
 412             else:                return _HtX.A1
 413
 414     # ---------------------------------------------------------
 415     def AdjointOperator(self, paire, **extraArgs ):
 416         """
 417         Calcul de l'adjoint à l'aide de la Jacobienne.
 418
 419         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 420         ne doivent pas être données ici à la fonction utilisateur.
 421         """
 422         if self.__mfEnabled:
 423             assert len(paire) == 1, "Incorrect lenght of arguments"
 424             _paire = paire[0]
 425             assert len(_paire) == 2, "Incorrect number of arguments"
 426         else:
 427             assert len(paire) == 2, "Incorrect number of arguments"
 428             _paire = paire
 429         X, Y = _paire
 430         _JacobienneT = self.TangentMatrix( X ).T
 431         if Y is None or len(Y) == 0:
 432             #
 433             # Calcul de la forme matricielle si le second argument est None
 434             # -------------------------------------------------------------
 435             if self.__mfEnabled: return [_JacobienneT,]
 436             else:                return _JacobienneT
 437         else:
 438             #
 439             # Calcul de la valeur de l'adjoint en X appliqué à Y
 440             # --------------------------------------------------
 441             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 442             _HaY = numpy.dot(_JacobienneT, _Y)
 443             if self.__mfEnabled: return [_HaY.A1,]
 444             else:                return _HaY.A1
 445
 446 # ==============================================================================
 447 def EnsembleOfCenteredPerturbations( _bgcenter, _bgcovariance, _nbmembers ):
 448     "Génération d'un ensemble de taille _nbmembers-1 d'états aléatoires centrés"
 449     #
 450     _bgcenter = numpy.ravel(_bgcenter)[:,None]
 451     if _nbmembers < 1:
 452         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
 453     #
 454     if _bgcovariance is None:
 455         BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
 456     else:
 457         _Z = numpy.random.multivariate_normal(numpy.zeros(_bgcenter.size), _bgcovariance, size=_nbmembers).T
 458         BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Z
 459     #
 460     return BackgroundEnsemble
 461
 462 # ==============================================================================
 463 def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
 464     "Génération d'un ensemble de taille _nbmembers-1 d'états aléatoires centrés"
 465     def __CenteredRandomAnomalies(Zr, N):
 466         """
 467         Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
 468         notes manuscrites de MB et conforme au code de PS avec eps = -1
 469         """
 470         eps = -1
 471         Q = numpy.identity(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
 472         Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
 473         R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
 474         Q = numpy.dot(Q,R)
 475         Zr = numpy.dot(Q,Zr)
 476         return Zr.T
 477     #
 478     _bgcenter = numpy.ravel(_bgcenter).reshape((-1,1))
 479     if _nbmembers < 1:
 480         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
 481     if _bgcovariance is None:
 482         BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
 483     else:
 484         if _withSVD:
 485             U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
 486             _nbctl = _bgcenter.size
 487             if _nbmembers > _nbctl:
 488                 _Z = numpy.concatenate((numpy.dot(
 489                     numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
 490                     numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1-_nbctl)), axis = 0)
 491             else:
 492                 _Z = numpy.dot(numpy.diag(numpy.sqrt(s[:_nbmembers-1])), V[:_nbmembers-1])
 493             _Zca = __CenteredRandomAnomalies(_Z, _nbmembers)
 494             BackgroundEnsemble = _bgcenter + _Zca
 495         else:
 496             if max(abs(_bgcovariance.flatten())) > 0:
 497                 _nbctl = _bgcenter.size
 498                 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
 499                 _Zca = __CenteredRandomAnomalies(_Z, _nbmembers)
 500                 BackgroundEnsemble = _bgcenter + _Zca
 501             else:
 502                 BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers)
 503     #
 504     return BackgroundEnsemble
 505
 506 # ==============================================================================
 507 def EnsembleMean( __Ensemble ):
 508     "Renvoie la moyenne empirique d'un ensemble"
 509     return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
 510
 511 # ==============================================================================
 512 def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.):
 513     "Renvoie les anomalies centrées à partir d'un ensemble"
 514     if __OptMean is None:
 515         __Em = EnsembleMean( __Ensemble )
 516     else:
 517         __Em = numpy.ravel( __OptMean ).reshape((-1,1))
 518     #
 519     return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
 520
 521 # ==============================================================================
 522 def EnsembleErrorCovariance( __Ensemble, __quick = False ):
 523     "Renvoie l'estimation empirique de la covariance d'ensemble"
 524     if __quick:
 525         # Covariance rapide mais rarement définie positive
 526         __Covariance = numpy.cov( __Ensemble )
 527     else:
 528         # Résultat souvent identique à numpy.cov, mais plus robuste
 529         __n, __m = numpy.asarray( __Ensemble ).shape
 530         __Anomalies = EnsembleOfAnomalies( __Ensemble )
 531         # Estimation empirique
 532         __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m-1)
 533         # Assure la symétrie
 534         __Covariance = ( __Covariance + __Covariance.T ) * 0.5
 535         # Assure la positivité
 536         __epsilon    = mpr*numpy.trace( __Covariance )
 537         __Covariance = __Covariance + __epsilon * numpy.identity(__n)
 538     #
 539     return __Covariance
 540
 541 # ==============================================================================
 542 def EnsemblePerturbationWithGivenCovariance( __Ensemble, __Covariance, __Seed=None ):
 543     "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
 544     if hasattr(__Covariance,"assparsematrix"):
 545         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
 546             # Traitement d'une covariance nulle ou presque
 547             return __Ensemble
 548         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
 549             # Traitement d'une covariance nulle ou presque
 550             return __Ensemble
 551     else:
 552         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
 553             # Traitement d'une covariance nulle ou presque
 554             return __Ensemble
 555         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
 556             # Traitement d'une covariance nulle ou presque
 557             return __Ensemble
 558     #
 559     __n, __m = __Ensemble.shape
 560     if __Seed is not None: numpy.random.seed(__Seed)
 561     #
 562     if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
 563         # Traitement d'une covariance multiple de l'identité
 564         __zero = 0.
 565         __std  = numpy.sqrt(__Covariance.assparsematrix())
 566         __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
 567     #
 568     elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
 569         # Traitement d'une covariance diagonale avec variances non identiques
 570         __zero = numpy.zeros(__n)
 571         __std  = numpy.sqrt(__Covariance.assparsematrix())
 572         __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
 573     #
 574     elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
 575         # Traitement d'une covariance pleine
 576         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
 577     #
 578     elif isinstance(__Covariance, numpy.ndarray):
 579         # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
 580         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
 581     #
 582     else:
 583         raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
 584     #
 585     return __Ensemble
 586
 587 # ==============================================================================
 588 def CovarianceInflation(
 589         InputCovOrEns,
 590         InflationType   = None,
 591         InflationFactor = None,
 592         BackgroundCov   = None,
 593         ):
 594     """
 595     Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
 596
 597     Synthèse : Hunt 2007, section 2.3.5
 598     """
 599     if InflationFactor is None:
 600         return InputCovOrEns
 601     else:
 602         InflationFactor = float(InflationFactor)
 603     #
 604     if InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
 605         if InflationFactor < 1.:
 606             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 607         if InflationFactor < 1.+mpr:
 608             return InputCovOrEns
 609         OutputCovOrEns = InflationFactor**2 * InputCovOrEns
 610     #
 611     elif InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
 612         if InflationFactor < 1.:
 613             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 614         if InflationFactor < 1.+mpr:
 615             return InputCovOrEns
 616         InputCovOrEnsMean = InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
 617         OutputCovOrEns = InputCovOrEnsMean[:,numpy.newaxis] \
 618             + InflationFactor * (InputCovOrEns - InputCovOrEnsMean[:,numpy.newaxis])
 619     #
 620     elif InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
 621         if InflationFactor < 0.:
 622             raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
 623         if InflationFactor < mpr:
 624             return InputCovOrEns
 625         __n, __m = numpy.asarray(InputCovOrEns).shape
 626         if __n != __m:
 627             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 628         OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * numpy.identity(__n)
 629     #
 630     elif InflationType == "HybridOnBackgroundCovariance":
 631         if InflationFactor < 0.:
 632             raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
 633         if InflationFactor < mpr:
 634             return InputCovOrEns
 635         __n, __m = numpy.asarray(InputCovOrEns).shape
 636         if __n != __m:
 637             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 638         if BackgroundCov is None:
 639             raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
 640         if InputCovOrEns.shape != BackgroundCov.shape:
 641             raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
 642         OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * BackgroundCov
 643     #
 644     elif InflationType == "Relaxation":
 645         raise NotImplementedError("InflationType Relaxation")
 646     #
 647     else:
 648         raise ValueError("Error in inflation type, '%s' is not a valid keyword."%InflationType)
 649     #
 650     return OutputCovOrEns
 651
 652 # ==============================================================================
 653 def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
 654     "Estimation des quantiles a posteriori (selfA est modifié)"
 655     nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
 656     #
 657     # Traitement des bornes
 658     if "StateBoundsForQuantiles" in selfA._parameters:
 659         LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
 660     elif "Bounds" in selfA._parameters:
 661         LBounds = selfA._parameters["Bounds"]  # Défaut raisonnable
 662     else:
 663         LBounds = None
 664     if LBounds is not None:
 665         LBounds = numpy.asmatrix(ForceNumericBounds( LBounds ))
 666     #
 667     # Échantillonnage des états
 668     YfQ  = None
 669     EXr  = None
 670     if selfA._parameters["SimulationForQuantiles"] == "Linear" and HXa is not None:
 671         HXa  = numpy.matrix(numpy.ravel( HXa )).T
 672     for i in range(nbsamples):
 673         if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None:
 674             dXr = numpy.matrix(numpy.random.multivariate_normal(numpy.ravel(Xa),A) - numpy.ravel(Xa)).T
 675             if LBounds is not None: # "EstimateProjection" par défaut
 676                 dXr = numpy.max(numpy.hstack((dXr,LBounds[:,0]) - Xa),axis=1)
 677                 dXr = numpy.min(numpy.hstack((dXr,LBounds[:,1]) - Xa),axis=1)
 678             dYr = numpy.matrix(numpy.ravel( HtM * dXr )).T
 679             Yr = HXa + dYr
 680             if selfA._toStore("SampledStateForQuantiles"): Xr = Xa + dXr
 681         elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
 682             Xr = numpy.matrix(numpy.random.multivariate_normal(numpy.ravel(Xa),A)).T
 683             if LBounds is not None: # "EstimateProjection" par défaut
 684                 Xr = numpy.max(numpy.hstack((Xr,LBounds[:,0])),axis=1)
 685                 Xr = numpy.min(numpy.hstack((Xr,LBounds[:,1])),axis=1)
 686             Yr = numpy.matrix(numpy.ravel( Hm( Xr ) )).T
 687         else:
 688             raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
 689         #
 690         if YfQ is None:
 691             YfQ = Yr
 692             if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.ravel(Xr)
 693         else:
 694             YfQ = numpy.hstack((YfQ,Yr))
 695             if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.vstack((EXr,numpy.ravel(Xr)))
 696     #
 697     # Extraction des quantiles
 698     YfQ.sort(axis=-1)
 699     YQ = None
 700     for quantile in selfA._parameters["Quantiles"]:
 701         if not (0. <= float(quantile) <= 1.): continue
 702         indice = int(nbsamples * float(quantile) - 1./nbsamples)
 703         if YQ is None: YQ = YfQ[:,indice]
 704         else:          YQ = numpy.hstack((YQ,YfQ[:,indice]))
 705     selfA.StoredVariables["SimulationQuantiles"].store( YQ )
 706     if selfA._toStore("SampledStateForQuantiles"):
 707         selfA.StoredVariables["SampledStateForQuantiles"].store( EXr.T )
 708     #
 709     return 0
 710
 711 # ==============================================================================
 712 def ForceNumericBounds( __Bounds ):
 713     "Force les bornes à être des valeurs numériques, sauf si globalement None"
 714     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
 715     if __Bounds is None: return None
 716     # Converti toutes les bornes individuelles None à +/- l'infini
 717     __Bounds = numpy.asarray( __Bounds, dtype=float )
 718     if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
 719         raise ValueError("Incorrectly defined bounds data")
 720     __Bounds[numpy.isnan(__Bounds)[:,0],0] = -sys.float_info.max
 721     __Bounds[numpy.isnan(__Bounds)[:,1],1] =  sys.float_info.max
 722     return __Bounds
 723
 724 # ==============================================================================
 725 def RecentredBounds( __Bounds, __Center):
 726     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
 727     if __Bounds is None: return None
 728     # Recentre les valeurs numériques de bornes
 729     return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).transpose((-1,1))
 730
 731 # ==============================================================================
 732 def cekf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
 733     """
 734     Contrained Extended Kalman Filter
 735     """
 736     if selfA._parameters["EstimationOf"] == "Parameters":
 737         selfA._parameters["StoreInternalVariables"] = True
 738     selfA._parameters["Bounds"] = ForceNumericBounds( selfA._parameters["Bounds"] )
 739     #
 740     # Opérateurs
 741     H = HO["Direct"].appliedControledFormTo
 742     #
 743     if selfA._parameters["EstimationOf"] == "State":
 744         M = EM["Direct"].appliedControledFormTo
 745     #
 746     if CM is not None and "Tangent" in CM and U is not None:
 747         Cm = CM["Tangent"].asMatrix(Xb)
 748     else:
 749         Cm = None
 750     #
 751     # Durée d'observation et tailles
 752     if hasattr(Y,"stepnumber"):
 753         duration = Y.stepnumber()
 754         __p = numpy.cumprod(Y.shape())[-1]
 755     else:
 756         duration = 2
 757         __p = numpy.array(Y).size
 758     #
 759     # Précalcul des inversions de B et R
 760     if selfA._parameters["StoreInternalVariables"] \
 761         or selfA._toStore("CostFunctionJ") \
 762         or selfA._toStore("CostFunctionJb") \
 763         or selfA._toStore("CostFunctionJo") \
 764         or selfA._toStore("CurrentOptimum") \
 765         or selfA._toStore("APosterioriCovariance"):
 766         BI = B.getI()
 767         RI = R.getI()
 768     #
 769     __n = Xb.size
 770     #
 771     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
 772         Xn = Xb
 773         Pn = B
 774         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 775         selfA.StoredVariables["Analysis"].store( Xb )
 776         if selfA._toStore("APosterioriCovariance"):
 777             if hasattr(B,"asfullmatrix"):
 778                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
 779             else:
 780                 selfA.StoredVariables["APosterioriCovariance"].store( B )
 781         selfA._setInternalState("seed", numpy.random.get_state())
 782     elif selfA._parameters["nextStep"]:
 783         Xn = selfA._getInternalState("Xn")
 784         Pn = selfA._getInternalState("Pn")
 785     #
 786     if selfA._parameters["EstimationOf"] == "Parameters":
 787         XaMin            = Xn
 788         previousJMinimum = numpy.finfo(float).max
 789     #
 790     for step in range(duration-1):
 791         if hasattr(Y,"store"):
 792             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
 793         else:
 794             Ynpu = numpy.ravel( Y ).reshape((__p,1))
 795         #
 796         Ht = HO["Tangent"].asMatrix(ValueForMethodForm = Xn)
 797         Ht = Ht.reshape(Ynpu.size,Xn.size) # ADAO & check shape
 798         Ha = HO["Adjoint"].asMatrix(ValueForMethodForm = Xn)
 799         Ha = Ha.reshape(Xn.size,Ynpu.size) # ADAO & check shape
 800         #
 801         if selfA._parameters["EstimationOf"] == "State":
 802             Mt = EM["Tangent"].asMatrix(ValueForMethodForm = Xn)
 803             Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
 804             Ma = EM["Adjoint"].asMatrix(ValueForMethodForm = Xn)
 805             Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
 806         #
 807         if U is not None:
 808             if hasattr(U,"store") and len(U)>1:
 809                 Un = numpy.ravel( U[step] ).reshape((-1,1))
 810             elif hasattr(U,"store") and len(U)==1:
 811                 Un = numpy.ravel( U[0] ).reshape((-1,1))
 812             else:
 813                 Un = numpy.ravel( U ).reshape((-1,1))
 814         else:
 815             Un = None
 816         #
 817         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
 818             Xn = numpy.max(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
 819             Xn = numpy.min(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
 820         #
 821         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
 822             Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).reshape((__n,1))
 823             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
 824                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
 825                 Xn_predicted = Xn_predicted + Cm * Un
 826             Pn_predicted = Q + Mt * Pn * Ma
 827         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
 828             # --- > Par principe, M = Id, Q = 0
 829             Xn_predicted = Xn
 830             Pn_predicted = Pn
 831         #
 832         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
 833             Xn_predicted = numpy.max(numpy.hstack((Xn_predicted,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
 834             Xn_predicted = numpy.min(numpy.hstack((Xn_predicted,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
 835         #
 836         if selfA._parameters["EstimationOf"] == "State":
 837             HX_predicted = numpy.ravel( H( (Xn_predicted, None) ) ).reshape((__p,1))
 838             _Innovation  = Ynpu - HX_predicted
 839         elif selfA._parameters["EstimationOf"] == "Parameters":
 840             HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1))
 841             _Innovation  = Ynpu - HX_predicted
 842             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
 843                 _Innovation = _Innovation - Cm * Un
 844         #
 845         Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
 846         Xn = Xn_predicted + Kn * _Innovation
 847         Pn = Pn_predicted - Kn * Ht * Pn_predicted
 848         #
 849         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
 850             Xn = numpy.max(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
 851             Xn = numpy.min(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
 852         #
 853         Xa = Xn # Pointeurs
 854         #--------------------------
 855         selfA._setInternalState("Xn", Xn)
 856         selfA._setInternalState("Pn", Pn)
 857         #--------------------------
 858         #
 859         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 860         # ---> avec analysis
 861         selfA.StoredVariables["Analysis"].store( Xa )
 862         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
 863             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( H((Xa, Un)) )
 864         if selfA._toStore("InnovationAtCurrentAnalysis"):
 865             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
 866         # ---> avec current state
 867         if selfA._parameters["StoreInternalVariables"] \
 868             or selfA._toStore("CurrentState"):
 869             selfA.StoredVariables["CurrentState"].store( Xn )
 870         if selfA._toStore("ForecastState"):
 871             selfA.StoredVariables["ForecastState"].store( Xn_predicted )
 872         if selfA._toStore("ForecastCovariance"):
 873             selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
 874         if selfA._toStore("BMA"):
 875             selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
 876         if selfA._toStore("InnovationAtCurrentState"):
 877             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
 878         if selfA._toStore("SimulatedObservationAtCurrentState") \
 879             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
 880             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
 881         # ---> autres
 882         if selfA._parameters["StoreInternalVariables"] \
 883             or selfA._toStore("CostFunctionJ") \
 884             or selfA._toStore("CostFunctionJb") \
 885             or selfA._toStore("CostFunctionJo") \
 886             or selfA._toStore("CurrentOptimum") \
 887             or selfA._toStore("APosterioriCovariance"):
 888             Jb  = float( 0.5 * (Xa - Xb).T @ (BI @ (Xa - Xb)) )
 889             Jo  = float( 0.5 * _Innovation.T @ (RI @ _Innovation) )
 890             J   = Jb + Jo
 891             selfA.StoredVariables["CostFunctionJb"].store( Jb )
 892             selfA.StoredVariables["CostFunctionJo"].store( Jo )
 893             selfA.StoredVariables["CostFunctionJ" ].store( J )
 894             #
 895             if selfA._toStore("IndexOfOptimum") \
 896                 or selfA._toStore("CurrentOptimum") \
 897                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
 898                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
 899                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
 900                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
 901                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
 902             if selfA._toStore("IndexOfOptimum"):
 903                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
 904             if selfA._toStore("CurrentOptimum"):
 905                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
 906             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
 907                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
 908             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
 909                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
 910             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
 911                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
 912             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
 913                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
 914         if selfA._toStore("APosterioriCovariance"):
 915             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
 916         if selfA._parameters["EstimationOf"] == "Parameters" \
 917             and J < previousJMinimum:
 918             previousJMinimum    = J
 919             XaMin               = Xa
 920             if selfA._toStore("APosterioriCovariance"):
 921                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
 922     #
 923     # Stockage final supplémentaire de l'optimum en estimation de paramètres
 924     # ----------------------------------------------------------------------
 925     if selfA._parameters["EstimationOf"] == "Parameters":
 926         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 927         selfA.StoredVariables["Analysis"].store( XaMin )
 928         if selfA._toStore("APosterioriCovariance"):
 929             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
 930         if selfA._toStore("BMA"):
 931             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
 932     #
 933     return 0
 934
 935 # ==============================================================================
 936 def enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterFormula"):
 937     """
 938     EnKS
 939     """
 940     #
 941     # Opérateurs
 942     H = HO["Direct"].appliedControledFormTo
 943     #
 944     if selfA._parameters["EstimationOf"] == "State":
 945         M = EM["Direct"].appliedControledFormTo
 946     #
 947     if CM is not None and "Tangent" in CM and U is not None:
 948         Cm = CM["Tangent"].asMatrix(Xb)
 949     else:
 950         Cm = None
 951     #
 952     # Précalcul des inversions de B et R
 953     RIdemi = R.sqrtmI()
 954     #
 955     # Durée d'observation et tailles
 956     LagL = selfA._parameters["SmootherLagL"]
 957     if (not hasattr(Y,"store")) or (not hasattr(Y,"stepnumber")):
 958         raise ValueError("Fixed-lag smoother requires a series of observation")
 959     if Y.stepnumber() < LagL:
 960         raise ValueError("Fixed-lag smoother requires a series of observation greater then the lag L")
 961     duration = Y.stepnumber()
 962     __p = numpy.cumprod(Y.shape())[-1]
 963     __n = Xb.size
 964     __m = selfA._parameters["NumberOfMembers"]
 965     #
 966     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
 967         selfA.StoredVariables["Analysis"].store( Xb )
 968         if selfA._toStore("APosterioriCovariance"):
 969             if hasattr(B,"asfullmatrix"):
 970                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
 971             else:
 972                 selfA.StoredVariables["APosterioriCovariance"].store( B )
 973     #
 974     # Calcul direct initial (on privilégie la mémorisation au recalcul)
 975     __seed = numpy.random.get_state()
 976     selfB = copy.deepcopy(selfA)
 977     selfB._parameters["StoreSupplementaryCalculations"] = ["CurrentEnsembleState"]
 978     if VariantM == "EnKS16-KalmanFilterFormula":
 979         etkf(selfB, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM = "KalmanFilterFormula")
 980     else:
 981         raise ValueError("VariantM has to be chosen in the authorized methods list.")
 982     if LagL > 0:
 983         EL  = selfB.StoredVariables["CurrentEnsembleState"][LagL-1]
 984     else:
 985         EL = EnsembleOfBackgroundPerturbations( Xb, None, __m ) # Cf. etkf
 986     selfA._parameters["SetSeed"] = numpy.random.set_state(__seed)
 987     #
 988     for step in range(LagL,duration-1):
 989         #
 990         sEL = selfB.StoredVariables["CurrentEnsembleState"][step+1-LagL:step+1]
 991         sEL.append(None)
 992         #
 993         if hasattr(Y,"store"):
 994             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
 995         else:
 996             Ynpu = numpy.ravel( Y ).reshape((__p,1))
 997         #
 998         if U is not None:
 999             if hasattr(U,"store") and len(U)>1:
1000                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
1001             elif hasattr(U,"store") and len(U)==1:
1002                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
1003             else:
1004                 Un = numpy.asmatrix(numpy.ravel( U )).T
1005         else:
1006             Un = None
1007         #
1008         #--------------------------
1009         if VariantM == "EnKS16-KalmanFilterFormula":
1010             if selfA._parameters["EstimationOf"] == "State": # Forecast
1011                 EL = M( [(EL[:,i], Un) for i in range(__m)],
1012                     argsAsSerie = True,
1013                     returnSerieAsArrayMatrix = True )
1014                 EL = EnsemblePerturbationWithGivenCovariance( EL, Q )
1015                 EZ = H( [(EL[:,i], Un) for i in range(__m)],
1016                     argsAsSerie = True,
1017                     returnSerieAsArrayMatrix = True )
1018                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
1019                     Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
1020                     EZ = EZ + Cm * Un
1021             elif selfA._parameters["EstimationOf"] == "Parameters":
1022                 # --- > Par principe, M = Id, Q = 0
1023                 EZ = H( [(EL[:,i], Un) for i in range(__m)],
1024                     argsAsSerie = True,
1025                     returnSerieAsArrayMatrix = True )
1026             #
1027             vEm   = EL.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1028             vZm   = EZ.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
1029             #
1030             mS    = RIdemi @ EnsembleOfAnomalies( EZ, vZm, 1./math.sqrt(__m-1) )
1031             mS    = mS.reshape((-1,__m)) # Pour dimension 1
1032             delta = RIdemi @ ( Ynpu - vZm )
1033             mT    = numpy.linalg.inv( numpy.identity(__m) + mS.T @ mS )
1034             vw    = mT @ mS.T @ delta
1035             #
1036             Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
1037             mU    = numpy.identity(__m)
1038             wTU   = (vw.reshape((__m,1)) + math.sqrt(__m-1) * Tdemi @ mU)
1039             #
1040             EX    = EnsembleOfAnomalies( EL, vEm, 1./math.sqrt(__m-1) )
1041             EL    = vEm + EX @ wTU
1042             #
1043             sEL[LagL] = EL
1044             for irl in range(LagL): # Lissage des L précédentes analysis
1045                 vEm = sEL[irl].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1046                 EX = EnsembleOfAnomalies( sEL[irl], vEm, 1./math.sqrt(__m-1) )
1047                 sEL[irl] = vEm + EX @ wTU
1048             #
1049             # Conservation de l'analyse retrospective d'ordre 0 avant rotation
1050             Xa = sEL[0].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1051             if selfA._toStore("APosterioriCovariance"):
1052                 EXn = sEL[0]
1053             #
1054             for irl in range(LagL):
1055                 sEL[irl] = sEL[irl+1]
1056             sEL[LagL] = None
1057         #--------------------------
1058         else:
1059             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1060         #
1061         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1062         # ---> avec analysis
1063         selfA.StoredVariables["Analysis"].store( Xa )
1064         if selfA._toStore("APosterioriCovariance"):
1065             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(EXn) )
1066     #
1067     # Stockage des dernières analyses incomplètement remises à jour
1068     for irl in range(LagL):
1069         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1070         Xa = sEL[irl].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1071         selfA.StoredVariables["Analysis"].store( Xa )
1072     #
1073     return 0
1074
1075 # ==============================================================================
1076 def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"):
1077     """
1078     Ensemble-Transform EnKF
1079     """
1080     if selfA._parameters["EstimationOf"] == "Parameters":
1081         selfA._parameters["StoreInternalVariables"] = True
1082     #
1083     # Opérateurs
1084     H = HO["Direct"].appliedControledFormTo
1085     #
1086     if selfA._parameters["EstimationOf"] == "State":
1087         M = EM["Direct"].appliedControledFormTo
1088     #
1089     if CM is not None and "Tangent" in CM and U is not None:
1090         Cm = CM["Tangent"].asMatrix(Xb)
1091     else:
1092         Cm = None
1093     #
1094     # Durée d'observation et tailles
1095     if hasattr(Y,"stepnumber"):
1096         duration = Y.stepnumber()
1097         __p = numpy.cumprod(Y.shape())[-1]
1098     else:
1099         duration = 2
1100         __p = numpy.array(Y).size
1101     #
1102     # Précalcul des inversions de B et R
1103     if selfA._parameters["StoreInternalVariables"] \
1104         or selfA._toStore("CostFunctionJ") \
1105         or selfA._toStore("CostFunctionJb") \
1106         or selfA._toStore("CostFunctionJo") \
1107         or selfA._toStore("CurrentOptimum") \
1108         or selfA._toStore("APosterioriCovariance"):
1109         BI = B.getI()
1110         RI = R.getI()
1111     elif VariantM != "KalmanFilterFormula":
1112         RI = R.getI()
1113     if VariantM == "KalmanFilterFormula":
1114         RIdemi = R.sqrtmI()
1115     #
1116     __n = Xb.size
1117     __m = selfA._parameters["NumberOfMembers"]
1118     #
1119     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
1120         Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
1121         selfA.StoredVariables["Analysis"].store( Xb )
1122         if selfA._toStore("APosterioriCovariance"):
1123             if hasattr(B,"asfullmatrix"):
1124                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
1125             else:
1126                 selfA.StoredVariables["APosterioriCovariance"].store( B )
1127         selfA._setInternalState("seed", numpy.random.get_state())
1128     elif selfA._parameters["nextStep"]:
1129         Xn = selfA._getInternalState("Xn")
1130     #
1131     previousJMinimum = numpy.finfo(float).max
1132     #
1133     for step in range(duration-1):
1134         numpy.random.set_state(selfA._getInternalState("seed"))
1135         if hasattr(Y,"store"):
1136             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
1137         else:
1138             Ynpu = numpy.ravel( Y ).reshape((__p,1))
1139         #
1140         if U is not None:
1141             if hasattr(U,"store") and len(U)>1:
1142                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
1143             elif hasattr(U,"store") and len(U)==1:
1144                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
1145             else:
1146                 Un = numpy.asmatrix(numpy.ravel( U )).T
1147         else:
1148             Un = None
1149         #
1150         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
1151             Xn = CovarianceInflation( Xn,
1152                 selfA._parameters["InflationType"],
1153                 selfA._parameters["InflationFactor"],
1154                 )
1155         #
1156         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
1157             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
1158                 argsAsSerie = True,
1159                 returnSerieAsArrayMatrix = True )
1160             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
1161             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
1162                 argsAsSerie = True,
1163                 returnSerieAsArrayMatrix = True )
1164             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
1165                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
1166                 Xn_predicted = Xn_predicted + Cm * Un
1167         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
1168             # --- > Par principe, M = Id, Q = 0
1169             Xn_predicted = EMX = Xn
1170             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
1171                 argsAsSerie = True,
1172                 returnSerieAsArrayMatrix = True )
1173         #
1174         # Mean of forecast and observation of forecast
1175         Xfm  = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1176         Hfm  = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
1177         #
1178         # Anomalies
1179         EaX   = EnsembleOfAnomalies( Xn_predicted, Xfm )
1180         EaHX  = EnsembleOfAnomalies( HX_predicted, Hfm)
1181         #
1182         #--------------------------
1183         if VariantM == "KalmanFilterFormula":
1184             mS    = RIdemi * EaHX / math.sqrt(__m-1)
1185             mS    = mS.reshape((-1,__m)) # Pour dimension 1
1186             delta = RIdemi * ( Ynpu - Hfm )
1187             mT    = numpy.linalg.inv( numpy.identity(__m) + mS.T @ mS )
1188             vw    = mT @ mS.T @ delta
1189             #
1190             Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
1191             mU    = numpy.identity(__m)
1192             #
1193             EaX   = EaX / math.sqrt(__m-1)
1194             Xn    = Xfm + EaX @ ( vw.reshape((__m,1)) + math.sqrt(__m-1) * Tdemi @ mU )
1195         #--------------------------
1196         elif VariantM == "Variational":
1197             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1198             def CostFunction(w):
1199                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1200                 _Jo = 0.5 * _A.T @ (RI * _A)
1201                 _Jb = 0.5 * (__m-1) * w.T @ w
1202                 _J  = _Jo + _Jb
1203                 return float(_J)
1204             def GradientOfCostFunction(w):
1205                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1206                 _GardJo = - EaHX.T @ (RI * _A)
1207                 _GradJb = (__m-1) * w.reshape((__m,1))
1208                 _GradJ  = _GardJo + _GradJb
1209                 return numpy.ravel(_GradJ)
1210             vw = scipy.optimize.fmin_cg(
1211                 f           = CostFunction,
1212                 x0          = numpy.zeros(__m),
1213                 fprime      = GradientOfCostFunction,
1214                 args        = (),
1215                 disp        = False,
1216                 )
1217             #
1218             Hto = EaHX.T @ (RI * EaHX).reshape((-1,__m))
1219             Htb = (__m-1) * numpy.identity(__m)
1220             Hta = Hto + Htb
1221             #
1222             Pta = numpy.linalg.inv( Hta )
1223             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1224             #
1225             Xn  = Xfm + EaX @ (vw[:,None] + EWa)
1226         #--------------------------
1227         elif VariantM == "FiniteSize11": # Jauge Boc2011
1228             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1229             def CostFunction(w):
1230                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1231                 _Jo = 0.5 * _A.T @ (RI * _A)
1232                 _Jb = 0.5 * __m * math.log(1 + 1/__m + w.T @ w)
1233                 _J  = _Jo + _Jb
1234                 return float(_J)
1235             def GradientOfCostFunction(w):
1236                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1237                 _GardJo = - EaHX.T @ (RI * _A)
1238                 _GradJb = __m * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
1239                 _GradJ  = _GardJo + _GradJb
1240                 return numpy.ravel(_GradJ)
1241             vw = scipy.optimize.fmin_cg(
1242                 f           = CostFunction,
1243                 x0          = numpy.zeros(__m),
1244                 fprime      = GradientOfCostFunction,
1245                 args        = (),
1246                 disp        = False,
1247                 )
1248             #
1249             Hto = EaHX.T @ (RI * EaHX).reshape((-1,__m))
1250             Htb = __m * \
1251                 ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
1252                 / (1 + 1/__m + vw.T @ vw)**2
1253             Hta = Hto + Htb
1254             #
1255             Pta = numpy.linalg.inv( Hta )
1256             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1257             #
1258             Xn  = Xfm + EaX @ (vw.reshape((__m,1)) + EWa)
1259         #--------------------------
1260         elif VariantM == "FiniteSize15": # Jauge Boc2015
1261             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1262             def CostFunction(w):
1263                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1264                 _Jo = 0.5 * _A.T * RI * _A
1265                 _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w)
1266                 _J  = _Jo + _Jb
1267                 return float(_J)
1268             def GradientOfCostFunction(w):
1269                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1270                 _GardJo = - EaHX.T @ (RI * _A)
1271                 _GradJb = (__m+1) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
1272                 _GradJ  = _GardJo + _GradJb
1273                 return numpy.ravel(_GradJ)
1274             vw = scipy.optimize.fmin_cg(
1275                 f           = CostFunction,
1276                 x0          = numpy.zeros(__m),
1277                 fprime      = GradientOfCostFunction,
1278                 args        = (),
1279                 disp        = False,
1280                 )
1281             #
1282             Hto = EaHX.T @ (RI * EaHX).reshape((-1,__m))
1283             Htb = (__m+1) * \
1284                 ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
1285                 / (1 + 1/__m + vw.T @ vw)**2
1286             Hta = Hto + Htb
1287             #
1288             Pta = numpy.linalg.inv( Hta )
1289             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1290             #
1291             Xn  = Xfm + EaX @ (vw.reshape((__m,1)) + EWa)
1292         #--------------------------
1293         elif VariantM == "FiniteSize16": # Jauge Boc2016
1294             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1295             def CostFunction(w):
1296                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1297                 _Jo = 0.5 * _A.T @ (RI * _A)
1298                 _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w / (__m-1))
1299                 _J  = _Jo + _Jb
1300                 return float(_J)
1301             def GradientOfCostFunction(w):
1302                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1303                 _GardJo = - EaHX.T @ (RI * _A)
1304                 _GradJb = ((__m+1) / (__m-1)) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w / (__m-1))
1305                 _GradJ  = _GardJo + _GradJb
1306                 return numpy.ravel(_GradJ)
1307             vw = scipy.optimize.fmin_cg(
1308                 f           = CostFunction,
1309                 x0          = numpy.zeros(__m),
1310                 fprime      = GradientOfCostFunction,
1311                 args        = (),
1312                 disp        = False,
1313                 )
1314             #
1315             Hto = EaHX.T @ (RI * EaHX).reshape((-1,__m))
1316             Htb = ((__m+1) / (__m-1)) * \
1317                 ( (1 + 1/__m + vw.T @ vw / (__m-1)) * numpy.identity(__m) - 2 * vw @ vw.T / (__m-1) ) \
1318                 / (1 + 1/__m + vw.T @ vw / (__m-1))**2
1319             Hta = Hto + Htb
1320             #
1321             Pta = numpy.linalg.inv( Hta )
1322             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1323             #
1324             Xn  = Xfm + EaX @ (vw[:,None] + EWa)
1325         #--------------------------
1326         else:
1327             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1328         #
1329         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
1330             Xn = CovarianceInflation( Xn,
1331                 selfA._parameters["InflationType"],
1332                 selfA._parameters["InflationFactor"],
1333                 )
1334         #
1335         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1336         #--------------------------
1337         selfA._setInternalState("Xn", Xn)
1338         selfA._setInternalState("seed", numpy.random.get_state())
1339         #--------------------------
1340         #
1341         if selfA._parameters["StoreInternalVariables"] \
1342             or selfA._toStore("CostFunctionJ") \
1343             or selfA._toStore("CostFunctionJb") \
1344             or selfA._toStore("CostFunctionJo") \
1345             or selfA._toStore("APosterioriCovariance") \
1346             or selfA._toStore("InnovationAtCurrentAnalysis") \
1347             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
1348             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1349             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
1350             _Innovation = Ynpu - _HXa
1351         #
1352         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1353         # ---> avec analysis
1354         selfA.StoredVariables["Analysis"].store( Xa )
1355         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1356             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
1357         if selfA._toStore("InnovationAtCurrentAnalysis"):
1358             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1359         # ---> avec current state
1360         if selfA._parameters["StoreInternalVariables"] \
1361             or selfA._toStore("CurrentState"):
1362             selfA.StoredVariables["CurrentState"].store( Xn )
1363         if selfA._toStore("ForecastState"):
1364             selfA.StoredVariables["ForecastState"].store( EMX )
1365         if selfA._toStore("ForecastCovariance"):
1366             selfA.StoredVariables["ForecastCovariance"].store( EnsembleErrorCovariance(EMX) )
1367         if selfA._toStore("BMA"):
1368             selfA.StoredVariables["BMA"].store( EMX - Xa )
1369         if selfA._toStore("InnovationAtCurrentState"):
1370             selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
1371         if selfA._toStore("SimulatedObservationAtCurrentState") \
1372             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1373             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
1374         # ---> autres
1375         if selfA._parameters["StoreInternalVariables"] \
1376             or selfA._toStore("CostFunctionJ") \
1377             or selfA._toStore("CostFunctionJb") \
1378             or selfA._toStore("CostFunctionJo") \
1379             or selfA._toStore("CurrentOptimum") \
1380             or selfA._toStore("APosterioriCovariance"):
1381             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
1382             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
1383             J   = Jb + Jo
1384             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1385             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1386             selfA.StoredVariables["CostFunctionJ" ].store( J )
1387             #
1388             if selfA._toStore("IndexOfOptimum") \
1389                 or selfA._toStore("CurrentOptimum") \
1390                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1391                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1392                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1393                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1394                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1395             if selfA._toStore("IndexOfOptimum"):
1396                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1397             if selfA._toStore("CurrentOptimum"):
1398                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1399             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1400                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1401             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1402                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1403             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1404                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1405             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1406                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1407         if selfA._toStore("APosterioriCovariance"):
1408             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
1409         if selfA._parameters["EstimationOf"] == "Parameters" \
1410             and J < previousJMinimum:
1411             previousJMinimum    = J
1412             XaMin               = Xa
1413             if selfA._toStore("APosterioriCovariance"):
1414                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
1415         # ---> Pour les smoothers
1416         if selfA._toStore("CurrentEnsembleState"):
1417             selfA.StoredVariables["CurrentEnsembleState"].store( Xn )
1418     #
1419     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1420     # ----------------------------------------------------------------------
1421     if selfA._parameters["EstimationOf"] == "Parameters":
1422         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1423         selfA.StoredVariables["Analysis"].store( XaMin )
1424         if selfA._toStore("APosterioriCovariance"):
1425             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1426         if selfA._toStore("BMA"):
1427             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1428     #
1429     return 0
1430
1431 # ==============================================================================
1432 def exkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
1433     """
1434     Extended Kalman Filter
1435     """
1436     if selfA._parameters["EstimationOf"] == "Parameters":
1437         selfA._parameters["StoreInternalVariables"] = True
1438     #
1439     # Opérateurs
1440     H = HO["Direct"].appliedControledFormTo
1441     #
1442     if selfA._parameters["EstimationOf"] == "State":
1443         M = EM["Direct"].appliedControledFormTo
1444     #
1445     if CM is not None and "Tangent" in CM and U is not None:
1446         Cm = CM["Tangent"].asMatrix(Xb)
1447     else:
1448         Cm = None
1449     #
1450     # Durée d'observation et tailles
1451     if hasattr(Y,"stepnumber"):
1452         duration = Y.stepnumber()
1453         __p = numpy.cumprod(Y.shape())[-1]
1454     else:
1455         duration = 2
1456         __p = numpy.array(Y).size
1457     #
1458     # Précalcul des inversions de B et R
1459     if selfA._parameters["StoreInternalVariables"] \
1460         or selfA._toStore("CostFunctionJ") \
1461         or selfA._toStore("CostFunctionJb") \
1462         or selfA._toStore("CostFunctionJo") \
1463         or selfA._toStore("CurrentOptimum") \
1464         or selfA._toStore("APosterioriCovariance"):
1465         BI = B.getI()
1466         RI = R.getI()
1467     #
1468     __n = Xb.size
1469     #
1470     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
1471         Xn = Xb
1472         Pn = B
1473         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1474         selfA.StoredVariables["Analysis"].store( Xb )
1475         if selfA._toStore("APosterioriCovariance"):
1476             if hasattr(B,"asfullmatrix"):
1477                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
1478             else:
1479                 selfA.StoredVariables["APosterioriCovariance"].store( B )
1480         selfA._setInternalState("seed", numpy.random.get_state())
1481     elif selfA._parameters["nextStep"]:
1482         Xn = selfA._getInternalState("Xn")
1483         Pn = selfA._getInternalState("Pn")
1484     #
1485     if selfA._parameters["EstimationOf"] == "Parameters":
1486         XaMin            = Xn
1487         previousJMinimum = numpy.finfo(float).max
1488     #
1489     for step in range(duration-1):
1490         if hasattr(Y,"store"):
1491             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
1492         else:
1493             Ynpu = numpy.ravel( Y ).reshape((__p,1))
1494         #
1495         Ht = HO["Tangent"].asMatrix(ValueForMethodForm = Xn)
1496         Ht = Ht.reshape(Ynpu.size,Xn.size) # ADAO & check shape
1497         Ha = HO["Adjoint"].asMatrix(ValueForMethodForm = Xn)
1498         Ha = Ha.reshape(Xn.size,Ynpu.size) # ADAO & check shape
1499         #
1500         if selfA._parameters["EstimationOf"] == "State":
1501             Mt = EM["Tangent"].asMatrix(ValueForMethodForm = Xn)
1502             Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
1503             Ma = EM["Adjoint"].asMatrix(ValueForMethodForm = Xn)
1504             Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
1505         #
1506         if U is not None:
1507             if hasattr(U,"store") and len(U)>1:
1508                 Un = numpy.ravel( U[step] ).reshape((-1,1))
1509             elif hasattr(U,"store") and len(U)==1:
1510                 Un = numpy.ravel( U[0] ).reshape((-1,1))
1511             else:
1512                 Un = numpy.ravel( U ).reshape((-1,1))
1513         else:
1514             Un = None
1515         #
1516         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
1517             Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).reshape((__n,1))
1518             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
1519                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
1520                 Xn_predicted = Xn_predicted + Cm * Un
1521             Pn_predicted = Q + Mt * Pn * Ma
1522         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
1523             # --- > Par principe, M = Id, Q = 0
1524             Xn_predicted = Xn
1525             Pn_predicted = Pn
1526         #
1527         if selfA._parameters["EstimationOf"] == "State":
1528             HX_predicted = numpy.ravel( H( (Xn_predicted, None) ) ).reshape((__p,1))
1529             _Innovation  = Ynpu - HX_predicted
1530         elif selfA._parameters["EstimationOf"] == "Parameters":
1531             HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1))
1532             _Innovation  = Ynpu - HX_predicted
1533             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
1534                 _Innovation = _Innovation - Cm * Un
1535         #
1536         Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
1537         Xn = Xn_predicted + Kn * _Innovation
1538         Pn = Pn_predicted - Kn * Ht * Pn_predicted
1539         #
1540         Xa = Xn # Pointeurs
1541         #--------------------------
1542         selfA._setInternalState("Xn", Xn)
1543         selfA._setInternalState("Pn", Pn)
1544         #--------------------------
1545         #
1546         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1547         # ---> avec analysis
1548         selfA.StoredVariables["Analysis"].store( Xa )
1549         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1550             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( H((Xa, Un)) )
1551         if selfA._toStore("InnovationAtCurrentAnalysis"):
1552             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1553         # ---> avec current state
1554         if selfA._parameters["StoreInternalVariables"] \
1555             or selfA._toStore("CurrentState"):
1556             selfA.StoredVariables["CurrentState"].store( Xn )
1557         if selfA._toStore("ForecastState"):
1558             selfA.StoredVariables["ForecastState"].store( Xn_predicted )
1559         if selfA._toStore("ForecastCovariance"):
1560             selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
1561         if selfA._toStore("BMA"):
1562             selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
1563         if selfA._toStore("InnovationAtCurrentState"):
1564             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
1565         if selfA._toStore("SimulatedObservationAtCurrentState") \
1566             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1567             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
1568         # ---> autres
1569         if selfA._parameters["StoreInternalVariables"] \
1570             or selfA._toStore("CostFunctionJ") \
1571             or selfA._toStore("CostFunctionJb") \
1572             or selfA._toStore("CostFunctionJo") \
1573             or selfA._toStore("CurrentOptimum") \
1574             or selfA._toStore("APosterioriCovariance"):
1575             Jb  = float( 0.5 * (Xa - Xb).T @ (BI @ (Xa - Xb)) )
1576             Jo  = float( 0.5 * _Innovation.T @ (RI @ _Innovation) )
1577             J   = Jb + Jo
1578             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1579             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1580             selfA.StoredVariables["CostFunctionJ" ].store( J )
1581             #
1582             if selfA._toStore("IndexOfOptimum") \
1583                 or selfA._toStore("CurrentOptimum") \
1584                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1585                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1586                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1587                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1588                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1589             if selfA._toStore("IndexOfOptimum"):
1590                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1591             if selfA._toStore("CurrentOptimum"):
1592                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1593             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1594                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1595             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1596                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1597             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1598                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1599             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1600                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1601         if selfA._toStore("APosterioriCovariance"):
1602             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
1603         if selfA._parameters["EstimationOf"] == "Parameters" \
1604             and J < previousJMinimum:
1605             previousJMinimum    = J
1606             XaMin               = Xa
1607             if selfA._toStore("APosterioriCovariance"):
1608                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
1609     #
1610     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1611     # ----------------------------------------------------------------------
1612     if selfA._parameters["EstimationOf"] == "Parameters":
1613         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1614         selfA.StoredVariables["Analysis"].store( XaMin )
1615         if selfA._toStore("APosterioriCovariance"):
1616             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1617         if selfA._toStore("BMA"):
1618             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1619     #
1620     return 0
1621
1622 # ==============================================================================
1623 def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12",
1624     BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000):
1625     """
1626     Iterative EnKF
1627     """
1628     if selfA._parameters["EstimationOf"] == "Parameters":
1629         selfA._parameters["StoreInternalVariables"] = True
1630     #
1631     # Opérateurs
1632     H = HO["Direct"].appliedControledFormTo
1633     #
1634     if selfA._parameters["EstimationOf"] == "State":
1635         M = EM["Direct"].appliedControledFormTo
1636     #
1637     if CM is not None and "Tangent" in CM and U is not None:
1638         Cm = CM["Tangent"].asMatrix(Xb)
1639     else:
1640         Cm = None
1641     #
1642     # Durée d'observation et tailles
1643     if hasattr(Y,"stepnumber"):
1644         duration = Y.stepnumber()
1645         __p = numpy.cumprod(Y.shape())[-1]
1646     else:
1647         duration = 2
1648         __p = numpy.array(Y).size
1649     #
1650     # Précalcul des inversions de B et R
1651     if selfA._parameters["StoreInternalVariables"] \
1652         or selfA._toStore("CostFunctionJ") \
1653         or selfA._toStore("CostFunctionJb") \
1654         or selfA._toStore("CostFunctionJo") \
1655         or selfA._toStore("CurrentOptimum") \
1656         or selfA._toStore("APosterioriCovariance"):
1657         BI = B.getI()
1658     RI = R.getI()
1659     #
1660     __n = Xb.size
1661     __m = selfA._parameters["NumberOfMembers"]
1662     #
1663     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
1664         if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
1665         else:                         Pn = B
1666         Xn = EnsembleOfBackgroundPerturbations( Xb, Pn, __m )
1667         selfA.StoredVariables["Analysis"].store( Xb )
1668         if selfA._toStore("APosterioriCovariance"):
1669             if hasattr(B,"asfullmatrix"):
1670                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
1671             else:
1672                 selfA.StoredVariables["APosterioriCovariance"].store( B )
1673         selfA._setInternalState("seed", numpy.random.get_state())
1674     elif selfA._parameters["nextStep"]:
1675         Xn = selfA._getInternalState("Xn")
1676     #
1677     previousJMinimum = numpy.finfo(float).max
1678     #
1679     for step in range(duration-1):
1680         numpy.random.set_state(selfA._getInternalState("seed"))
1681         if hasattr(Y,"store"):
1682             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
1683         else:
1684             Ynpu = numpy.ravel( Y ).reshape((__p,1))
1685         #
1686         if U is not None:
1687             if hasattr(U,"store") and len(U)>1:
1688                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
1689             elif hasattr(U,"store") and len(U)==1:
1690                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
1691             else:
1692                 Un = numpy.asmatrix(numpy.ravel( U )).T
1693         else:
1694             Un = None
1695         #
1696         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
1697             Xn = CovarianceInflation( Xn,
1698                 selfA._parameters["InflationType"],
1699                 selfA._parameters["InflationFactor"],
1700                 )
1701         #
1702         #--------------------------
1703         if VariantM == "IEnKF12":
1704             Xfm = numpy.ravel(Xn.mean(axis=1, dtype=mfp).astype('float'))
1705             EaX = EnsembleOfAnomalies( Xn ) / math.sqrt(__m-1)
1706             __j = 0
1707             Deltaw = 1
1708             if not BnotT:
1709                 Ta  = numpy.identity(__m)
1710             vw  = numpy.zeros(__m)
1711             while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
1712                 vx1 = (Xfm + EaX @ vw).reshape((__n,1))
1713                 #
1714                 if BnotT:
1715                     E1 = vx1 + _epsilon * EaX
1716                 else:
1717                     E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
1718                 #
1719                 if selfA._parameters["EstimationOf"] == "State": # Forecast + Q
1720                     E2 = M( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
1721                         argsAsSerie = True,
1722                         returnSerieAsArrayMatrix = True )
1723                 elif selfA._parameters["EstimationOf"] == "Parameters":
1724                     # --- > Par principe, M = Id
1725                     E2 = Xn
1726                 vx2 = E2.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1727                 vy1 = H((vx2, Un)).reshape((__p,1))
1728                 #
1729                 HE2 = H( [(E2[:,i,numpy.newaxis], Un) for i in range(__m)],
1730                     argsAsSerie = True,
1731                     returnSerieAsArrayMatrix = True )
1732                 vy2 = HE2.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
1733                 #
1734                 if BnotT:
1735                     EaY = (HE2 - vy2) / _epsilon
1736                 else:
1737                     EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
1738                 #
1739                 GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy1 )))
1740                 mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY).reshape((-1,__m))
1741                 Deltaw = - numpy.linalg.solve(mH,GradJ)
1742                 #
1743                 vw = vw + Deltaw
1744                 #
1745                 if not BnotT:
1746                     Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1747                 #
1748                 __j = __j + 1
1749             #
1750             A2 = EnsembleOfAnomalies( E2 )
1751             #
1752             if BnotT:
1753                 Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1754                 A2 = math.sqrt(__m-1) * A2 @ Ta / _epsilon
1755             #
1756             Xn = vx2 + A2
1757         #--------------------------
1758         else:
1759             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1760         #
1761         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
1762             Xn = CovarianceInflation( Xn,
1763                 selfA._parameters["InflationType"],
1764                 selfA._parameters["InflationFactor"],
1765                 )
1766         #
1767         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1768         #--------------------------
1769         selfA._setInternalState("Xn", Xn)
1770         selfA._setInternalState("seed", numpy.random.get_state())
1771         #--------------------------
1772         #
1773         if selfA._parameters["StoreInternalVariables"] \
1774             or selfA._toStore("CostFunctionJ") \
1775             or selfA._toStore("CostFunctionJb") \
1776             or selfA._toStore("CostFunctionJo") \
1777             or selfA._toStore("APosterioriCovariance") \
1778             or selfA._toStore("InnovationAtCurrentAnalysis") \
1779             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
1780             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1781             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
1782             _Innovation = Ynpu - _HXa
1783         #
1784         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1785         # ---> avec analysis
1786         selfA.StoredVariables["Analysis"].store( Xa )
1787         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1788             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
1789         if selfA._toStore("InnovationAtCurrentAnalysis"):
1790             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1791         # ---> avec current state
1792         if selfA._parameters["StoreInternalVariables"] \
1793             or selfA._toStore("CurrentState"):
1794             selfA.StoredVariables["CurrentState"].store( Xn )
1795         if selfA._toStore("ForecastState"):
1796             selfA.StoredVariables["ForecastState"].store( E2 )
1797         if selfA._toStore("ForecastCovariance"):
1798             selfA.StoredVariables["ForecastCovariance"].store( EnsembleErrorCovariance(E2) )
1799         if selfA._toStore("BMA"):
1800             selfA.StoredVariables["BMA"].store( E2 - Xa )
1801         if selfA._toStore("InnovationAtCurrentState"):
1802             selfA.StoredVariables["InnovationAtCurrentState"].store( - HE2 + Ynpu )
1803         if selfA._toStore("SimulatedObservationAtCurrentState") \
1804             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1805             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HE2 )
1806         # ---> autres
1807         if selfA._parameters["StoreInternalVariables"] \
1808             or selfA._toStore("CostFunctionJ") \
1809             or selfA._toStore("CostFunctionJb") \
1810             or selfA._toStore("CostFunctionJo") \
1811             or selfA._toStore("CurrentOptimum") \
1812             or selfA._toStore("APosterioriCovariance"):
1813             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
1814             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
1815             J   = Jb + Jo
1816             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1817             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1818             selfA.StoredVariables["CostFunctionJ" ].store( J )
1819             #
1820             if selfA._toStore("IndexOfOptimum") \
1821                 or selfA._toStore("CurrentOptimum") \
1822                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1823                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1824                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1825                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1826                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1827             if selfA._toStore("IndexOfOptimum"):
1828                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1829             if selfA._toStore("CurrentOptimum"):
1830                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1831             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1832                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1833             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1834                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1835             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1836                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1837             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1838                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1839         if selfA._toStore("APosterioriCovariance"):
1840             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
1841         if selfA._parameters["EstimationOf"] == "Parameters" \
1842             and J < previousJMinimum:
1843             previousJMinimum    = J
1844             XaMin               = Xa
1845             if selfA._toStore("APosterioriCovariance"):
1846                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
1847     #
1848     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1849     # ----------------------------------------------------------------------
1850     if selfA._parameters["EstimationOf"] == "Parameters":
1851         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1852         selfA.StoredVariables["Analysis"].store( XaMin )
1853         if selfA._toStore("APosterioriCovariance"):
1854             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1855         if selfA._toStore("BMA"):
1856             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1857     #
1858     return 0
1859
1860 # ==============================================================================
1861 def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
1862     """
1863     3DVAR incrémental
1864     """
1865     #
1866     # Initialisations
1867     # ---------------
1868     #
1869     # Opérateur non-linéaire pour la boucle externe
1870     Hm = HO["Direct"].appliedTo
1871     #
1872     # Précalcul des inversions de B et R
1873     BI = B.getI()
1874     RI = R.getI()
1875     #
1876     # Point de démarrage de l'optimisation
1877     Xini = selfA._parameters["InitializationPoint"]
1878     #
1879     HXb = numpy.asmatrix(numpy.ravel( Hm( Xb ) )).T
1880     Innovation = Y - HXb
1881     #
1882     # Outer Loop
1883     # ----------
1884     iOuter = 0
1885     J      = 1./mpr
1886     DeltaJ = 1./mpr
1887     Xr     = Xini.reshape((-1,1))
1888     while abs(DeltaJ) >= selfA._parameters["CostDecrementTolerance"] and iOuter <= selfA._parameters["MaximumNumberOfSteps"]:
1889         #
1890         # Inner Loop
1891         # ----------
1892         Ht = HO["Tangent"].asMatrix(Xr)
1893         Ht = Ht.reshape(Y.size,Xr.size) # ADAO & check shape
1894         #
1895         # Définition de la fonction-coût
1896         # ------------------------------
1897         def CostFunction(dx):
1898             _dX  = numpy.asmatrix(numpy.ravel( dx )).T
1899             if selfA._parameters["StoreInternalVariables"] or \
1900                 selfA._toStore("CurrentState") or \
1901                 selfA._toStore("CurrentOptimum"):
1902                 selfA.StoredVariables["CurrentState"].store( Xb + _dX )
1903             _HdX = Ht * _dX
1904             _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T
1905             _dInnovation = Innovation - _HdX
1906             if selfA._toStore("SimulatedObservationAtCurrentState") or \
1907                 selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1908                 selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HXb + _HdX )
1909             if selfA._toStore("InnovationAtCurrentState"):
1910                 selfA.StoredVariables["InnovationAtCurrentState"].store( _dInnovation )
1911             #
1912             Jb  = float( 0.5 * _dX.T * BI * _dX )
1913             Jo  = float( 0.5 * _dInnovation.T * RI * _dInnovation )
1914             J   = Jb + Jo
1915             #
1916             selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
1917             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1918             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1919             selfA.StoredVariables["CostFunctionJ" ].store( J )
1920             if selfA._toStore("IndexOfOptimum") or \
1921                 selfA._toStore("CurrentOptimum") or \
1922                 selfA._toStore("CostFunctionJAtCurrentOptimum") or \
1923                 selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
1924                 selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
1925                 selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1926                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1927             if selfA._toStore("IndexOfOptimum"):
1928                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1929             if selfA._toStore("CurrentOptimum"):
1930                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
1931             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1932                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
1933             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1934                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1935             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1936                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1937             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1938                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1939             return J
1940         #
1941         def GradientOfCostFunction(dx):
1942             _dX          = numpy.asmatrix(numpy.ravel( dx )).T
1943             _HdX         = Ht * _dX
1944             _HdX         = numpy.asmatrix(numpy.ravel( _HdX )).T
1945             _dInnovation = Innovation - _HdX
1946             GradJb       = BI * _dX
1947             GradJo       = - Ht.T @ (RI * _dInnovation)
1948             GradJ        = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
1949             return GradJ
1950         #
1951         # Minimisation de la fonctionnelle
1952         # --------------------------------
1953         nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
1954         #
1955         if selfA._parameters["Minimizer"] == "LBFGSB":
1956             # Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
1957             if "0.19" <= scipy.version.version <= "1.1.0":
1958                 import lbfgsbhlt as optimiseur
1959             else:
1960                 import scipy.optimize as optimiseur
1961             Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
1962                 func        = CostFunction,
1963                 x0          = numpy.zeros(Xini.size),
1964                 fprime      = GradientOfCostFunction,
1965                 args        = (),
1966                 bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
1967                 maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
1968                 factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
1969                 pgtol       = selfA._parameters["ProjectedGradientTolerance"],
1970                 iprint      = selfA._parameters["optiprint"],
1971                 )
1972             nfeval = Informations['funcalls']
1973             rc     = Informations['warnflag']
1974         elif selfA._parameters["Minimizer"] == "TNC":
1975             Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
1976                 func        = CostFunction,
1977                 x0          = numpy.zeros(Xini.size),
1978                 fprime      = GradientOfCostFunction,
1979                 args        = (),
1980                 bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
1981                 maxfun      = selfA._parameters["MaximumNumberOfSteps"],
1982                 pgtol       = selfA._parameters["ProjectedGradientTolerance"],
1983                 ftol        = selfA._parameters["CostDecrementTolerance"],
1984                 messages    = selfA._parameters["optmessages"],
1985                 )
1986         elif selfA._parameters["Minimizer"] == "CG":
1987             Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
1988                 f           = CostFunction,
1989                 x0          = numpy.zeros(Xini.size),
1990                 fprime      = GradientOfCostFunction,
1991                 args        = (),
1992                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
1993                 gtol        = selfA._parameters["GradientNormTolerance"],
1994                 disp        = selfA._parameters["optdisp"],
1995                 full_output = True,
1996                 )
1997         elif selfA._parameters["Minimizer"] == "NCG":
1998             Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
1999                 f           = CostFunction,
2000                 x0          = numpy.zeros(Xini.size),
2001                 fprime      = GradientOfCostFunction,
2002                 args        = (),
2003                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2004                 avextol     = selfA._parameters["CostDecrementTolerance"],
2005                 disp        = selfA._parameters["optdisp"],
2006                 full_output = True,
2007                 )
2008         elif selfA._parameters["Minimizer"] == "BFGS":
2009             Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
2010                 f           = CostFunction,
2011                 x0          = numpy.zeros(Xini.size),
2012                 fprime      = GradientOfCostFunction,
2013                 args        = (),
2014                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2015                 gtol        = selfA._parameters["GradientNormTolerance"],
2016                 disp        = selfA._parameters["optdisp"],
2017                 full_output = True,
2018                 )
2019         else:
2020             raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
2021         #
2022         IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2023         MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
2024         #
2025         if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
2026             Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
2027             Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
2028         else:
2029             Minimum = Xb + numpy.asmatrix(numpy.ravel( Minimum )).T
2030         #
2031         Xr     = Minimum
2032         DeltaJ = selfA.StoredVariables["CostFunctionJ" ][-1] - J
2033         iOuter = selfA.StoredVariables["CurrentIterationNumber"][-1]
2034     #
2035     # Obtention de l'analyse
2036     # ----------------------
2037     Xa = Xr
2038     #
2039     selfA.StoredVariables["Analysis"].store( Xa )
2040     #
2041     if selfA._toStore("OMA") or \
2042         selfA._toStore("SigmaObs2") or \
2043         selfA._toStore("SimulationQuantiles") or \
2044         selfA._toStore("SimulatedObservationAtOptimum"):
2045         if selfA._toStore("SimulatedObservationAtCurrentState"):
2046             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
2047         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2048             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
2049         else:
2050             HXa = Hm( Xa )
2051     #
2052     # Calcul de la covariance d'analyse
2053     # ---------------------------------
2054     if selfA._toStore("APosterioriCovariance") or \
2055         selfA._toStore("SimulationQuantiles") or \
2056         selfA._toStore("JacobianMatrixAtOptimum") or \
2057         selfA._toStore("KalmanGainAtOptimum"):
2058         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
2059         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
2060     if selfA._toStore("APosterioriCovariance") or \
2061         selfA._toStore("SimulationQuantiles") or \
2062         selfA._toStore("KalmanGainAtOptimum"):
2063         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
2064         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
2065     if selfA._toStore("APosterioriCovariance") or \
2066         selfA._toStore("SimulationQuantiles"):
2067         HessienneI = []
2068         nb = Xa.size
2069         for i in range(nb):
2070             _ee    = numpy.matrix(numpy.zeros(nb)).T
2071             _ee[i] = 1.
2072             _HtEE  = numpy.dot(HtM,_ee)
2073             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
2074             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
2075         HessienneI = numpy.matrix( HessienneI )
2076         A = HessienneI.I
2077         if min(A.shape) != max(A.shape):
2078             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)))
2079         if (numpy.diag(A) < 0).any():
2080             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,))
2081         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
2082             try:
2083                 L = numpy.linalg.cholesky( A )
2084             except:
2085                 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,))
2086     if selfA._toStore("APosterioriCovariance"):
2087         selfA.StoredVariables["APosterioriCovariance"].store( A )
2088     if selfA._toStore("JacobianMatrixAtOptimum"):
2089         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
2090     if selfA._toStore("KalmanGainAtOptimum"):
2091         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
2092         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
2093         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
2094     #
2095     # Calculs et/ou stockages supplémentaires
2096     # ---------------------------------------
2097     if selfA._toStore("Innovation") or \
2098         selfA._toStore("SigmaObs2") or \
2099         selfA._toStore("MahalanobisConsistency") or \
2100         selfA._toStore("OMB"):
2101         d  = Y - HXb
2102     if selfA._toStore("Innovation"):
2103         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
2104     if selfA._toStore("BMA"):
2105         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
2106     if selfA._toStore("OMA"):
2107         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
2108     if selfA._toStore("OMB"):
2109         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
2110     if selfA._toStore("SigmaObs2"):
2111         TraceR = R.trace(Y.size)
2112         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
2113     if selfA._toStore("MahalanobisConsistency"):
2114         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
2115     if selfA._toStore("SimulationQuantiles"):
2116         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
2117     if selfA._toStore("SimulatedObservationAtBackground"):
2118         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
2119     if selfA._toStore("SimulatedObservationAtOptimum"):
2120         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
2121     #
2122     return 0
2123
2124 # ==============================================================================
2125 def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13",
2126     BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000):
2127     """
2128     Maximum Likelihood Ensemble Filter
2129     """
2130     if selfA._parameters["EstimationOf"] == "Parameters":
2131         selfA._parameters["StoreInternalVariables"] = True
2132     #
2133     # Opérateurs
2134     H = HO["Direct"].appliedControledFormTo
2135     #
2136     if selfA._parameters["EstimationOf"] == "State":
2137         M = EM["Direct"].appliedControledFormTo
2138     #
2139     if CM is not None and "Tangent" in CM and U is not None:
2140         Cm = CM["Tangent"].asMatrix(Xb)
2141     else:
2142         Cm = None
2143     #
2144     # Durée d'observation et tailles
2145     if hasattr(Y,"stepnumber"):
2146         duration = Y.stepnumber()
2147         __p = numpy.cumprod(Y.shape())[-1]
2148     else:
2149         duration = 2
2150         __p = numpy.array(Y).size
2151     #
2152     # Précalcul des inversions de B et R
2153     if selfA._parameters["StoreInternalVariables"] \
2154         or selfA._toStore("CostFunctionJ") \
2155         or selfA._toStore("CostFunctionJb") \
2156         or selfA._toStore("CostFunctionJo") \
2157         or selfA._toStore("CurrentOptimum") \
2158         or selfA._toStore("APosterioriCovariance"):
2159         BI = B.getI()
2160     RI = R.getI()
2161     #
2162     __n = Xb.size
2163     __m = selfA._parameters["NumberOfMembers"]
2164     #
2165     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
2166         Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
2167         selfA.StoredVariables["Analysis"].store( Xb )
2168         if selfA._toStore("APosterioriCovariance"):
2169             if hasattr(B,"asfullmatrix"):
2170                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
2171             else:
2172                 selfA.StoredVariables["APosterioriCovariance"].store( B )
2173         selfA._setInternalState("seed", numpy.random.get_state())
2174     elif selfA._parameters["nextStep"]:
2175         Xn = selfA._getInternalState("Xn")
2176     #
2177     previousJMinimum = numpy.finfo(float).max
2178     #
2179     for step in range(duration-1):
2180         numpy.random.set_state(selfA._getInternalState("seed"))
2181         if hasattr(Y,"store"):
2182             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
2183         else:
2184             Ynpu = numpy.ravel( Y ).reshape((__p,1))
2185         #
2186         if U is not None:
2187             if hasattr(U,"store") and len(U)>1:
2188                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
2189             elif hasattr(U,"store") and len(U)==1:
2190                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
2191             else:
2192                 Un = numpy.asmatrix(numpy.ravel( U )).T
2193         else:
2194             Un = None
2195         #
2196         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
2197             Xn = CovarianceInflation( Xn,
2198                 selfA._parameters["InflationType"],
2199                 selfA._parameters["InflationFactor"],
2200                 )
2201         #
2202         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
2203             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
2204                 argsAsSerie = True,
2205                 returnSerieAsArrayMatrix = True )
2206             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
2207             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
2208                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
2209                 Xn_predicted = Xn_predicted + Cm * Un
2210         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
2211             # --- > Par principe, M = Id, Q = 0
2212             Xn_predicted = EMX = Xn
2213         #
2214         #--------------------------
2215         if VariantM == "MLEF13":
2216             Xfm = numpy.ravel(Xn_predicted.mean(axis=1, dtype=mfp).astype('float'))
2217             EaX = EnsembleOfAnomalies( Xn_predicted, Xfm, 1./math.sqrt(__m-1) )
2218             Ua  = numpy.identity(__m)
2219             __j = 0
2220             Deltaw = 1
2221             if not BnotT:
2222                 Ta  = numpy.identity(__m)
2223             vw  = numpy.zeros(__m)
2224             while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
2225                 vx1 = (Xfm + EaX @ vw).reshape((__n,1))
2226                 #
2227                 if BnotT:
2228                     E1 = vx1 + _epsilon * EaX
2229                 else:
2230                     E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
2231                 #
2232                 HE2 = H( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
2233                     argsAsSerie = True,
2234                     returnSerieAsArrayMatrix = True )
2235                 vy2 = HE2.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
2236                 #
2237                 if BnotT:
2238                     EaY = (HE2 - vy2) / _epsilon
2239                 else:
2240                     EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
2241                 #
2242                 GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy2 )))
2243                 mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY).reshape((-1,__m))
2244                 Deltaw = - numpy.linalg.solve(mH,GradJ)
2245                 #
2246                 vw = vw + Deltaw
2247                 #
2248                 if not BnotT:
2249                     Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
2250                 #
2251                 __j = __j + 1
2252             #
2253             if BnotT:
2254                 Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
2255             #
2256             Xn = vx1 + math.sqrt(__m-1) * EaX @ Ta @ Ua
2257         #--------------------------
2258         else:
2259             raise ValueError("VariantM has to be chosen in the authorized methods list.")
2260         #
2261         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
2262             Xn = CovarianceInflation( Xn,
2263                 selfA._parameters["InflationType"],
2264                 selfA._parameters["InflationFactor"],
2265                 )
2266         #
2267         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
2268         #--------------------------
2269         selfA._setInternalState("Xn", Xn)
2270         selfA._setInternalState("seed", numpy.random.get_state())
2271         #--------------------------
2272         #
2273         if selfA._parameters["StoreInternalVariables"] \
2274             or selfA._toStore("CostFunctionJ") \
2275             or selfA._toStore("CostFunctionJb") \
2276             or selfA._toStore("CostFunctionJo") \
2277             or selfA._toStore("APosterioriCovariance") \
2278             or selfA._toStore("InnovationAtCurrentAnalysis") \
2279             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
2280             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2281             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
2282             _Innovation = Ynpu - _HXa
2283         #
2284         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2285         # ---> avec analysis
2286         selfA.StoredVariables["Analysis"].store( Xa )
2287         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
2288             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
2289         if selfA._toStore("InnovationAtCurrentAnalysis"):
2290             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
2291         # ---> avec current state
2292         if selfA._parameters["StoreInternalVariables"] \
2293             or selfA._toStore("CurrentState"):
2294             selfA.StoredVariables["CurrentState"].store( Xn )
2295         if selfA._toStore("ForecastState"):
2296             selfA.StoredVariables["ForecastState"].store( EMX )
2297         if selfA._toStore("ForecastCovariance"):
2298             selfA.StoredVariables["ForecastCovariance"].store( EnsembleErrorCovariance(EMX) )
2299         if selfA._toStore("BMA"):
2300             selfA.StoredVariables["BMA"].store( EMX - Xa )
2301         if selfA._toStore("InnovationAtCurrentState"):
2302             selfA.StoredVariables["InnovationAtCurrentState"].store( - HE2 + Ynpu )
2303         if selfA._toStore("SimulatedObservationAtCurrentState") \
2304             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2305             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HE2 )
2306         # ---> autres
2307         if selfA._parameters["StoreInternalVariables"] \
2308             or selfA._toStore("CostFunctionJ") \
2309             or selfA._toStore("CostFunctionJb") \
2310             or selfA._toStore("CostFunctionJo") \
2311             or selfA._toStore("CurrentOptimum") \
2312             or selfA._toStore("APosterioriCovariance"):
2313             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
2314             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
2315             J   = Jb + Jo
2316             selfA.StoredVariables["CostFunctionJb"].store( Jb )
2317             selfA.StoredVariables["CostFunctionJo"].store( Jo )
2318             selfA.StoredVariables["CostFunctionJ" ].store( J )
2319             #
2320             if selfA._toStore("IndexOfOptimum") \
2321                 or selfA._toStore("CurrentOptimum") \
2322                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
2323                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
2324                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
2325                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2326                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2327             if selfA._toStore("IndexOfOptimum"):
2328                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2329             if selfA._toStore("CurrentOptimum"):
2330                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
2331             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2332                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
2333             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2334                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2335             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2336                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2337             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2338                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2339         if selfA._toStore("APosterioriCovariance"):
2340             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
2341         if selfA._parameters["EstimationOf"] == "Parameters" \
2342             and J < previousJMinimum:
2343             previousJMinimum    = J
2344             XaMin               = Xa
2345             if selfA._toStore("APosterioriCovariance"):
2346                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
2347     #
2348     # Stockage final supplémentaire de l'optimum en estimation de paramètres
2349     # ----------------------------------------------------------------------
2350     if selfA._parameters["EstimationOf"] == "Parameters":
2351         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2352         selfA.StoredVariables["Analysis"].store( XaMin )
2353         if selfA._toStore("APosterioriCovariance"):
2354             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
2355         if selfA._toStore("BMA"):
2356             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
2357     #
2358     return 0
2359
2360 # ==============================================================================
2361 def mmqr(
2362         func     = None,
2363         x0       = None,
2364         fprime   = None,
2365         bounds   = None,
2366         quantile = 0.5,
2367         maxfun   = 15000,
2368         toler    = 1.e-06,
2369         y        = None,
2370         ):
2371     """
2372     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
2373     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
2374     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
2375     """
2376     #
2377     # Recuperation des donnees et informations initiales
2378     # --------------------------------------------------
2379     variables = numpy.ravel( x0 )
2380     mesures   = numpy.ravel( y )
2381     increment = sys.float_info[0]
2382     p         = variables.size
2383     n         = mesures.size
2384     quantile  = float(quantile)
2385     #
2386     # Calcul des parametres du MM
2387     # ---------------------------
2388     tn      = float(toler) / n
2389     e0      = -tn / math.log(tn)
2390     epsilon = (e0-tn)/(1+math.log(e0))
2391     #
2392     # Calculs d'initialisation
2393     # ------------------------
2394     residus  = mesures - numpy.ravel( func( variables ) )
2395     poids    = 1./(epsilon+numpy.abs(residus))
2396     veps     = 1. - 2. * quantile - residus * poids
2397     lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
2398     iteration = 0
2399     #
2400     # Recherche iterative
2401     # -------------------
2402     while (increment > toler) and (iteration < maxfun) :
2403         iteration += 1
2404         #
2405         Derivees  = numpy.array(fprime(variables))
2406         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
2407         DeriveesT = Derivees.transpose()
2408         M         =   numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
2409         SM        =   numpy.transpose(numpy.dot( DeriveesT , veps ))
2410         step      = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
2411         #
2412         variables = variables + step
2413         if bounds is not None:
2414             # Attention : boucle infinie à éviter si un intervalle est trop petit
2415             while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
2416                 step      = step/2.
2417                 variables = variables - step
2418         residus   = mesures - numpy.ravel( func(variables) )
2419         surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
2420         #
2421         while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
2422             step      = step/2.
2423             variables = variables - step
2424             residus   = mesures - numpy.ravel( func(variables) )
2425             surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
2426         #
2427         increment     = lastsurrogate-surrogate
2428         poids         = 1./(epsilon+numpy.abs(residus))
2429         veps          = 1. - 2. * quantile - residus * poids
2430         lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
2431     #
2432     # Mesure d'écart
2433     # --------------
2434     Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
2435     #
2436     return variables, Ecart, [n,p,iteration,increment,0]
2437
2438 # ==============================================================================
2439 def multi3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle):
2440     """
2441     3DVAR multi-pas et multi-méthodes
2442     """
2443     #
2444     # Initialisation
2445     if selfA._parameters["EstimationOf"] == "State":
2446         M = EM["Direct"].appliedTo
2447         #
2448         if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
2449             Xn = numpy.ravel(Xb).reshape((-1,1))
2450             selfA.StoredVariables["Analysis"].store( Xn )
2451             if selfA._toStore("APosterioriCovariance"):
2452                 if hasattr(B,"asfullmatrix"):
2453                     selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
2454                 else:
2455                     selfA.StoredVariables["APosterioriCovariance"].store( B )
2456             if selfA._toStore("ForecastState"):
2457                 selfA.StoredVariables["ForecastState"].store( Xn )
2458         elif selfA._parameters["nextStep"]:
2459             Xn = selfA._getInternalState("Xn")
2460     else:
2461         Xn = numpy.ravel(Xb).reshape((-1,1))
2462     #
2463     if hasattr(Y,"stepnumber"):
2464         duration = Y.stepnumber()
2465     else:
2466         duration = 2
2467     #
2468     # Multi-pas
2469     for step in range(duration-1):
2470         if hasattr(Y,"store"):
2471             Ynpu = numpy.ravel( Y[step+1] ).reshape((-1,1))
2472         else:
2473             Ynpu = numpy.ravel( Y ).reshape((-1,1))
2474         #
2475         if selfA._parameters["EstimationOf"] == "State": # Forecast
2476             Xn_predicted = M( Xn )
2477             if selfA._toStore("ForecastState"):
2478                 selfA.StoredVariables["ForecastState"].store( Xn_predicted )
2479         elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
2480             # --- > Par principe, M = Id, Q = 0
2481             Xn_predicted = Xn
2482         Xn_predicted = numpy.ravel(Xn_predicted).reshape((-1,1))
2483         #
2484         oneCycle(selfA, Xn_predicted, Ynpu, U, HO, None, None, R, B, None)
2485         #
2486         Xn = selfA.StoredVariables["Analysis"][-1]
2487         #--------------------------
2488         selfA._setInternalState("Xn", Xn)
2489     #
2490     return 0
2491
2492 # ==============================================================================
2493 def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
2494     """
2495     3DVAR PSAS
2496     """
2497     #
2498     # Initialisations
2499     # ---------------
2500     #
2501     # Opérateurs
2502     Hm = HO["Direct"].appliedTo
2503     #
2504     # Utilisation éventuelle d'un vecteur H(Xb) précalculé
2505     if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
2506         HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
2507     else:
2508         HXb = Hm( Xb )
2509     HXb = numpy.asmatrix(numpy.ravel( HXb )).T
2510     if Y.size != HXb.size:
2511         raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
2512     if max(Y.shape) != max(HXb.shape):
2513         raise ValueError("The shapes %s of observations Y and %s of observed calculation H(X) are different, they have to be identical."%(Y.shape,HXb.shape))
2514     #
2515     if selfA._toStore("JacobianMatrixAtBackground"):
2516         HtMb = HO["Tangent"].asMatrix(ValueForMethodForm = Xb)
2517         HtMb = HtMb.reshape(Y.size,Xb.size) # ADAO & check shape
2518         selfA.StoredVariables["JacobianMatrixAtBackground"].store( HtMb )
2519     #
2520     Ht = HO["Tangent"].asMatrix(Xb)
2521     BHT = B * Ht.T
2522     HBHTpR = R + Ht * BHT
2523     Innovation = Y - HXb
2524     #
2525     # Point de démarrage de l'optimisation
2526     Xini = numpy.zeros(Xb.shape)
2527     #
2528     # Définition de la fonction-coût
2529     # ------------------------------
2530     def CostFunction(w):
2531         _W = numpy.asmatrix(numpy.ravel( w )).T
2532         if selfA._parameters["StoreInternalVariables"] or \
2533             selfA._toStore("CurrentState") or \
2534             selfA._toStore("CurrentOptimum"):
2535             selfA.StoredVariables["CurrentState"].store( Xb + BHT * _W )
2536         if selfA._toStore("SimulatedObservationAtCurrentState") or \
2537             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2538             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( Hm( Xb + BHT * _W ) )
2539         if selfA._toStore("InnovationAtCurrentState"):
2540             selfA.StoredVariables["InnovationAtCurrentState"].store( Innovation )
2541         #
2542         Jb  = float( 0.5 * _W.T * HBHTpR * _W )
2543         Jo  = float( - _W.T * Innovation )
2544         J   = Jb + Jo
2545         #
2546         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
2547         selfA.StoredVariables["CostFunctionJb"].store( Jb )
2548         selfA.StoredVariables["CostFunctionJo"].store( Jo )
2549         selfA.StoredVariables["CostFunctionJ" ].store( J )
2550         if selfA._toStore("IndexOfOptimum") or \
2551             selfA._toStore("CurrentOptimum") or \
2552             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
2553             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
2554             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
2555             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2556             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2557         if selfA._toStore("IndexOfOptimum"):
2558             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2559         if selfA._toStore("CurrentOptimum"):
2560             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
2561         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2562             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
2563         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2564             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2565         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2566             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2567         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2568             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2569         return J
2570     #
2571     def GradientOfCostFunction(w):
2572         _W = numpy.asmatrix(numpy.ravel( w )).T
2573         GradJb  = HBHTpR * _W
2574         GradJo  = - Innovation
2575         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
2576         return GradJ
2577     #
2578     # Minimisation de la fonctionnelle
2579     # --------------------------------
2580     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
2581     #
2582     if selfA._parameters["Minimizer"] == "LBFGSB":
2583         if "0.19" <= scipy.version.version <= "1.1.0":
2584             import lbfgsbhlt as optimiseur
2585         else:
2586             import scipy.optimize as optimiseur
2587         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
2588             func        = CostFunction,
2589             x0          = Xini,
2590             fprime      = GradientOfCostFunction,
2591             args        = (),
2592             bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
2593             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
2594             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
2595             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2596             iprint      = selfA._parameters["optiprint"],
2597             )
2598         nfeval = Informations['funcalls']
2599         rc     = Informations['warnflag']
2600     elif selfA._parameters["Minimizer"] == "TNC":
2601         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
2602             func        = CostFunction,
2603             x0          = Xini,
2604             fprime      = GradientOfCostFunction,
2605             args        = (),
2606             bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
2607             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
2608             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2609             ftol        = selfA._parameters["CostDecrementTolerance"],
2610             messages    = selfA._parameters["optmessages"],
2611             )
2612     elif selfA._parameters["Minimizer"] == "CG":
2613         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
2614             f           = CostFunction,
2615             x0          = Xini,
2616             fprime      = GradientOfCostFunction,
2617             args        = (),
2618             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2619             gtol        = selfA._parameters["GradientNormTolerance"],
2620             disp        = selfA._parameters["optdisp"],
2621             full_output = True,
2622             )
2623     elif selfA._parameters["Minimizer"] == "NCG":
2624         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
2625             f           = CostFunction,
2626             x0          = Xini,
2627             fprime      = GradientOfCostFunction,
2628             args        = (),
2629             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2630             avextol     = selfA._parameters["CostDecrementTolerance"],
2631             disp        = selfA._parameters["optdisp"],
2632             full_output = True,
2633             )
2634     elif selfA._parameters["Minimizer"] == "BFGS":
2635         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
2636             f           = CostFunction,
2637             x0          = Xini,
2638             fprime      = GradientOfCostFunction,
2639             args        = (),
2640             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2641             gtol        = selfA._parameters["GradientNormTolerance"],
2642             disp        = selfA._parameters["optdisp"],
2643             full_output = True,
2644             )
2645     else:
2646         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
2647     #
2648     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2649     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
2650     #
2651     # Correction pour pallier a un bug de TNC sur le retour du Minimum
2652     # ----------------------------------------------------------------
2653     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
2654         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
2655         Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
2656     else:
2657         Minimum = Xb + BHT * numpy.asmatrix(numpy.ravel( Minimum )).T
2658     #
2659     # Obtention de l'analyse
2660     # ----------------------
2661     Xa = Minimum
2662     #
2663     selfA.StoredVariables["Analysis"].store( Xa )
2664     #
2665     if selfA._toStore("OMA") or \
2666         selfA._toStore("SigmaObs2") or \
2667         selfA._toStore("SimulationQuantiles") or \
2668         selfA._toStore("SimulatedObservationAtOptimum"):
2669         if selfA._toStore("SimulatedObservationAtCurrentState"):
2670             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
2671         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2672             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
2673         else:
2674             HXa = Hm( Xa )
2675     #
2676     # Calcul de la covariance d'analyse
2677     # ---------------------------------
2678     if selfA._toStore("APosterioriCovariance") or \
2679         selfA._toStore("SimulationQuantiles") or \
2680         selfA._toStore("JacobianMatrixAtOptimum") or \
2681         selfA._toStore("KalmanGainAtOptimum"):
2682         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
2683         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
2684     if selfA._toStore("APosterioriCovariance") or \
2685         selfA._toStore("SimulationQuantiles") or \
2686         selfA._toStore("KalmanGainAtOptimum"):
2687         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
2688         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
2689     if selfA._toStore("APosterioriCovariance") or \
2690         selfA._toStore("SimulationQuantiles"):
2691         BI = B.getI()
2692         RI = R.getI()
2693         HessienneI = []
2694         nb = Xa.size
2695         for i in range(nb):
2696             _ee    = numpy.matrix(numpy.zeros(nb)).T
2697             _ee[i] = 1.
2698             _HtEE  = numpy.dot(HtM,_ee)
2699             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
2700             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
2701         HessienneI = numpy.matrix( HessienneI )
2702         A = HessienneI.I
2703         if min(A.shape) != max(A.shape):
2704             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)))
2705         if (numpy.diag(A) < 0).any():
2706             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,))
2707         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
2708             try:
2709                 L = numpy.linalg.cholesky( A )
2710             except:
2711                 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,))
2712     if selfA._toStore("APosterioriCovariance"):
2713         selfA.StoredVariables["APosterioriCovariance"].store( A )
2714     if selfA._toStore("JacobianMatrixAtOptimum"):
2715         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
2716     if selfA._toStore("KalmanGainAtOptimum"):
2717         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
2718         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
2719         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
2720     #
2721     # Calculs et/ou stockages supplémentaires
2722     # ---------------------------------------
2723     if selfA._toStore("Innovation") or \
2724         selfA._toStore("SigmaObs2") or \
2725         selfA._toStore("MahalanobisConsistency") or \
2726         selfA._toStore("OMB"):
2727         d  = Y - HXb
2728     if selfA._toStore("Innovation"):
2729         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
2730     if selfA._toStore("BMA"):
2731         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
2732     if selfA._toStore("OMA"):
2733         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
2734     if selfA._toStore("OMB"):
2735         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
2736     if selfA._toStore("SigmaObs2"):
2737         TraceR = R.trace(Y.size)
2738         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
2739     if selfA._toStore("MahalanobisConsistency"):
2740         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
2741     if selfA._toStore("SimulationQuantiles"):
2742         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
2743     if selfA._toStore("SimulatedObservationAtBackground"):
2744         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
2745     if selfA._toStore("SimulatedObservationAtOptimum"):
2746         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
2747     #
2748     return 0
2749
2750 # ==============================================================================
2751 def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16"):
2752     """
2753     Stochastic EnKF
2754     """
2755     if selfA._parameters["EstimationOf"] == "Parameters":
2756         selfA._parameters["StoreInternalVariables"] = True
2757     #
2758     # Opérateurs
2759     H = HO["Direct"].appliedControledFormTo
2760     #
2761     if selfA._parameters["EstimationOf"] == "State":
2762         M = EM["Direct"].appliedControledFormTo
2763     #
2764     if CM is not None and "Tangent" in CM and U is not None:
2765         Cm = CM["Tangent"].asMatrix(Xb)
2766     else:
2767         Cm = None
2768     #
2769     # Durée d'observation et tailles
2770     if hasattr(Y,"stepnumber"):
2771         duration = Y.stepnumber()
2772         __p = numpy.cumprod(Y.shape())[-1]
2773     else:
2774         duration = 2
2775         __p = numpy.array(Y).size
2776     #
2777     # Précalcul des inversions de B et R
2778     if selfA._parameters["StoreInternalVariables"] \
2779         or selfA._toStore("CostFunctionJ") \
2780         or selfA._toStore("CostFunctionJb") \
2781         or selfA._toStore("CostFunctionJo") \
2782         or selfA._toStore("CurrentOptimum") \
2783         or selfA._toStore("APosterioriCovariance"):
2784         BI = B.getI()
2785         RI = R.getI()
2786     #
2787     __n = Xb.size
2788     __m = selfA._parameters["NumberOfMembers"]
2789     #
2790     if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
2791     else:                         Rn = R
2792     #
2793     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
2794         Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
2795         selfA.StoredVariables["Analysis"].store( Xb )
2796         if selfA._toStore("APosterioriCovariance"):
2797             if hasattr(B,"asfullmatrix"):
2798                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
2799             else:
2800                 selfA.StoredVariables["APosterioriCovariance"].store( B )
2801         selfA._setInternalState("seed", numpy.random.get_state())
2802     elif selfA._parameters["nextStep"]:
2803         Xn = selfA._getInternalState("Xn")
2804     #
2805     previousJMinimum = numpy.finfo(float).max
2806     #
2807     for step in range(duration-1):
2808         numpy.random.set_state(selfA._getInternalState("seed"))
2809         if hasattr(Y,"store"):
2810             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
2811         else:
2812             Ynpu = numpy.ravel( Y ).reshape((__p,1))
2813         #
2814         if U is not None:
2815             if hasattr(U,"store") and len(U)>1:
2816                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
2817             elif hasattr(U,"store") and len(U)==1:
2818                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
2819             else:
2820                 Un = numpy.asmatrix(numpy.ravel( U )).T
2821         else:
2822             Un = None
2823         #
2824         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
2825             Xn = CovarianceInflation( Xn,
2826                 selfA._parameters["InflationType"],
2827                 selfA._parameters["InflationFactor"],
2828                 )
2829         #
2830         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
2831             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
2832                 argsAsSerie = True,
2833                 returnSerieAsArrayMatrix = True )
2834             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
2835             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
2836                 argsAsSerie = True,
2837                 returnSerieAsArrayMatrix = True )
2838             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
2839                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
2840                 Xn_predicted = Xn_predicted + Cm * Un
2841         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
2842             # --- > Par principe, M = Id, Q = 0
2843             Xn_predicted = EMX = Xn
2844             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
2845                 argsAsSerie = True,
2846                 returnSerieAsArrayMatrix = True )
2847         #
2848         # Mean of forecast and observation of forecast
2849         Xfm  = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
2850         Hfm  = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
2851         #
2852         #--------------------------
2853         if VariantM == "KalmanFilterFormula05":
2854             PfHT, HPfHT = 0., 0.
2855             for i in range(__m):
2856                 Exfi = Xn_predicted[:,i].reshape((__n,1)) - Xfm
2857                 Eyfi = HX_predicted[:,i].reshape((__p,1)) - Hfm
2858                 PfHT  += Exfi * Eyfi.T
2859                 HPfHT += Eyfi * Eyfi.T
2860             PfHT  = (1./(__m-1)) * PfHT
2861             HPfHT = (1./(__m-1)) * HPfHT
2862             Kn     = PfHT * ( R + HPfHT ).I
2863             del PfHT, HPfHT
2864             #
2865             for i in range(__m):
2866                 ri = numpy.random.multivariate_normal(numpy.zeros(__p), Rn)
2867                 Xn[:,i] = numpy.ravel(Xn_predicted[:,i]) + Kn @ (numpy.ravel(Ynpu) + ri - HX_predicted[:,i])
2868         #--------------------------
2869         elif VariantM == "KalmanFilterFormula16":
2870             EpY   = EnsembleOfCenteredPerturbations(Ynpu, Rn, __m)
2871             EpYm  = EpY.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
2872             #
2873             EaX   = EnsembleOfAnomalies( Xn_predicted ) / math.sqrt(__m-1)
2874             EaY = (HX_predicted - Hfm - EpY + EpYm) / math.sqrt(__m-1)
2875             #
2876             Kn = EaX @ EaY.T @ numpy.linalg.inv( EaY @ EaY.T)
2877             #
2878             for i in range(__m):
2879                 Xn[:,i] = numpy.ravel(Xn_predicted[:,i]) + Kn @ (numpy.ravel(EpY[:,i]) - HX_predicted[:,i])
2880         #--------------------------
2881         else:
2882             raise ValueError("VariantM has to be chosen in the authorized methods list.")
2883         #
2884         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
2885             Xn = CovarianceInflation( Xn,
2886                 selfA._parameters["InflationType"],
2887                 selfA._parameters["InflationFactor"],
2888                 )
2889         #
2890         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
2891         #--------------------------
2892         selfA._setInternalState("Xn", Xn)
2893         selfA._setInternalState("seed", numpy.random.get_state())
2894         #--------------------------
2895         #
2896         if selfA._parameters["StoreInternalVariables"] \
2897             or selfA._toStore("CostFunctionJ") \
2898             or selfA._toStore("CostFunctionJb") \
2899             or selfA._toStore("CostFunctionJo") \
2900             or selfA._toStore("APosterioriCovariance") \
2901             or selfA._toStore("InnovationAtCurrentAnalysis") \
2902             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
2903             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2904             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
2905             _Innovation = Ynpu - _HXa
2906         #
2907         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2908         # ---> avec analysis
2909         selfA.StoredVariables["Analysis"].store( Xa )
2910         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
2911             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
2912         if selfA._toStore("InnovationAtCurrentAnalysis"):
2913             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
2914         # ---> avec current state
2915         if selfA._parameters["StoreInternalVariables"] \
2916             or selfA._toStore("CurrentState"):
2917             selfA.StoredVariables["CurrentState"].store( Xn )
2918         if selfA._toStore("ForecastState"):
2919             selfA.StoredVariables["ForecastState"].store( EMX )
2920         if selfA._toStore("ForecastCovariance"):
2921             selfA.StoredVariables["ForecastCovariance"].store( EnsembleErrorCovariance(EMX) )
2922         if selfA._toStore("BMA"):
2923             selfA.StoredVariables["BMA"].store( EMX - Xa )
2924         if selfA._toStore("InnovationAtCurrentState"):
2925             selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
2926         if selfA._toStore("SimulatedObservationAtCurrentState") \
2927             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2928             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
2929         # ---> autres
2930         if selfA._parameters["StoreInternalVariables"] \
2931             or selfA._toStore("CostFunctionJ") \
2932             or selfA._toStore("CostFunctionJb") \
2933             or selfA._toStore("CostFunctionJo") \
2934             or selfA._toStore("CurrentOptimum") \
2935             or selfA._toStore("APosterioriCovariance"):
2936             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
2937             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
2938             J   = Jb + Jo
2939             selfA.StoredVariables["CostFunctionJb"].store( Jb )
2940             selfA.StoredVariables["CostFunctionJo"].store( Jo )
2941             selfA.StoredVariables["CostFunctionJ" ].store( J )
2942             #
2943             if selfA._toStore("IndexOfOptimum") \
2944                 or selfA._toStore("CurrentOptimum") \
2945                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
2946                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
2947                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
2948                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2949                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2950             if selfA._toStore("IndexOfOptimum"):
2951                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2952             if selfA._toStore("CurrentOptimum"):
2953                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
2954             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2955                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
2956             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2957                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2958             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2959                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2960             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2961                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2962         if selfA._toStore("APosterioriCovariance"):
2963             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
2964         if selfA._parameters["EstimationOf"] == "Parameters" \
2965             and J < previousJMinimum:
2966             previousJMinimum    = J
2967             XaMin               = Xa
2968             if selfA._toStore("APosterioriCovariance"):
2969                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
2970     #
2971     # Stockage final supplémentaire de l'optimum en estimation de paramètres
2972     # ----------------------------------------------------------------------
2973     if selfA._parameters["EstimationOf"] == "Parameters":
2974         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2975         selfA.StoredVariables["Analysis"].store( XaMin )
2976         if selfA._toStore("APosterioriCovariance"):
2977             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
2978         if selfA._toStore("BMA"):
2979             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
2980     #
2981     return 0
2982
2983 # ==============================================================================
2984 def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
2985     """
2986     3DVAR
2987     """
2988     #
2989     # Initialisations
2990     # ---------------
2991     #
2992     # Opérateurs
2993     Hm = HO["Direct"].appliedTo
2994     Ha = HO["Adjoint"].appliedInXTo
2995     #
2996     # Utilisation éventuelle d'un vecteur H(Xb) précalculé
2997     if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
2998         HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
2999     else:
3000         HXb = Hm( Xb )
3001     HXb = numpy.asmatrix(numpy.ravel( HXb )).T
3002     if Y.size != HXb.size:
3003         raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
3004     if max(Y.shape) != max(HXb.shape):
3005         raise ValueError("The shapes %s of observations Y and %s of observed calculation H(X) are different, they have to be identical."%(Y.shape,HXb.shape))
3006     #
3007     if selfA._toStore("JacobianMatrixAtBackground"):
3008         HtMb = HO["Tangent"].asMatrix(ValueForMethodForm = Xb)
3009         HtMb = HtMb.reshape(Y.size,Xb.size) # ADAO & check shape
3010         selfA.StoredVariables["JacobianMatrixAtBackground"].store( HtMb )
3011     #
3012     # Précalcul des inversions de B et R
3013     BI = B.getI()
3014     RI = R.getI()
3015     #
3016     # Point de démarrage de l'optimisation
3017     Xini = selfA._parameters["InitializationPoint"]
3018     #
3019     # Définition de la fonction-coût
3020     # ------------------------------
3021     def CostFunction(x):
3022         _X  = numpy.asmatrix(numpy.ravel( x )).T
3023         if selfA._parameters["StoreInternalVariables"] or \
3024             selfA._toStore("CurrentState") or \
3025             selfA._toStore("CurrentOptimum"):
3026             selfA.StoredVariables["CurrentState"].store( _X )
3027         _HX = Hm( _X )
3028         _HX = numpy.asmatrix(numpy.ravel( _HX )).T
3029         _Innovation = Y - _HX
3030         if selfA._toStore("SimulatedObservationAtCurrentState") or \
3031             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3032             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
3033         if selfA._toStore("InnovationAtCurrentState"):
3034             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
3035         #
3036         Jb  = float( 0.5 * (_X - Xb).T * BI * (_X - Xb) )
3037         Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
3038         J   = Jb + Jo
3039         #
3040         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
3041         selfA.StoredVariables["CostFunctionJb"].store( Jb )
3042         selfA.StoredVariables["CostFunctionJo"].store( Jo )
3043         selfA.StoredVariables["CostFunctionJ" ].store( J )
3044         if selfA._toStore("IndexOfOptimum") or \
3045             selfA._toStore("CurrentOptimum") or \
3046             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
3047             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
3048             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
3049             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3050             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3051         if selfA._toStore("IndexOfOptimum"):
3052             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3053         if selfA._toStore("CurrentOptimum"):
3054             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
3055         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3056             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
3057         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3058             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3059         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3060             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3061         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3062             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3063         return J
3064     #
3065     def GradientOfCostFunction(x):
3066         _X      = numpy.asmatrix(numpy.ravel( x )).T
3067         _HX     = Hm( _X )
3068         _HX     = numpy.asmatrix(numpy.ravel( _HX )).T
3069         GradJb  = BI * (_X - Xb)
3070         GradJo  = - Ha( (_X, RI * (Y - _HX)) )
3071         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
3072         return GradJ
3073     #
3074     # Minimisation de la fonctionnelle
3075     # --------------------------------
3076     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
3077     #
3078     if selfA._parameters["Minimizer"] == "LBFGSB":
3079         if "0.19" <= scipy.version.version <= "1.1.0":
3080             import lbfgsbhlt as optimiseur
3081         else:
3082             import scipy.optimize as optimiseur
3083         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
3084             func        = CostFunction,
3085             x0          = Xini,
3086             fprime      = GradientOfCostFunction,
3087             args        = (),
3088             bounds      = selfA._parameters["Bounds"],
3089             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
3090             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
3091             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3092             iprint      = selfA._parameters["optiprint"],
3093             )
3094         nfeval = Informations['funcalls']
3095         rc     = Informations['warnflag']
3096     elif selfA._parameters["Minimizer"] == "TNC":
3097         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
3098             func        = CostFunction,
3099             x0          = Xini,
3100             fprime      = GradientOfCostFunction,
3101             args        = (),
3102             bounds      = selfA._parameters["Bounds"],
3103             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
3104             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3105             ftol        = selfA._parameters["CostDecrementTolerance"],
3106             messages    = selfA._parameters["optmessages"],
3107             )
3108     elif selfA._parameters["Minimizer"] == "CG":
3109         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
3110             f           = CostFunction,
3111             x0          = Xini,
3112             fprime      = GradientOfCostFunction,
3113             args        = (),
3114             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3115             gtol        = selfA._parameters["GradientNormTolerance"],
3116             disp        = selfA._parameters["optdisp"],
3117             full_output = True,
3118             )
3119     elif selfA._parameters["Minimizer"] == "NCG":
3120         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
3121             f           = CostFunction,
3122             x0          = Xini,
3123             fprime      = GradientOfCostFunction,
3124             args        = (),
3125             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3126             avextol     = selfA._parameters["CostDecrementTolerance"],
3127             disp        = selfA._parameters["optdisp"],
3128             full_output = True,
3129             )
3130     elif selfA._parameters["Minimizer"] == "BFGS":
3131         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
3132             f           = CostFunction,
3133             x0          = Xini,
3134             fprime      = GradientOfCostFunction,
3135             args        = (),
3136             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3137             gtol        = selfA._parameters["GradientNormTolerance"],
3138             disp        = selfA._parameters["optdisp"],
3139             full_output = True,
3140             )
3141     else:
3142         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
3143     #
3144     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3145     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
3146     #
3147     # Correction pour pallier a un bug de TNC sur le retour du Minimum
3148     # ----------------------------------------------------------------
3149     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
3150         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
3151     #
3152     # Obtention de l'analyse
3153     # ----------------------
3154     Xa = numpy.asmatrix(numpy.ravel( Minimum )).T
3155     #
3156     selfA.StoredVariables["Analysis"].store( Xa )
3157     #
3158     if selfA._toStore("OMA") or \
3159         selfA._toStore("SigmaObs2") or \
3160         selfA._toStore("SimulationQuantiles") or \
3161         selfA._toStore("SimulatedObservationAtOptimum"):
3162         if selfA._toStore("SimulatedObservationAtCurrentState"):
3163             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
3164         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3165             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
3166         else:
3167             HXa = Hm( Xa )
3168     #
3169     # Calcul de la covariance d'analyse
3170     # ---------------------------------
3171     if selfA._toStore("APosterioriCovariance") or \
3172         selfA._toStore("SimulationQuantiles") or \
3173         selfA._toStore("JacobianMatrixAtOptimum") or \
3174         selfA._toStore("KalmanGainAtOptimum"):
3175         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
3176         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
3177     if selfA._toStore("APosterioriCovariance") or \
3178         selfA._toStore("SimulationQuantiles") or \
3179         selfA._toStore("KalmanGainAtOptimum"):
3180         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
3181         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
3182     if selfA._toStore("APosterioriCovariance") or \
3183         selfA._toStore("SimulationQuantiles"):
3184         HessienneI = []
3185         nb = Xa.size
3186         for i in range(nb):
3187             _ee    = numpy.matrix(numpy.zeros(nb)).T
3188             _ee[i] = 1.
3189             _HtEE  = numpy.dot(HtM,_ee)
3190             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
3191             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
3192         HessienneI = numpy.matrix( HessienneI )
3193         A = HessienneI.I
3194         if min(A.shape) != max(A.shape):
3195             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)))
3196         if (numpy.diag(A) < 0).any():
3197             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,))
3198         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
3199             try:
3200                 L = numpy.linalg.cholesky( A )
3201             except:
3202                 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,))
3203     if selfA._toStore("APosterioriCovariance"):
3204         selfA.StoredVariables["APosterioriCovariance"].store( A )
3205     if selfA._toStore("JacobianMatrixAtOptimum"):
3206         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
3207     if selfA._toStore("KalmanGainAtOptimum"):
3208         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
3209         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
3210         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
3211     #
3212     # Calculs et/ou stockages supplémentaires
3213     # ---------------------------------------
3214     if selfA._toStore("Innovation") or \
3215         selfA._toStore("SigmaObs2") or \
3216         selfA._toStore("MahalanobisConsistency") or \
3217         selfA._toStore("OMB"):
3218         d  = Y - HXb
3219     if selfA._toStore("Innovation"):
3220         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
3221     if selfA._toStore("BMA"):
3222         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
3223     if selfA._toStore("OMA"):
3224         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
3225     if selfA._toStore("OMB"):
3226         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
3227     if selfA._toStore("SigmaObs2"):
3228         TraceR = R.trace(Y.size)
3229         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
3230     if selfA._toStore("MahalanobisConsistency"):
3231         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
3232     if selfA._toStore("SimulationQuantiles"):
3233         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
3234     if selfA._toStore("SimulatedObservationAtBackground"):
3235         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
3236     if selfA._toStore("SimulatedObservationAtOptimum"):
3237         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
3238     #
3239     return 0
3240
3241 # ==============================================================================
3242 def std4dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
3243     """
3244     4DVAR
3245     """
3246     #
3247     # Initialisations
3248     # ---------------
3249     #
3250     # Opérateurs
3251     Hm = HO["Direct"].appliedControledFormTo
3252     Mm = EM["Direct"].appliedControledFormTo
3253     #
3254     if CM is not None and "Tangent" in CM and U is not None:
3255         Cm = CM["Tangent"].asMatrix(Xb)
3256     else:
3257         Cm = None
3258     #
3259     def Un(_step):
3260         if U is not None:
3261             if hasattr(U,"store") and 1<=_step<len(U) :
3262                 _Un = numpy.asmatrix(numpy.ravel( U[_step] )).T
3263             elif hasattr(U,"store") and len(U)==1:
3264                 _Un = numpy.asmatrix(numpy.ravel( U[0] )).T
3265             else:
3266                 _Un = numpy.asmatrix(numpy.ravel( U )).T
3267         else:
3268             _Un = None
3269         return _Un
3270     def CmUn(_xn,_un):
3271         if Cm is not None and _un is not None: # Attention : si Cm est aussi dans M, doublon !
3272             _Cm   = Cm.reshape(_xn.size,_un.size) # ADAO & check shape
3273             _CmUn = _Cm * _un
3274         else:
3275             _CmUn = 0.
3276         return _CmUn
3277     #
3278     # Remarque : les observations sont exploitées à partir du pas de temps
3279     # numéro 1, et sont utilisées dans Yo comme rangées selon ces indices.
3280     # Donc le pas 0 n'est pas utilisé puisque la première étape commence
3281     # avec l'observation du pas 1.
3282     #
3283     # Nombre de pas identique au nombre de pas d'observations
3284     if hasattr(Y,"stepnumber"):
3285         duration = Y.stepnumber()
3286     else:
3287         duration = 2
3288     #
3289     # Précalcul des inversions de B et R
3290     BI = B.getI()
3291     RI = R.getI()
3292     #
3293     # Point de démarrage de l'optimisation
3294     Xini = selfA._parameters["InitializationPoint"]
3295     #
3296     # Définition de la fonction-coût
3297     # ------------------------------
3298     selfA.DirectCalculation = [None,] # Le pas 0 n'est pas observé
3299     selfA.DirectInnovation  = [None,] # Le pas 0 n'est pas observé
3300     def CostFunction(x):
3301         _X  = numpy.asmatrix(numpy.ravel( x )).T
3302         if selfA._parameters["StoreInternalVariables"] or \
3303             selfA._toStore("CurrentState") or \
3304             selfA._toStore("CurrentOptimum"):
3305             selfA.StoredVariables["CurrentState"].store( _X )
3306         Jb  = float( 0.5 * (_X - Xb).T * BI * (_X - Xb) )
3307         selfA.DirectCalculation = [None,]
3308         selfA.DirectInnovation  = [None,]
3309         Jo  = 0.
3310         _Xn = _X
3311         for step in range(0,duration-1):
3312             if hasattr(Y,"store"):
3313                 _Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
3314             else:
3315                 _Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
3316             _Un = Un(step)
3317             #
3318             # Etape d'évolution
3319             if selfA._parameters["EstimationOf"] == "State":
3320                 _Xn = Mm( (_Xn, _Un) ) + CmUn(_Xn, _Un)
3321             elif selfA._parameters["EstimationOf"] == "Parameters":
3322                 pass
3323             #
3324             if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3325                 __Bounds = ForceNumericBounds( selfA._parameters["Bounds"] )
3326                 _Xn = numpy.max(numpy.hstack((_Xn,numpy.asmatrix(__Bounds)[:,0])),axis=1)
3327                 _Xn = numpy.min(numpy.hstack((_Xn,numpy.asmatrix(__Bounds)[:,1])),axis=1)
3328             #
3329             # Etape de différence aux observations
3330             if selfA._parameters["EstimationOf"] == "State":
3331                 _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, None) ) )).T
3332             elif selfA._parameters["EstimationOf"] == "Parameters":
3333                 _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, _Un) ) )).T - CmUn(_Xn, _Un)
3334             #
3335             # Stockage de l'état
3336             selfA.DirectCalculation.append( _Xn )
3337             selfA.DirectInnovation.append( _YmHMX )
3338             #
3339             # Ajout dans la fonctionnelle d'observation
3340             Jo = Jo + 0.5 * float( _YmHMX.T * RI * _YmHMX )
3341         J = Jb + Jo
3342         #
3343         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
3344         selfA.StoredVariables["CostFunctionJb"].store( Jb )
3345         selfA.StoredVariables["CostFunctionJo"].store( Jo )
3346         selfA.StoredVariables["CostFunctionJ" ].store( J )
3347         if selfA._toStore("IndexOfOptimum") or \
3348             selfA._toStore("CurrentOptimum") or \
3349             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
3350             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
3351             selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3352             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3353         if selfA._toStore("IndexOfOptimum"):
3354             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3355         if selfA._toStore("CurrentOptimum"):
3356             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
3357         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3358             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3359         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3360             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3361         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3362             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3363         return J
3364     #
3365     def GradientOfCostFunction(x):
3366         _X      = numpy.asmatrix(numpy.ravel( x )).T
3367         GradJb  = BI * (_X - Xb)
3368         GradJo  = 0.
3369         for step in range(duration-1,0,-1):
3370             # Étape de récupération du dernier stockage de l'évolution
3371             _Xn = selfA.DirectCalculation.pop()
3372             # Étape de récupération du dernier stockage de l'innovation
3373             _YmHMX = selfA.DirectInnovation.pop()
3374             # Calcul des adjoints
3375             Ha = HO["Adjoint"].asMatrix(ValueForMethodForm = _Xn)
3376             Ha = Ha.reshape(_Xn.size,_YmHMX.size) # ADAO & check shape
3377             Ma = EM["Adjoint"].asMatrix(ValueForMethodForm = _Xn)
3378             Ma = Ma.reshape(_Xn.size,_Xn.size) # ADAO & check shape
3379             # Calcul du gradient par état adjoint
3380             GradJo = GradJo + Ha * RI * _YmHMX # Équivaut pour Ha linéaire à : Ha( (_Xn, RI * _YmHMX) )
3381             GradJo = Ma * GradJo               # Équivaut pour Ma linéaire à : Ma( (_Xn, GradJo) )
3382         GradJ = numpy.ravel( GradJb ) - numpy.ravel( GradJo )
3383         return GradJ
3384     #
3385     # Minimisation de la fonctionnelle
3386     # --------------------------------
3387     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
3388     #
3389     if selfA._parameters["Minimizer"] == "LBFGSB":
3390         if "0.19" <= scipy.version.version <= "1.1.0":
3391             import lbfgsbhlt as optimiseur
3392         else:
3393             import scipy.optimize as optimiseur
3394         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
3395             func        = CostFunction,
3396             x0          = Xini,
3397             fprime      = GradientOfCostFunction,
3398             args        = (),
3399             bounds      = selfA._parameters["Bounds"],
3400             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
3401             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
3402             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3403             iprint      = selfA._parameters["optiprint"],
3404             )
3405         nfeval = Informations['funcalls']
3406         rc     = Informations['warnflag']
3407     elif selfA._parameters["Minimizer"] == "TNC":
3408         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
3409             func        = CostFunction,
3410             x0          = Xini,
3411             fprime      = GradientOfCostFunction,
3412             args        = (),
3413             bounds      = selfA._parameters["Bounds"],
3414             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
3415             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3416             ftol        = selfA._parameters["CostDecrementTolerance"],
3417             messages    = selfA._parameters["optmessages"],
3418             )
3419     elif selfA._parameters["Minimizer"] == "CG":
3420         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
3421             f           = CostFunction,
3422             x0          = Xini,
3423             fprime      = GradientOfCostFunction,
3424             args        = (),
3425             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3426             gtol        = selfA._parameters["GradientNormTolerance"],
3427             disp        = selfA._parameters["optdisp"],
3428             full_output = True,
3429             )
3430     elif selfA._parameters["Minimizer"] == "NCG":
3431         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
3432             f           = CostFunction,
3433             x0          = Xini,
3434             fprime      = GradientOfCostFunction,
3435             args        = (),
3436             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3437             avextol     = selfA._parameters["CostDecrementTolerance"],
3438             disp        = selfA._parameters["optdisp"],
3439             full_output = True,
3440             )
3441     elif selfA._parameters["Minimizer"] == "BFGS":
3442         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
3443             f           = CostFunction,
3444             x0          = Xini,
3445             fprime      = GradientOfCostFunction,
3446             args        = (),
3447             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3448             gtol        = selfA._parameters["GradientNormTolerance"],
3449             disp        = selfA._parameters["optdisp"],
3450             full_output = True,
3451             )
3452     else:
3453         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
3454     #
3455     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3456     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
3457     #
3458     # Correction pour pallier a un bug de TNC sur le retour du Minimum
3459     # ----------------------------------------------------------------
3460     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
3461         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
3462     #
3463     # Obtention de l'analyse
3464     # ----------------------
3465     Xa = numpy.asmatrix(numpy.ravel( Minimum )).T
3466     #
3467     selfA.StoredVariables["Analysis"].store( Xa )
3468     #
3469     # Calculs et/ou stockages supplémentaires
3470     # ---------------------------------------
3471     if selfA._toStore("BMA"):
3472         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
3473     #
3474     return 0
3475
3476 # ==============================================================================
3477 def stdkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
3478     """
3479     Standard Kalman Filter
3480     """
3481     if selfA._parameters["EstimationOf"] == "Parameters":
3482         selfA._parameters["StoreInternalVariables"] = True
3483     #
3484     # Opérateurs
3485     # ----------
3486     Ht = HO["Tangent"].asMatrix(Xb)
3487     Ha = HO["Adjoint"].asMatrix(Xb)
3488     #
3489     if selfA._parameters["EstimationOf"] == "State":
3490         Mt = EM["Tangent"].asMatrix(Xb)
3491         Ma = EM["Adjoint"].asMatrix(Xb)
3492     #
3493     if CM is not None and "Tangent" in CM and U is not None:
3494         Cm = CM["Tangent"].asMatrix(Xb)
3495     else:
3496         Cm = None
3497     #
3498     # Durée d'observation et tailles
3499     if hasattr(Y,"stepnumber"):
3500         duration = Y.stepnumber()
3501         __p = numpy.cumprod(Y.shape())[-1]
3502     else:
3503         duration = 2
3504         __p = numpy.array(Y).size
3505     #
3506     # Précalcul des inversions de B et R
3507     if selfA._parameters["StoreInternalVariables"] \
3508         or selfA._toStore("CostFunctionJ") \
3509         or selfA._toStore("CostFunctionJb") \
3510         or selfA._toStore("CostFunctionJo") \
3511         or selfA._toStore("CurrentOptimum") \
3512         or selfA._toStore("APosterioriCovariance"):
3513         BI = B.getI()
3514         RI = R.getI()
3515     #
3516     __n = Xb.size
3517     #
3518     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
3519         Xn = Xb
3520         Pn = B
3521         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3522         selfA.StoredVariables["Analysis"].store( Xb )
3523         if selfA._toStore("APosterioriCovariance"):
3524             if hasattr(B,"asfullmatrix"):
3525                 selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(__n) )
3526             else:
3527                 selfA.StoredVariables["APosterioriCovariance"].store( B )
3528         selfA._setInternalState("seed", numpy.random.get_state())
3529     elif selfA._parameters["nextStep"]:
3530         Xn = selfA._getInternalState("Xn")
3531         Pn = selfA._getInternalState("Pn")
3532     #
3533     if selfA._parameters["EstimationOf"] == "Parameters":
3534         XaMin            = Xn
3535         previousJMinimum = numpy.finfo(float).max
3536     #
3537     for step in range(duration-1):
3538         if hasattr(Y,"store"):
3539             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
3540         else:
3541             Ynpu = numpy.ravel( Y ).reshape((__p,1))
3542         #
3543         if U is not None:
3544             if hasattr(U,"store") and len(U)>1:
3545                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
3546             elif hasattr(U,"store") and len(U)==1:
3547                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
3548             else:
3549                 Un = numpy.asmatrix(numpy.ravel( U )).T
3550         else:
3551             Un = None
3552         #
3553         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
3554             Xn_predicted = Mt * Xn
3555             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
3556                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
3557                 Xn_predicted = Xn_predicted + Cm * Un
3558             Pn_predicted = Q + Mt * Pn * Ma
3559         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
3560             # --- > Par principe, M = Id, Q = 0
3561             Xn_predicted = Xn
3562             Pn_predicted = Pn
3563         #
3564         if selfA._parameters["EstimationOf"] == "State":
3565             HX_predicted = Ht * Xn_predicted
3566             _Innovation  = Ynpu - HX_predicted
3567         elif selfA._parameters["EstimationOf"] == "Parameters":
3568             HX_predicted = Ht * Xn_predicted
3569             _Innovation  = Ynpu - HX_predicted
3570             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
3571                 _Innovation = _Innovation - Cm * Un
3572         #
3573         Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha))
3574         Xn = Xn_predicted + Kn * _Innovation
3575         Pn = Pn_predicted - Kn * Ht * Pn_predicted
3576         #
3577         Xa = Xn # Pointeurs
3578         #--------------------------
3579         selfA._setInternalState("Xn", Xn)
3580         selfA._setInternalState("Pn", Pn)
3581         #--------------------------
3582         #
3583         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3584         # ---> avec analysis
3585         selfA.StoredVariables["Analysis"].store( Xa )
3586         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
3587             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( Ht * Xa )
3588         if selfA._toStore("InnovationAtCurrentAnalysis"):
3589             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
3590         # ---> avec current state
3591         if selfA._parameters["StoreInternalVariables"] \
3592             or selfA._toStore("CurrentState"):
3593             selfA.StoredVariables["CurrentState"].store( Xn )
3594         if selfA._toStore("ForecastState"):
3595             selfA.StoredVariables["ForecastState"].store( Xn_predicted )
3596         if selfA._toStore("ForecastCovariance"):
3597             selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
3598         if selfA._toStore("BMA"):
3599             selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
3600         if selfA._toStore("InnovationAtCurrentState"):
3601             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
3602         if selfA._toStore("SimulatedObservationAtCurrentState") \
3603             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3604             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
3605         # ---> autres
3606         if selfA._parameters["StoreInternalVariables"] \
3607             or selfA._toStore("CostFunctionJ") \
3608             or selfA._toStore("CostFunctionJb") \
3609             or selfA._toStore("CostFunctionJo") \
3610             or selfA._toStore("CurrentOptimum") \
3611             or selfA._toStore("APosterioriCovariance"):
3612             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
3613             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
3614             J   = Jb + Jo
3615             selfA.StoredVariables["CostFunctionJb"].store( Jb )
3616             selfA.StoredVariables["CostFunctionJo"].store( Jo )
3617             selfA.StoredVariables["CostFunctionJ" ].store( J )
3618             #
3619             if selfA._toStore("IndexOfOptimum") \
3620                 or selfA._toStore("CurrentOptimum") \
3621                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
3622                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
3623                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
3624                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3625                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3626             if selfA._toStore("IndexOfOptimum"):
3627                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3628             if selfA._toStore("CurrentOptimum"):
3629                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
3630             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3631                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
3632             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3633                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3634             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3635                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3636             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3637                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3638         if selfA._toStore("APosterioriCovariance"):
3639             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
3640         if selfA._parameters["EstimationOf"] == "Parameters" \
3641             and J < previousJMinimum:
3642             previousJMinimum    = J
3643             XaMin               = Xa
3644             if selfA._toStore("APosterioriCovariance"):
3645                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
3646     #
3647     # Stockage final supplémentaire de l'optimum en estimation de paramètres
3648     # ----------------------------------------------------------------------
3649     if selfA._parameters["EstimationOf"] == "Parameters":
3650         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3651         selfA.StoredVariables["Analysis"].store( XaMin )
3652         if selfA._toStore("APosterioriCovariance"):
3653             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
3654         if selfA._toStore("BMA"):
3655             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
3656     #
3657     return 0
3658
3659 # ==============================================================================
3660 def uckf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
3661     """
3662     Unscented Kalman Filter
3663     """
3664     if selfA._parameters["EstimationOf"] == "Parameters":
3665         selfA._parameters["StoreInternalVariables"] = True
3666     selfA._parameters["Bounds"] = ForceNumericBounds( selfA._parameters["Bounds"] )
3667     #
3668     L     = Xb.size
3669     Alpha = selfA._parameters["Alpha"]
3670     Beta  = selfA._parameters["Beta"]
3671     if selfA._parameters["Kappa"] == 0:
3672         if selfA._parameters["EstimationOf"] == "State":
3673             Kappa = 0
3674         elif selfA._parameters["EstimationOf"] == "Parameters":
3675             Kappa = 3 - L
3676     else:
3677         Kappa = selfA._parameters["Kappa"]
3678     Lambda = float( Alpha**2 ) * ( L + Kappa ) - L
3679     Gamma  = math.sqrt( L + Lambda )
3680     #
3681     Ww = []
3682     Ww.append( 0. )
3683     for i in range(2*L):
3684         Ww.append( 1. / (2.*(L + Lambda)) )
3685     #
3686     Wm = numpy.array( Ww )
3687     Wm[0] = Lambda / (L + Lambda)
3688     Wc = numpy.array( Ww )
3689     Wc[0] = Lambda / (L + Lambda) + (1. - Alpha**2 + Beta)
3690     #
3691     # Opérateurs
3692     Hm = HO["Direct"].appliedControledFormTo
3693     #
3694     if selfA._parameters["EstimationOf"] == "State":
3695         Mm = EM["Direct"].appliedControledFormTo
3696     #
3697     if CM is not None and "Tangent" in CM and U is not None:
3698         Cm = CM["Tangent"].asMatrix(Xb)
3699     else:
3700         Cm = None
3701     #
3702     # Durée d'observation et tailles
3703     if hasattr(Y,"stepnumber"):
3704         duration = Y.stepnumber()
3705         __p = numpy.cumprod(Y.shape())[-1]
3706     else:
3707         duration = 2
3708         __p = numpy.array(Y).size
3709     #
3710     # Précalcul des inversions de B et R
3711     if selfA._parameters["StoreInternalVariables"] \
3712         or selfA._toStore("CostFunctionJ") \
3713         or selfA._toStore("CostFunctionJb") \
3714         or selfA._toStore("CostFunctionJo") \
3715         or selfA._toStore("CurrentOptimum") \
3716         or selfA._toStore("APosterioriCovariance"):
3717         BI = B.getI()
3718         RI = R.getI()
3719     #
3720     __n = Xb.size
3721     #
3722     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
3723         Xn = Xb
3724         if hasattr(B,"asfullmatrix"):
3725             Pn = B.asfullmatrix(__n)
3726         else:
3727             Pn = B
3728         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3729         selfA.StoredVariables["Analysis"].store( Xb )
3730         if selfA._toStore("APosterioriCovariance"):
3731             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
3732     elif selfA._parameters["nextStep"]:
3733         Xn = selfA._getInternalState("Xn")
3734         Pn = selfA._getInternalState("Pn")
3735     #
3736     if selfA._parameters["EstimationOf"] == "Parameters":
3737         XaMin            = Xn
3738         previousJMinimum = numpy.finfo(float).max
3739     #
3740     for step in range(duration-1):
3741         if hasattr(Y,"store"):
3742             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
3743         else:
3744             Ynpu = numpy.ravel( Y ).reshape((__p,1))
3745         #
3746         if U is not None:
3747             if hasattr(U,"store") and len(U)>1:
3748                 Un = numpy.ravel( U[step] ).reshape((-1,1))
3749             elif hasattr(U,"store") and len(U)==1:
3750                 Un = numpy.ravel( U[0] ).reshape((-1,1))
3751             else:
3752                 Un = numpy.ravel( U ).reshape((-1,1))
3753         else:
3754             Un = None
3755         #
3756         Pndemi = numpy.linalg.cholesky(Pn)
3757         Xnp = numpy.hstack([Xn, Xn+Gamma*Pndemi, Xn-Gamma*Pndemi])
3758         nbSpts = 2*Xn.size+1
3759         #
3760         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3761             for point in range(nbSpts):
3762                 Xnp[:,point] = numpy.max(numpy.hstack((Xnp[:,point],numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
3763                 Xnp[:,point] = numpy.min(numpy.hstack((Xnp[:,point],numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
3764         #
3765         XEtnnp = []
3766         for point in range(nbSpts):
3767             if selfA._parameters["EstimationOf"] == "State":
3768                 XEtnnpi = numpy.asmatrix(numpy.ravel( Mm( (Xnp[:,point], Un) ) )).T
3769                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
3770                     Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
3771                     XEtnnpi = XEtnnpi + Cm * Un
3772                 if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3773                     XEtnnpi = numpy.max(numpy.hstack((XEtnnpi,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
3774                     XEtnnpi = numpy.min(numpy.hstack((XEtnnpi,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
3775             elif selfA._parameters["EstimationOf"] == "Parameters":
3776                 # --- > Par principe, M = Id, Q = 0
3777                 XEtnnpi = Xnp[:,point]
3778             XEtnnp.append( XEtnnpi )
3779         XEtnnp = numpy.hstack( XEtnnp )
3780         #
3781         Xncm = numpy.matrix( XEtnnp.getA()*numpy.array(Wm) ).sum(axis=1)
3782         #
3783         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3784             Xncm = numpy.max(numpy.hstack((Xncm,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
3785             Xncm = numpy.min(numpy.hstack((Xncm,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
3786         #
3787         if selfA._parameters["EstimationOf"] == "State":        Pnm = Q
3788         elif selfA._parameters["EstimationOf"] == "Parameters": Pnm = 0.
3789         for point in range(nbSpts):
3790             Pnm += Wc[i] * (XEtnnp[:,point]-Xncm) * (XEtnnp[:,point]-Xncm).T
3791         #
3792         if selfA._parameters["EstimationOf"] == "Parameters" and selfA._parameters["Bounds"] is not None:
3793             Pnmdemi = selfA._parameters["Reconditioner"] * numpy.linalg.cholesky(Pnm)
3794         else:
3795             Pnmdemi = numpy.linalg.cholesky(Pnm)
3796         #
3797         Xnnp = numpy.hstack([Xncm, Xncm+Gamma*Pnmdemi, Xncm-Gamma*Pnmdemi])
3798         #
3799         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3800             for point in range(nbSpts):
3801                 Xnnp[:,point] = numpy.max(numpy.hstack((Xnnp[:,point],numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
3802                 Xnnp[:,point] = numpy.min(numpy.hstack((Xnnp[:,point],numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
3803         #
3804         Ynnp = []
3805         for point in range(nbSpts):
3806             if selfA._parameters["EstimationOf"] == "State":
3807                 Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], None) ) )).T
3808             elif selfA._parameters["EstimationOf"] == "Parameters":
3809                 Ynnpi = numpy.asmatrix(numpy.ravel( Hm( (Xnnp[:,point], Un) ) )).T
3810             Ynnp.append( Ynnpi )
3811         Ynnp = numpy.hstack( Ynnp )
3812         #
3813         Yncm = numpy.matrix( Ynnp.getA()*numpy.array(Wm) ).sum(axis=1)
3814         #
3815         Pyyn = R
3816         Pxyn = 0.
3817         for point in range(nbSpts):
3818             Pyyn += Wc[i] * (Ynnp[:,point]-Yncm) * (Ynnp[:,point]-Yncm).T
3819             Pxyn += Wc[i] * (Xnnp[:,point]-Xncm) * (Ynnp[:,point]-Yncm).T
3820         #
3821         _Innovation  = Ynpu - Yncm
3822         if selfA._parameters["EstimationOf"] == "Parameters":
3823             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon !
3824                 _Innovation = _Innovation - Cm * Un
3825         #
3826         Kn = Pxyn * Pyyn.I
3827         Xn = Xncm + Kn * _Innovation
3828         Pn = Pnm - Kn * Pyyn * Kn.T
3829         #
3830         if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
3831             Xn = numpy.max(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
3832             Xn = numpy.min(numpy.hstack((Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
3833         #
3834         Xa = Xn # Pointeurs
3835         #--------------------------
3836         selfA._setInternalState("Xn", Xn)
3837         selfA._setInternalState("Pn", Pn)
3838         #--------------------------
3839         #
3840         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3841         # ---> avec analysis
3842         selfA.StoredVariables["Analysis"].store( Xa )
3843         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
3844             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( Hm((Xa, Un)) )
3845         if selfA._toStore("InnovationAtCurrentAnalysis"):
3846             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
3847         # ---> avec current state
3848         if selfA._parameters["StoreInternalVariables"] \
3849             or selfA._toStore("CurrentState"):
3850             selfA.StoredVariables["CurrentState"].store( Xn )
3851         if selfA._toStore("ForecastState"):
3852             selfA.StoredVariables["ForecastState"].store( Xncm )
3853         if selfA._toStore("ForecastCovariance"):
3854             selfA.StoredVariables["ForecastCovariance"].store( Pnm )
3855         if selfA._toStore("BMA"):
3856             selfA.StoredVariables["BMA"].store( Xncm - Xa )
3857         if selfA._toStore("InnovationAtCurrentState"):
3858             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
3859         if selfA._toStore("SimulatedObservationAtCurrentState") \
3860             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3861             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( Yncm )
3862         # ---> autres
3863         if selfA._parameters["StoreInternalVariables"] \
3864             or selfA._toStore("CostFunctionJ") \
3865             or selfA._toStore("CostFunctionJb") \
3866             or selfA._toStore("CostFunctionJo") \
3867             or selfA._toStore("CurrentOptimum") \
3868             or selfA._toStore("APosterioriCovariance"):
3869             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
3870             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
3871             J   = Jb + Jo
3872             selfA.StoredVariables["CostFunctionJb"].store( Jb )
3873             selfA.StoredVariables["CostFunctionJo"].store( Jo )
3874             selfA.StoredVariables["CostFunctionJ" ].store( J )
3875             #
3876             if selfA._toStore("IndexOfOptimum") \
3877                 or selfA._toStore("CurrentOptimum") \
3878                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
3879                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
3880                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
3881                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3882                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3883             if selfA._toStore("IndexOfOptimum"):
3884                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3885             if selfA._toStore("CurrentOptimum"):
3886                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
3887             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3888                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
3889             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3890                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3891             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3892                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3893             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3894                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3895         if selfA._toStore("APosterioriCovariance"):
3896             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
3897         if selfA._parameters["EstimationOf"] == "Parameters" \
3898             and J < previousJMinimum:
3899             previousJMinimum    = J
3900             XaMin               = Xa
3901             if selfA._toStore("APosterioriCovariance"):
3902                 covarianceXaMin = selfA.StoredVariables["APosterioriCovariance"][-1]
3903     #
3904     # Stockage final supplémentaire de l'optimum en estimation de paramètres
3905     # ----------------------------------------------------------------------
3906     if selfA._parameters["EstimationOf"] == "Parameters":
3907         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
3908         selfA.StoredVariables["Analysis"].store( XaMin )
3909         if selfA._toStore("APosterioriCovariance"):
3910             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
3911         if selfA._toStore("BMA"):
3912             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
3913     #
3914     return 0
3915
3916 # ==============================================================================
3917 def van3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
3918     """
3919     3DVAR variational analysis with no inversion of B
3920     """
3921     #
3922     # Initialisations
3923     # ---------------
3924     #
3925     # Opérateurs
3926     Hm = HO["Direct"].appliedTo
3927     Ha = HO["Adjoint"].appliedInXTo
3928     #
3929     # Précalcul des inversions de B et R
3930     BT = B.getT()
3931     RI = R.getI()
3932     #
3933     # Point de démarrage de l'optimisation
3934     Xini = numpy.zeros(Xb.shape)
3935     #
3936     # Définition de la fonction-coût
3937     # ------------------------------
3938     def CostFunction(v):
3939         _V = numpy.asmatrix(numpy.ravel( v )).T
3940         _X = Xb + B * _V
3941         if selfA._parameters["StoreInternalVariables"] or \
3942             selfA._toStore("CurrentState") or \
3943             selfA._toStore("CurrentOptimum"):
3944             selfA.StoredVariables["CurrentState"].store( _X )
3945         _HX = Hm( _X )
3946         _HX = numpy.asmatrix(numpy.ravel( _HX )).T
3947         _Innovation = Y - _HX
3948         if selfA._toStore("SimulatedObservationAtCurrentState") or \
3949             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3950             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
3951         if selfA._toStore("InnovationAtCurrentState"):
3952             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
3953         #
3954         Jb  = float( 0.5 * _V.T * BT * _V )
3955         Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
3956         J   = Jb + Jo
3957         #
3958         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
3959         selfA.StoredVariables["CostFunctionJb"].store( Jb )
3960         selfA.StoredVariables["CostFunctionJo"].store( Jo )
3961         selfA.StoredVariables["CostFunctionJ" ].store( J )
3962         if selfA._toStore("IndexOfOptimum") or \
3963             selfA._toStore("CurrentOptimum") or \
3964             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
3965             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
3966             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
3967             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3968             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3969         if selfA._toStore("IndexOfOptimum"):
3970             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3971         if selfA._toStore("CurrentOptimum"):
3972             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
3973         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3974             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
3975         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3976             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3977         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3978             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3979         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3980             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3981         return J
3982     #
3983     def GradientOfCostFunction(v):
3984         _V = numpy.asmatrix(numpy.ravel( v )).T
3985         _X = Xb + B * _V
3986         _HX     = Hm( _X )
3987         _HX     = numpy.asmatrix(numpy.ravel( _HX )).T
3988         GradJb  = BT * _V
3989         GradJo  = - Ha( (_X, RI * (Y - _HX)) )
3990         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
3991         return GradJ
3992     #
3993     # Minimisation de la fonctionnelle
3994     # --------------------------------
3995     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
3996     #
3997     if selfA._parameters["Minimizer"] == "LBFGSB":
3998         if "0.19" <= scipy.version.version <= "1.1.0":
3999             import lbfgsbhlt as optimiseur
4000         else:
4001             import scipy.optimize as optimiseur
4002         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
4003             func        = CostFunction,
4004             x0          = Xini,
4005             fprime      = GradientOfCostFunction,
4006             args        = (),
4007             bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
4008             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
4009             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
4010             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
4011             iprint      = selfA._parameters["optiprint"],
4012             )
4013         nfeval = Informations['funcalls']
4014         rc     = Informations['warnflag']
4015     elif selfA._parameters["Minimizer"] == "TNC":
4016         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
4017             func        = CostFunction,
4018             x0          = Xini,
4019             fprime      = GradientOfCostFunction,
4020             args        = (),
4021             bounds      = RecentredBounds(selfA._parameters["Bounds"], Xb),
4022             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
4023             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
4024             ftol        = selfA._parameters["CostDecrementTolerance"],
4025             messages    = selfA._parameters["optmessages"],
4026             )
4027     elif selfA._parameters["Minimizer"] == "CG":
4028         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
4029             f           = CostFunction,
4030             x0          = Xini,
4031             fprime      = GradientOfCostFunction,
4032             args        = (),
4033             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
4034             gtol        = selfA._parameters["GradientNormTolerance"],
4035             disp        = selfA._parameters["optdisp"],
4036             full_output = True,
4037             )
4038     elif selfA._parameters["Minimizer"] == "NCG":
4039         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
4040             f           = CostFunction,
4041             x0          = Xini,
4042             fprime      = GradientOfCostFunction,
4043             args        = (),
4044             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
4045             avextol     = selfA._parameters["CostDecrementTolerance"],
4046             disp        = selfA._parameters["optdisp"],
4047             full_output = True,
4048             )
4049     elif selfA._parameters["Minimizer"] == "BFGS":
4050         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
4051             f           = CostFunction,
4052             x0          = Xini,
4053             fprime      = GradientOfCostFunction,
4054             args        = (),
4055             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
4056             gtol        = selfA._parameters["GradientNormTolerance"],
4057             disp        = selfA._parameters["optdisp"],
4058             full_output = True,
4059             )
4060     else:
4061         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
4062     #
4063     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
4064     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
4065     #
4066     # Correction pour pallier a un bug de TNC sur le retour du Minimum
4067     # ----------------------------------------------------------------
4068     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
4069         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
4070         Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
4071     else:
4072         Minimum = Xb + B * numpy.asmatrix(numpy.ravel( Minimum )).T
4073     #
4074     # Obtention de l'analyse
4075     # ----------------------
4076     Xa = Minimum
4077     #
4078     selfA.StoredVariables["Analysis"].store( Xa )
4079     #
4080     if selfA._toStore("OMA") or \
4081         selfA._toStore("SigmaObs2") or \
4082         selfA._toStore("SimulationQuantiles") or \
4083         selfA._toStore("SimulatedObservationAtOptimum"):
4084         if selfA._toStore("SimulatedObservationAtCurrentState"):
4085             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
4086         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
4087             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
4088         else:
4089             HXa = Hm( Xa )
4090     #
4091     # Calcul de la covariance d'analyse
4092     # ---------------------------------
4093     if selfA._toStore("APosterioriCovariance") or \
4094         selfA._toStore("SimulationQuantiles") or \
4095         selfA._toStore("JacobianMatrixAtOptimum") or \
4096         selfA._toStore("KalmanGainAtOptimum"):
4097         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
4098         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
4099     if selfA._toStore("APosterioriCovariance") or \
4100         selfA._toStore("SimulationQuantiles") or \
4101         selfA._toStore("KalmanGainAtOptimum"):
4102         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
4103         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
4104     if selfA._toStore("APosterioriCovariance") or \
4105         selfA._toStore("SimulationQuantiles"):
4106         BI = B.getI()
4107         HessienneI = []
4108         nb = Xa.size
4109         for i in range(nb):
4110             _ee    = numpy.matrix(numpy.zeros(nb)).T
4111             _ee[i] = 1.
4112             _HtEE  = numpy.dot(HtM,_ee)
4113             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
4114             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
4115         HessienneI = numpy.matrix( HessienneI )
4116         A = HessienneI.I
4117         if min(A.shape) != max(A.shape):
4118             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)))
4119         if (numpy.diag(A) < 0).any():
4120             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,))
4121         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
4122             try:
4123                 L = numpy.linalg.cholesky( A )
4124             except:
4125                 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,))
4126     if selfA._toStore("APosterioriCovariance"):
4127         selfA.StoredVariables["APosterioriCovariance"].store( A )
4128     if selfA._toStore("JacobianMatrixAtOptimum"):
4129         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
4130     if selfA._toStore("KalmanGainAtOptimum"):
4131         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
4132         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
4133         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
4134     #
4135     # Calculs et/ou stockages supplémentaires
4136     # ---------------------------------------
4137     if selfA._toStore("Innovation") or \
4138         selfA._toStore("SigmaObs2") or \
4139         selfA._toStore("MahalanobisConsistency") or \
4140         selfA._toStore("OMB"):
4141         d  = Y - HXb
4142     if selfA._toStore("Innovation"):
4143         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
4144     if selfA._toStore("BMA"):
4145         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
4146     if selfA._toStore("OMA"):
4147         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
4148     if selfA._toStore("OMB"):
4149         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
4150     if selfA._toStore("SigmaObs2"):
4151         TraceR = R.trace(Y.size)
4152         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
4153     if selfA._toStore("MahalanobisConsistency"):
4154         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
4155     if selfA._toStore("SimulationQuantiles"):
4156         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
4157     if selfA._toStore("SimulatedObservationAtBackground"):
4158         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
4159     if selfA._toStore("SimulatedObservationAtOptimum"):
4160         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
4161     #
4162     return 0
4163
4164 # ==============================================================================
4165 if __name__ == "__main__":
4166     print('\n AUTODIAGNOSTIC\n')