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 EnsembleOfAnomalies( Ensemble, OptMean = None, Normalisation = 1.):
 508     "Renvoie les anomalies centrées à partir d'un ensemble TailleEtat*NbMembres"
 509     if OptMean is None:
 510         __Em = numpy.asarray(Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
 511     else:
 512         __Em = numpy.ravel(OptMean).reshape((-1,1))
 513     #
 514     return Normalisation * (numpy.asarray(Ensemble) - __Em)
 515
 516 # ==============================================================================
 517 def EnsembleErrorCovariance( Ensemble, __quick = False ):
 518     "Renvoie l'estimation empirique de la covariance d'ensemble"
 519     if __quick:
 520         # Covariance rapide mais rarement définie positive
 521         __Covariance = numpy.cov(Ensemble)
 522     else:
 523         # Résultat souvent identique à numpy.cov, mais plus robuste
 524         __n, __m = numpy.asarray(Ensemble).shape
 525         __Anomalies = EnsembleOfAnomalies( Ensemble )
 526         # Estimation empirique
 527         __Covariance = (__Anomalies @ __Anomalies.T) / (__m-1)
 528         # Assure la symétrie
 529         __Covariance = (__Covariance + __Covariance.T) * 0.5
 530         # Assure la positivité
 531         __epsilon    = mpr*numpy.trace(__Covariance)
 532         __Covariance = __Covariance + __epsilon * numpy.identity(__n)
 533     #
 534     return __Covariance
 535
 536 # ==============================================================================
 537 def EnsemblePerturbationWithGivenCovariance( __Ensemble, __Covariance, __Seed=None ):
 538     "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
 539     if hasattr(__Covariance,"assparsematrix"):
 540         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
 541             # Traitement d'une covariance nulle ou presque
 542             return __Ensemble
 543         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
 544             # Traitement d'une covariance nulle ou presque
 545             return __Ensemble
 546     else:
 547         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
 548             # Traitement d'une covariance nulle ou presque
 549             return __Ensemble
 550         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
 551             # Traitement d'une covariance nulle ou presque
 552             return __Ensemble
 553     #
 554     __n, __m = __Ensemble.shape
 555     if __Seed is not None: numpy.random.seed(__Seed)
 556     #
 557     if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
 558         # Traitement d'une covariance multiple de l'identité
 559         __zero = 0.
 560         __std  = numpy.sqrt(__Covariance.assparsematrix())
 561         __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
 562     #
 563     elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
 564         # Traitement d'une covariance diagonale avec variances non identiques
 565         __zero = numpy.zeros(__n)
 566         __std  = numpy.sqrt(__Covariance.assparsematrix())
 567         __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
 568     #
 569     elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
 570         # Traitement d'une covariance pleine
 571         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
 572     #
 573     elif isinstance(__Covariance, numpy.ndarray):
 574         # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
 575         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
 576     #
 577     else:
 578         raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
 579     #
 580     return __Ensemble
 581
 582 # ==============================================================================
 583 def CovarianceInflation(
 584         InputCovOrEns,
 585         InflationType   = None,
 586         InflationFactor = None,
 587         BackgroundCov   = None,
 588         ):
 589     """
 590     Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
 591
 592     Synthèse : Hunt 2007, section 2.3.5
 593     """
 594     if InflationFactor is None:
 595         return InputCovOrEns
 596     else:
 597         InflationFactor = float(InflationFactor)
 598     #
 599     if InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
 600         if InflationFactor < 1.:
 601             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 602         if InflationFactor < 1.+mpr:
 603             return InputCovOrEns
 604         OutputCovOrEns = InflationFactor**2 * InputCovOrEns
 605     #
 606     elif InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
 607         if InflationFactor < 1.:
 608             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 609         if InflationFactor < 1.+mpr:
 610             return InputCovOrEns
 611         InputCovOrEnsMean = InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
 612         OutputCovOrEns = InputCovOrEnsMean[:,numpy.newaxis] \
 613             + InflationFactor * (InputCovOrEns - InputCovOrEnsMean[:,numpy.newaxis])
 614     #
 615     elif InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
 616         if InflationFactor < 0.:
 617             raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
 618         if InflationFactor < mpr:
 619             return InputCovOrEns
 620         __n, __m = numpy.asarray(InputCovOrEns).shape
 621         if __n != __m:
 622             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 623         OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * numpy.identity(__n)
 624     #
 625     elif InflationType == "HybridOnBackgroundCovariance":
 626         if InflationFactor < 0.:
 627             raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
 628         if InflationFactor < mpr:
 629             return InputCovOrEns
 630         __n, __m = numpy.asarray(InputCovOrEns).shape
 631         if __n != __m:
 632             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 633         if BackgroundCov is None:
 634             raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
 635         if InputCovOrEns.shape != BackgroundCov.shape:
 636             raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
 637         OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * BackgroundCov
 638     #
 639     elif InflationType == "Relaxation":
 640         raise NotImplementedError("InflationType Relaxation")
 641     #
 642     else:
 643         raise ValueError("Error in inflation type, '%s' is not a valid keyword."%InflationType)
 644     #
 645     return OutputCovOrEns
 646
 647 # ==============================================================================
 648 def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
 649     "Estimation des quantiles a posteriori (selfA est modifié)"
 650     nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
 651     #
 652     # Traitement des bornes
 653     if "StateBoundsForQuantiles" in selfA._parameters:
 654         LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
 655     elif "Bounds" in selfA._parameters:
 656         LBounds = selfA._parameters["Bounds"]  # Défaut raisonnable
 657     else:
 658         LBounds = None
 659     if LBounds is not None:
 660         def NoneRemove(paire):
 661             bmin, bmax = paire
 662             if bmin is None: bmin = numpy.finfo('float').min
 663             if bmax is None: bmax = numpy.finfo('float').max
 664             return [bmin, bmax]
 665         LBounds = numpy.matrix( [NoneRemove(paire) for paire in 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 enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterFormula"):
 713     """
 714     EnKS
 715     """
 716     #
 717     # Opérateurs
 718     H = HO["Direct"].appliedControledFormTo
 719     #
 720     if selfA._parameters["EstimationOf"] == "State":
 721         M = EM["Direct"].appliedControledFormTo
 722     #
 723     if CM is not None and "Tangent" in CM and U is not None:
 724         Cm = CM["Tangent"].asMatrix(Xb)
 725     else:
 726         Cm = None
 727     #
 728     # Précalcul des inversions de B et R
 729     RIdemi = R.sqrtmI()
 730     #
 731     # Durée d'observation et tailles
 732     LagL = selfA._parameters["SmootherLagL"]
 733     if (not hasattr(Y,"store")) or (not hasattr(Y,"stepnumber")):
 734         raise ValueError("Fixed-lag smoother requires a series of observation")
 735     if Y.stepnumber() < LagL:
 736         raise ValueError("Fixed-lag smoother requires a series of observation greater then the lag L")
 737     duration = Y.stepnumber()
 738     __p = numpy.cumprod(Y.shape())[-1]
 739     __n = Xb.size
 740     __m = selfA._parameters["NumberOfMembers"]
 741     #
 742     if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
 743     else:                         Pn = B
 744     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
 745         selfA.StoredVariables["Analysis"].store( Xb )
 746         if selfA._toStore("APosterioriCovariance"):
 747             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
 748             covarianceXa = Pn
 749     #
 750     # Calcul direct initial (on privilégie la mémorisation au recalcul)
 751     __seed = numpy.random.get_state()
 752     selfB = copy.deepcopy(selfA)
 753     selfB._parameters["StoreSupplementaryCalculations"] = ["CurrentEnsembleState"]
 754     if VariantM == "EnKS16-KalmanFilterFormula":
 755         etkf(selfB, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM = "KalmanFilterFormula")
 756     else:
 757         raise ValueError("VariantM has to be chosen in the authorized methods list.")
 758     if LagL > 0:
 759         EL  = selfB.StoredVariables["CurrentEnsembleState"][LagL-1]
 760     else:
 761         EL = EnsembleOfBackgroundPerturbations( Xb, None, __m ) # Cf. etkf
 762     selfA._parameters["SetSeed"] = numpy.random.set_state(__seed)
 763     #
 764     for step in range(LagL,duration-1):
 765         #
 766         sEL = selfB.StoredVariables["CurrentEnsembleState"][step+1-LagL:step+1]
 767         sEL.append(None)
 768         #
 769         if hasattr(Y,"store"):
 770             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
 771         else:
 772             Ynpu = numpy.ravel( Y ).reshape((__p,1))
 773         #
 774         if U is not None:
 775             if hasattr(U,"store") and len(U)>1:
 776                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
 777             elif hasattr(U,"store") and len(U)==1:
 778                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
 779             else:
 780                 Un = numpy.asmatrix(numpy.ravel( U )).T
 781         else:
 782             Un = None
 783         #
 784         #--------------------------
 785         if VariantM == "EnKS16-KalmanFilterFormula":
 786             if selfA._parameters["EstimationOf"] == "State": # Forecast
 787                 EL = M( [(EL[:,i], Un) for i in range(__m)],
 788                     argsAsSerie = True,
 789                     returnSerieAsArrayMatrix = True )
 790                 EL = EnsemblePerturbationWithGivenCovariance( EL, Q )
 791                 EZ = H( [(EL[:,i], Un) for i in range(__m)],
 792                     argsAsSerie = True,
 793                     returnSerieAsArrayMatrix = True )
 794                 if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
 795                     Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
 796                     EZ = EZ + Cm * Un
 797             elif selfA._parameters["EstimationOf"] == "Parameters":
 798                 # --- > Par principe, M = Id, Q = 0
 799                 EZ = H( [(EL[:,i], Un) for i in range(__m)],
 800                     argsAsSerie = True,
 801                     returnSerieAsArrayMatrix = True )
 802             #
 803             vEm   = EL.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
 804             vZm   = EZ.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
 805             #
 806             mS    = RIdemi @ EnsembleOfAnomalies( EZ, vZm, 1./math.sqrt(__m-1) )
 807             delta = RIdemi @ ( Ynpu - vZm )
 808             mT    = numpy.linalg.inv( numpy.identity(__m) + mS.T @ mS )
 809             vw    = mT @ mS.T @ delta
 810             #
 811             Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
 812             mU    = numpy.identity(__m)
 813             wTU   = (vw.reshape((__m,1)) + math.sqrt(__m-1) * Tdemi @ mU)
 814             #
 815             EX    = EnsembleOfAnomalies( EL, vEm, 1./math.sqrt(__m-1) )
 816             EL    = vEm + EX @ wTU
 817             #
 818             sEL[LagL] = EL
 819             for irl in range(LagL): # Lissage des L précédentes analysis
 820                 vEm = sEL[irl].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
 821                 EX = EnsembleOfAnomalies( sEL[irl], vEm, 1./math.sqrt(__m-1) )
 822                 sEL[irl] = vEm + EX @ wTU
 823             #
 824             # Conservation de l'analyse retrospective d'ordre 0 avant rotation
 825             Xa = sEL[0].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
 826             if selfA._toStore("APosterioriCovariance"):
 827                 EXn = sEL[0]
 828             #
 829             for irl in range(LagL):
 830                 sEL[irl] = sEL[irl+1]
 831             sEL[LagL] = None
 832         #--------------------------
 833         else:
 834             raise ValueError("VariantM has to be chosen in the authorized methods list.")
 835         #
 836         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 837         # ---> avec analysis
 838         selfA.StoredVariables["Analysis"].store( Xa )
 839         if selfA._toStore("APosterioriCovariance"):
 840             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(EXn) )
 841     #
 842     # Stockage des dernières analyses incomplètement remises à jour
 843     for irl in range(LagL):
 844         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
 845         Xa = sEL[irl].mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
 846         selfA.StoredVariables["Analysis"].store( Xa )
 847     #
 848     return 0
 849
 850 # ==============================================================================
 851 def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"):
 852     """
 853     Ensemble-Transform EnKF
 854     """
 855     if selfA._parameters["EstimationOf"] == "Parameters":
 856         selfA._parameters["StoreInternalVariables"] = True
 857     #
 858     # Opérateurs
 859     # ----------
 860     H = HO["Direct"].appliedControledFormTo
 861     #
 862     if selfA._parameters["EstimationOf"] == "State":
 863         M = EM["Direct"].appliedControledFormTo
 864     #
 865     if CM is not None and "Tangent" in CM and U is not None:
 866         Cm = CM["Tangent"].asMatrix(Xb)
 867     else:
 868         Cm = None
 869     #
 870     # Nombre de pas identique au nombre de pas d'observations
 871     # -------------------------------------------------------
 872     if hasattr(Y,"stepnumber"):
 873         duration = Y.stepnumber()
 874         __p = numpy.cumprod(Y.shape())[-1]
 875     else:
 876         duration = 2
 877         __p = numpy.array(Y).size
 878     #
 879     # Précalcul des inversions de B et R
 880     # ----------------------------------
 881     if selfA._parameters["StoreInternalVariables"] \
 882         or selfA._toStore("CostFunctionJ") \
 883         or selfA._toStore("CostFunctionJb") \
 884         or selfA._toStore("CostFunctionJo") \
 885         or selfA._toStore("CurrentOptimum") \
 886         or selfA._toStore("APosterioriCovariance"):
 887         BI = B.getI()
 888         RI = R.getI()
 889     elif VariantM != "KalmanFilterFormula":
 890         RI = R.getI()
 891     if VariantM == "KalmanFilterFormula":
 892         RIdemi = R.sqrtmI()
 893     #
 894     # Initialisation
 895     # --------------
 896     __n = Xb.size
 897     __m = selfA._parameters["NumberOfMembers"]
 898     if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
 899     else:                         Pn = B
 900     Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
 901     #~ Xn = EnsembleOfBackgroundPerturbations( Xb, Pn, __m )
 902     #
 903     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
 904         selfA.StoredVariables["Analysis"].store( Xb )
 905         if selfA._toStore("APosterioriCovariance"):
 906             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
 907             covarianceXa = Pn
 908     #
 909     previousJMinimum = numpy.finfo(float).max
 910     #
 911     for step in range(duration-1):
 912         if hasattr(Y,"store"):
 913             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
 914         else:
 915             Ynpu = numpy.ravel( Y ).reshape((__p,1))
 916         #
 917         if U is not None:
 918             if hasattr(U,"store") and len(U)>1:
 919                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
 920             elif hasattr(U,"store") and len(U)==1:
 921                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
 922             else:
 923                 Un = numpy.asmatrix(numpy.ravel( U )).T
 924         else:
 925             Un = None
 926         #
 927         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
 928             Xn = CovarianceInflation( Xn,
 929                 selfA._parameters["InflationType"],
 930                 selfA._parameters["InflationFactor"],
 931                 )
 932         #
 933         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
 934             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
 935                 argsAsSerie = True,
 936                 returnSerieAsArrayMatrix = True )
 937             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
 938             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
 939                 argsAsSerie = True,
 940                 returnSerieAsArrayMatrix = True )
 941             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
 942                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
 943                 Xn_predicted = Xn_predicted + Cm * Un
 944         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
 945             # --- > Par principe, M = Id, Q = 0
 946             Xn_predicted = Xn
 947             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
 948                 argsAsSerie = True,
 949                 returnSerieAsArrayMatrix = True )
 950         #
 951         # Mean of forecast and observation of forecast
 952         Xfm  = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
 953         Hfm  = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
 954         #
 955         # Anomalies
 956         EaX   = EnsembleOfAnomalies( Xn_predicted, Xfm )
 957         EaHX  = EnsembleOfAnomalies( HX_predicted, Hfm)
 958         #
 959         #--------------------------
 960         if VariantM == "KalmanFilterFormula":
 961             mS    = RIdemi * EaHX / math.sqrt(__m-1)
 962             delta = RIdemi * ( Ynpu - Hfm )
 963             mT    = numpy.linalg.inv( numpy.identity(__m) + mS.T @ mS )
 964             vw    = mT @ mS.T @ delta
 965             #
 966             Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
 967             mU    = numpy.identity(__m)
 968             #
 969             EaX   = EaX / math.sqrt(__m-1)
 970             Xn    = Xfm + EaX @ ( vw.reshape((__m,1)) + math.sqrt(__m-1) * Tdemi @ mU )
 971         #--------------------------
 972         elif VariantM == "Variational":
 973             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
 974             def CostFunction(w):
 975                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
 976                 _Jo = 0.5 * _A.T @ (RI * _A)
 977                 _Jb = 0.5 * (__m-1) * w.T @ w
 978                 _J  = _Jo + _Jb
 979                 return float(_J)
 980             def GradientOfCostFunction(w):
 981                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
 982                 _GardJo = - EaHX.T @ (RI * _A)
 983                 _GradJb = (__m-1) * w.reshape((__m,1))
 984                 _GradJ  = _GardJo + _GradJb
 985                 return numpy.ravel(_GradJ)
 986             vw = scipy.optimize.fmin_cg(
 987                 f           = CostFunction,
 988                 x0          = numpy.zeros(__m),
 989                 fprime      = GradientOfCostFunction,
 990                 args        = (),
 991                 disp        = False,
 992                 )
 993             #
 994             Hto = EaHX.T @ (RI * EaHX)
 995             Htb = (__m-1) * numpy.identity(__m)
 996             Hta = Hto + Htb
 997             #
 998             Pta = numpy.linalg.inv( Hta )
 999             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1000             #
1001             Xn  = Xfm + EaX @ (vw[:,None] + EWa)
1002         #--------------------------
1003         elif VariantM == "FiniteSize11": # Jauge Boc2011
1004             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1005             def CostFunction(w):
1006                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1007                 _Jo = 0.5 * _A.T @ (RI * _A)
1008                 _Jb = 0.5 * __m * math.log(1 + 1/__m + w.T @ w)
1009                 _J  = _Jo + _Jb
1010                 return float(_J)
1011             def GradientOfCostFunction(w):
1012                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1013                 _GardJo = - EaHX.T @ (RI * _A)
1014                 _GradJb = __m * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
1015                 _GradJ  = _GardJo + _GradJb
1016                 return numpy.ravel(_GradJ)
1017             vw = scipy.optimize.fmin_cg(
1018                 f           = CostFunction,
1019                 x0          = numpy.zeros(__m),
1020                 fprime      = GradientOfCostFunction,
1021                 args        = (),
1022                 disp        = False,
1023                 )
1024             #
1025             Hto = EaHX.T @ (RI * EaHX)
1026             Htb = __m * \
1027                 ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
1028                 / (1 + 1/__m + vw.T @ vw)**2
1029             Hta = Hto + Htb
1030             #
1031             Pta = numpy.linalg.inv( Hta )
1032             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1033             #
1034             Xn  = Xfm + EaX @ (vw.reshape((__m,1)) + EWa)
1035         #--------------------------
1036         elif VariantM == "FiniteSize15": # Jauge Boc2015
1037             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1038             def CostFunction(w):
1039                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1040                 _Jo = 0.5 * _A.T * RI * _A
1041                 _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w)
1042                 _J  = _Jo + _Jb
1043                 return float(_J)
1044             def GradientOfCostFunction(w):
1045                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1046                 _GardJo = - EaHX.T @ (RI * _A)
1047                 _GradJb = (__m+1) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
1048                 _GradJ  = _GardJo + _GradJb
1049                 return numpy.ravel(_GradJ)
1050             vw = scipy.optimize.fmin_cg(
1051                 f           = CostFunction,
1052                 x0          = numpy.zeros(__m),
1053                 fprime      = GradientOfCostFunction,
1054                 args        = (),
1055                 disp        = False,
1056                 )
1057             #
1058             Hto = EaHX.T @ (RI * EaHX)
1059             Htb = (__m+1) * \
1060                 ( (1 + 1/__m + vw.T @ vw) * numpy.identity(__m) - 2 * vw @ vw.T ) \
1061                 / (1 + 1/__m + vw.T @ vw)**2
1062             Hta = Hto + Htb
1063             #
1064             Pta = numpy.linalg.inv( Hta )
1065             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1066             #
1067             Xn  = Xfm + EaX @ (vw.reshape((__m,1)) + EWa)
1068         #--------------------------
1069         elif VariantM == "FiniteSize16": # Jauge Boc2016
1070             HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm
1071             def CostFunction(w):
1072                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1073                 _Jo = 0.5 * _A.T @ (RI * _A)
1074                 _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w / (__m-1))
1075                 _J  = _Jo + _Jb
1076                 return float(_J)
1077             def GradientOfCostFunction(w):
1078                 _A  = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1))
1079                 _GardJo = - EaHX.T @ (RI * _A)
1080                 _GradJb = ((__m+1) / (__m-1)) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w / (__m-1))
1081                 _GradJ  = _GardJo + _GradJb
1082                 return numpy.ravel(_GradJ)
1083             vw = scipy.optimize.fmin_cg(
1084                 f           = CostFunction,
1085                 x0          = numpy.zeros(__m),
1086                 fprime      = GradientOfCostFunction,
1087                 args        = (),
1088                 disp        = False,
1089                 )
1090             #
1091             Hto = EaHX.T @ (RI * EaHX)
1092             Htb = ((__m+1) / (__m-1)) * \
1093                 ( (1 + 1/__m + vw.T @ vw / (__m-1)) * numpy.identity(__m) - 2 * vw @ vw.T / (__m-1) ) \
1094                 / (1 + 1/__m + vw.T @ vw / (__m-1))**2
1095             Hta = Hto + Htb
1096             #
1097             Pta = numpy.linalg.inv( Hta )
1098             EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
1099             #
1100             Xn  = Xfm + EaX @ (vw[:,None] + EWa)
1101         #--------------------------
1102         else:
1103             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1104         #
1105         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
1106             Xn = CovarianceInflation( Xn,
1107                 selfA._parameters["InflationType"],
1108                 selfA._parameters["InflationFactor"],
1109                 )
1110         #
1111         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1112         #--------------------------
1113         #
1114         if selfA._parameters["StoreInternalVariables"] \
1115             or selfA._toStore("CostFunctionJ") \
1116             or selfA._toStore("CostFunctionJb") \
1117             or selfA._toStore("CostFunctionJo") \
1118             or selfA._toStore("APosterioriCovariance") \
1119             or selfA._toStore("InnovationAtCurrentAnalysis") \
1120             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
1121             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1122             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
1123             _Innovation = Ynpu - _HXa
1124         #
1125         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1126         # ---> avec analysis
1127         selfA.StoredVariables["Analysis"].store( Xa )
1128         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1129             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
1130         if selfA._toStore("InnovationAtCurrentAnalysis"):
1131             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1132         # ---> avec current state
1133         if selfA._parameters["StoreInternalVariables"] \
1134             or selfA._toStore("CurrentState"):
1135             selfA.StoredVariables["CurrentState"].store( Xn )
1136         if selfA._toStore("ForecastState"):
1137             selfA.StoredVariables["ForecastState"].store( EMX )
1138         if selfA._toStore("BMA"):
1139             selfA.StoredVariables["BMA"].store( EMX - Xa.reshape((__n,1)) )
1140         if selfA._toStore("InnovationAtCurrentState"):
1141             selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
1142         if selfA._toStore("SimulatedObservationAtCurrentState") \
1143             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1144             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
1145         # ---> autres
1146         if selfA._parameters["StoreInternalVariables"] \
1147             or selfA._toStore("CostFunctionJ") \
1148             or selfA._toStore("CostFunctionJb") \
1149             or selfA._toStore("CostFunctionJo") \
1150             or selfA._toStore("CurrentOptimum") \
1151             or selfA._toStore("APosterioriCovariance"):
1152             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
1153             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
1154             J   = Jb + Jo
1155             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1156             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1157             selfA.StoredVariables["CostFunctionJ" ].store( J )
1158             #
1159             if selfA._toStore("IndexOfOptimum") \
1160                 or selfA._toStore("CurrentOptimum") \
1161                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1162                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1163                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1164                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1165                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1166             if selfA._toStore("IndexOfOptimum"):
1167                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1168             if selfA._toStore("CurrentOptimum"):
1169                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1170             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1171                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1172             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1173                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1174             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1175                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1176             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1177                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1178         if selfA._toStore("APosterioriCovariance"):
1179             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
1180         if selfA._parameters["EstimationOf"] == "Parameters" \
1181             and J < previousJMinimum:
1182             previousJMinimum    = J
1183             XaMin               = Xa
1184             if selfA._toStore("APosterioriCovariance"):
1185                 covarianceXaMin = Pn
1186         # ---> Pour les smoothers
1187         if selfA._toStore("CurrentEnsembleState"):
1188             selfA.StoredVariables["CurrentEnsembleState"].store( Xn )
1189     #
1190     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1191     # ----------------------------------------------------------------------
1192     if selfA._parameters["EstimationOf"] == "Parameters":
1193         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1194         selfA.StoredVariables["Analysis"].store( XaMin )
1195         if selfA._toStore("APosterioriCovariance"):
1196             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1197         if selfA._toStore("BMA"):
1198             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1199     #
1200     return 0
1201
1202 # ==============================================================================
1203 def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12",
1204     BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000):
1205     """
1206     Iterative EnKF
1207     """
1208     if selfA._parameters["EstimationOf"] == "Parameters":
1209         selfA._parameters["StoreInternalVariables"] = True
1210     #
1211     # Opérateurs
1212     # ----------
1213     H = HO["Direct"].appliedControledFormTo
1214     #
1215     if selfA._parameters["EstimationOf"] == "State":
1216         M = EM["Direct"].appliedControledFormTo
1217     #
1218     if CM is not None and "Tangent" in CM and U is not None:
1219         Cm = CM["Tangent"].asMatrix(Xb)
1220     else:
1221         Cm = None
1222     #
1223     # Nombre de pas identique au nombre de pas d'observations
1224     # -------------------------------------------------------
1225     if hasattr(Y,"stepnumber"):
1226         duration = Y.stepnumber()
1227         __p = numpy.cumprod(Y.shape())[-1]
1228     else:
1229         duration = 2
1230         __p = numpy.array(Y).size
1231     #
1232     # Précalcul des inversions de B et R
1233     # ----------------------------------
1234     if selfA._parameters["StoreInternalVariables"] \
1235         or selfA._toStore("CostFunctionJ") \
1236         or selfA._toStore("CostFunctionJb") \
1237         or selfA._toStore("CostFunctionJo") \
1238         or selfA._toStore("CurrentOptimum") \
1239         or selfA._toStore("APosterioriCovariance"):
1240         BI = B.getI()
1241     RI = R.getI()
1242     #
1243     # Initialisation
1244     # --------------
1245     __n = Xb.size
1246     __m = selfA._parameters["NumberOfMembers"]
1247     if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
1248     else:                         Pn = B
1249     if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
1250     else:                         Rn = R
1251     if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
1252     else:                         Qn = Q
1253     Xn = EnsembleOfBackgroundPerturbations( Xb, Pn, __m )
1254     #
1255     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
1256         selfA.StoredVariables["Analysis"].store( Xb )
1257         if selfA._toStore("APosterioriCovariance"):
1258             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
1259             covarianceXa = Pn
1260     #
1261     previousJMinimum = numpy.finfo(float).max
1262     #
1263     for step in range(duration-1):
1264         if hasattr(Y,"store"):
1265             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
1266         else:
1267             Ynpu = numpy.ravel( Y ).reshape((__p,1))
1268         #
1269         if U is not None:
1270             if hasattr(U,"store") and len(U)>1:
1271                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
1272             elif hasattr(U,"store") and len(U)==1:
1273                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
1274             else:
1275                 Un = numpy.asmatrix(numpy.ravel( U )).T
1276         else:
1277             Un = None
1278         #
1279         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
1280             Xn = CovarianceInflation( Xn,
1281                 selfA._parameters["InflationType"],
1282                 selfA._parameters["InflationFactor"],
1283                 )
1284         #
1285         #--------------------------
1286         if VariantM == "IEnKF12":
1287             Xfm = numpy.ravel(Xn.mean(axis=1, dtype=mfp).astype('float'))
1288             EaX = EnsembleOfAnomalies( Xn ) / math.sqrt(__m-1)
1289             __j = 0
1290             Deltaw = 1
1291             if not BnotT:
1292                 Ta  = numpy.identity(__m)
1293             vw  = numpy.zeros(__m)
1294             while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
1295                 vx1 = (Xfm + EaX @ vw).reshape((__n,1))
1296                 #
1297                 if BnotT:
1298                     E1 = vx1 + _epsilon * EaX
1299                 else:
1300                     E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
1301                 #
1302                 if selfA._parameters["EstimationOf"] == "State": # Forecast + Q
1303                     E2 = M( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
1304                         argsAsSerie = True,
1305                         returnSerieAsArrayMatrix = True )
1306                 elif selfA._parameters["EstimationOf"] == "Parameters":
1307                     # --- > Par principe, M = Id
1308                     E2 = Xn
1309                 vx2 = E2.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1310                 vy1 = H((vx2, Un)).reshape((__p,1))
1311                 #
1312                 HE2 = H( [(E2[:,i,numpy.newaxis], Un) for i in range(__m)],
1313                     argsAsSerie = True,
1314                     returnSerieAsArrayMatrix = True )
1315                 vy2 = HE2.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
1316                 #
1317                 if BnotT:
1318                     EaY = (HE2 - vy2) / _epsilon
1319                 else:
1320                     EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
1321                 #
1322                 GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy1 )))
1323                 mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY)
1324                 Deltaw = - numpy.linalg.solve(mH,GradJ)
1325                 #
1326                 vw = vw + Deltaw
1327                 #
1328                 if not BnotT:
1329                     Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1330                 #
1331                 __j = __j + 1
1332             #
1333             A2 = EnsembleOfAnomalies( E2 )
1334             #
1335             if BnotT:
1336                 Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1337                 A2 = math.sqrt(__m-1) * A2 @ Ta / _epsilon
1338             #
1339             Xn = vx2 + A2
1340         #--------------------------
1341         else:
1342             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1343         #
1344         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
1345             Xn = CovarianceInflation( Xn,
1346                 selfA._parameters["InflationType"],
1347                 selfA._parameters["InflationFactor"],
1348                 )
1349         #
1350         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1351         #--------------------------
1352         #
1353         if selfA._parameters["StoreInternalVariables"] \
1354             or selfA._toStore("CostFunctionJ") \
1355             or selfA._toStore("CostFunctionJb") \
1356             or selfA._toStore("CostFunctionJo") \
1357             or selfA._toStore("APosterioriCovariance") \
1358             or selfA._toStore("InnovationAtCurrentAnalysis") \
1359             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
1360             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1361             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
1362             _Innovation = Ynpu - _HXa
1363         #
1364         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1365         # ---> avec analysis
1366         selfA.StoredVariables["Analysis"].store( Xa )
1367         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1368             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
1369         if selfA._toStore("InnovationAtCurrentAnalysis"):
1370             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1371         # ---> avec current state
1372         if selfA._parameters["StoreInternalVariables"] \
1373             or selfA._toStore("CurrentState"):
1374             selfA.StoredVariables["CurrentState"].store( Xn )
1375         if selfA._toStore("ForecastState"):
1376             selfA.StoredVariables["ForecastState"].store( E2 )
1377         if selfA._toStore("BMA"):
1378             selfA.StoredVariables["BMA"].store( E2 - Xa )
1379         if selfA._toStore("InnovationAtCurrentState"):
1380             selfA.StoredVariables["InnovationAtCurrentState"].store( - HE2 + Ynpu )
1381         if selfA._toStore("SimulatedObservationAtCurrentState") \
1382             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1383             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HE2 )
1384         # ---> autres
1385         if selfA._parameters["StoreInternalVariables"] \
1386             or selfA._toStore("CostFunctionJ") \
1387             or selfA._toStore("CostFunctionJb") \
1388             or selfA._toStore("CostFunctionJo") \
1389             or selfA._toStore("CurrentOptimum") \
1390             or selfA._toStore("APosterioriCovariance"):
1391             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
1392             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
1393             J   = Jb + Jo
1394             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1395             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1396             selfA.StoredVariables["CostFunctionJ" ].store( J )
1397             #
1398             if selfA._toStore("IndexOfOptimum") \
1399                 or selfA._toStore("CurrentOptimum") \
1400                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1401                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1402                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1403                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1404                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1405             if selfA._toStore("IndexOfOptimum"):
1406                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1407             if selfA._toStore("CurrentOptimum"):
1408                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1409             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1410                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1411             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1412                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1413             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1414                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1415             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1416                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1417         if selfA._toStore("APosterioriCovariance"):
1418             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
1419         if selfA._parameters["EstimationOf"] == "Parameters" \
1420             and J < previousJMinimum:
1421             previousJMinimum    = J
1422             XaMin               = Xa
1423             if selfA._toStore("APosterioriCovariance"):
1424                 covarianceXaMin = Pn
1425     #
1426     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1427     # ----------------------------------------------------------------------
1428     if selfA._parameters["EstimationOf"] == "Parameters":
1429         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1430         selfA.StoredVariables["Analysis"].store( XaMin )
1431         if selfA._toStore("APosterioriCovariance"):
1432             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1433         if selfA._toStore("BMA"):
1434             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1435     #
1436     return 0
1437
1438 # ==============================================================================
1439 def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
1440     """
1441     3DVAR incrémental
1442     """
1443     #
1444     # Initialisations
1445     # ---------------
1446     #
1447     # Opérateur non-linéaire pour la boucle externe
1448     Hm = HO["Direct"].appliedTo
1449     #
1450     # Précalcul des inversions de B et R
1451     BI = B.getI()
1452     RI = R.getI()
1453     #
1454     # Point de démarrage de l'optimisation
1455     Xini = selfA._parameters["InitializationPoint"]
1456     #
1457     HXb = numpy.asmatrix(numpy.ravel( Hm( Xb ) )).T
1458     Innovation = Y - HXb
1459     #
1460     # Outer Loop
1461     # ----------
1462     iOuter = 0
1463     J      = 1./mpr
1464     DeltaJ = 1./mpr
1465     Xr     = Xini.reshape((-1,1))
1466     while abs(DeltaJ) >= selfA._parameters["CostDecrementTolerance"] and iOuter <= selfA._parameters["MaximumNumberOfSteps"]:
1467         #
1468         # Inner Loop
1469         # ----------
1470         Ht = HO["Tangent"].asMatrix(Xr)
1471         Ht = Ht.reshape(Y.size,Xr.size) # ADAO & check shape
1472         #
1473         # Définition de la fonction-coût
1474         # ------------------------------
1475         def CostFunction(dx):
1476             _dX  = numpy.asmatrix(numpy.ravel( dx )).T
1477             if selfA._parameters["StoreInternalVariables"] or \
1478                 selfA._toStore("CurrentState") or \
1479                 selfA._toStore("CurrentOptimum"):
1480                 selfA.StoredVariables["CurrentState"].store( Xb + _dX )
1481             _HdX = Ht * _dX
1482             _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T
1483             _dInnovation = Innovation - _HdX
1484             if selfA._toStore("SimulatedObservationAtCurrentState") or \
1485                 selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1486                 selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HXb + _HdX )
1487             if selfA._toStore("InnovationAtCurrentState"):
1488                 selfA.StoredVariables["InnovationAtCurrentState"].store( _dInnovation )
1489             #
1490             Jb  = float( 0.5 * _dX.T * BI * _dX )
1491             Jo  = float( 0.5 * _dInnovation.T * RI * _dInnovation )
1492             J   = Jb + Jo
1493             #
1494             selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
1495             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1496             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1497             selfA.StoredVariables["CostFunctionJ" ].store( J )
1498             if selfA._toStore("IndexOfOptimum") or \
1499                 selfA._toStore("CurrentOptimum") or \
1500                 selfA._toStore("CostFunctionJAtCurrentOptimum") or \
1501                 selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
1502                 selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
1503                 selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1504                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1505             if selfA._toStore("IndexOfOptimum"):
1506                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1507             if selfA._toStore("CurrentOptimum"):
1508                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
1509             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1510                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
1511             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1512                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1513             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1514                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1515             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1516                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1517             return J
1518         #
1519         def GradientOfCostFunction(dx):
1520             _dX          = numpy.asmatrix(numpy.ravel( dx )).T
1521             _HdX         = Ht * _dX
1522             _HdX         = numpy.asmatrix(numpy.ravel( _HdX )).T
1523             _dInnovation = Innovation - _HdX
1524             GradJb       = BI * _dX
1525             GradJo       = - Ht.T @ (RI * _dInnovation)
1526             GradJ        = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
1527             return GradJ
1528         #
1529         # Minimisation de la fonctionnelle
1530         # --------------------------------
1531         nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
1532         #
1533         if selfA._parameters["Minimizer"] == "LBFGSB":
1534             # Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
1535             if "0.19" <= scipy.version.version <= "1.1.0":
1536                 import lbfgsbhlt as optimiseur
1537             else:
1538                 import scipy.optimize as optimiseur
1539             Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
1540                 func        = CostFunction,
1541                 x0          = numpy.zeros(Xini.size),
1542                 fprime      = GradientOfCostFunction,
1543                 args        = (),
1544                 bounds      = selfA._parameters["Bounds"],
1545                 maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
1546                 factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
1547                 pgtol       = selfA._parameters["ProjectedGradientTolerance"],
1548                 iprint      = selfA._parameters["optiprint"],
1549                 )
1550             nfeval = Informations['funcalls']
1551             rc     = Informations['warnflag']
1552         elif selfA._parameters["Minimizer"] == "TNC":
1553             Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
1554                 func        = CostFunction,
1555                 x0          = numpy.zeros(Xini.size),
1556                 fprime      = GradientOfCostFunction,
1557                 args        = (),
1558                 bounds      = selfA._parameters["Bounds"],
1559                 maxfun      = selfA._parameters["MaximumNumberOfSteps"],
1560                 pgtol       = selfA._parameters["ProjectedGradientTolerance"],
1561                 ftol        = selfA._parameters["CostDecrementTolerance"],
1562                 messages    = selfA._parameters["optmessages"],
1563                 )
1564         elif selfA._parameters["Minimizer"] == "CG":
1565             Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
1566                 f           = CostFunction,
1567                 x0          = numpy.zeros(Xini.size),
1568                 fprime      = GradientOfCostFunction,
1569                 args        = (),
1570                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
1571                 gtol        = selfA._parameters["GradientNormTolerance"],
1572                 disp        = selfA._parameters["optdisp"],
1573                 full_output = True,
1574                 )
1575         elif selfA._parameters["Minimizer"] == "NCG":
1576             Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
1577                 f           = CostFunction,
1578                 x0          = numpy.zeros(Xini.size),
1579                 fprime      = GradientOfCostFunction,
1580                 args        = (),
1581                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
1582                 avextol     = selfA._parameters["CostDecrementTolerance"],
1583                 disp        = selfA._parameters["optdisp"],
1584                 full_output = True,
1585                 )
1586         elif selfA._parameters["Minimizer"] == "BFGS":
1587             Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
1588                 f           = CostFunction,
1589                 x0          = numpy.zeros(Xini.size),
1590                 fprime      = GradientOfCostFunction,
1591                 args        = (),
1592                 maxiter     = selfA._parameters["MaximumNumberOfSteps"],
1593                 gtol        = selfA._parameters["GradientNormTolerance"],
1594                 disp        = selfA._parameters["optdisp"],
1595                 full_output = True,
1596                 )
1597         else:
1598             raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
1599         #
1600         IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1601         MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
1602         #
1603         if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
1604             Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
1605             Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
1606         else:
1607             Minimum = Xb + numpy.asmatrix(numpy.ravel( Minimum )).T
1608         #
1609         Xr     = Minimum
1610         DeltaJ = selfA.StoredVariables["CostFunctionJ" ][-1] - J
1611         iOuter = selfA.StoredVariables["CurrentIterationNumber"][-1]
1612     #
1613     # Obtention de l'analyse
1614     # ----------------------
1615     Xa = Xr
1616     #
1617     selfA.StoredVariables["Analysis"].store( Xa )
1618     #
1619     if selfA._toStore("OMA") or \
1620         selfA._toStore("SigmaObs2") or \
1621         selfA._toStore("SimulationQuantiles") or \
1622         selfA._toStore("SimulatedObservationAtOptimum"):
1623         if selfA._toStore("SimulatedObservationAtCurrentState"):
1624             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
1625         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1626             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
1627         else:
1628             HXa = Hm( Xa )
1629     #
1630     # Calcul de la covariance d'analyse
1631     # ---------------------------------
1632     if selfA._toStore("APosterioriCovariance") or \
1633         selfA._toStore("SimulationQuantiles") or \
1634         selfA._toStore("JacobianMatrixAtOptimum") or \
1635         selfA._toStore("KalmanGainAtOptimum"):
1636         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
1637         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
1638     if selfA._toStore("APosterioriCovariance") or \
1639         selfA._toStore("SimulationQuantiles") or \
1640         selfA._toStore("KalmanGainAtOptimum"):
1641         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
1642         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
1643     if selfA._toStore("APosterioriCovariance") or \
1644         selfA._toStore("SimulationQuantiles"):
1645         HessienneI = []
1646         nb = Xa.size
1647         for i in range(nb):
1648             _ee    = numpy.matrix(numpy.zeros(nb)).T
1649             _ee[i] = 1.
1650             _HtEE  = numpy.dot(HtM,_ee)
1651             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
1652             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
1653         HessienneI = numpy.matrix( HessienneI )
1654         A = HessienneI.I
1655         if min(A.shape) != max(A.shape):
1656             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)))
1657         if (numpy.diag(A) < 0).any():
1658             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,))
1659         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
1660             try:
1661                 L = numpy.linalg.cholesky( A )
1662             except:
1663                 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,))
1664     if selfA._toStore("APosterioriCovariance"):
1665         selfA.StoredVariables["APosterioriCovariance"].store( A )
1666     if selfA._toStore("JacobianMatrixAtOptimum"):
1667         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
1668     if selfA._toStore("KalmanGainAtOptimum"):
1669         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
1670         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
1671         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
1672     #
1673     # Calculs et/ou stockages supplémentaires
1674     # ---------------------------------------
1675     if selfA._toStore("Innovation") or \
1676         selfA._toStore("SigmaObs2") or \
1677         selfA._toStore("MahalanobisConsistency") or \
1678         selfA._toStore("OMB"):
1679         d  = Y - HXb
1680     if selfA._toStore("Innovation"):
1681         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
1682     if selfA._toStore("BMA"):
1683         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
1684     if selfA._toStore("OMA"):
1685         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
1686     if selfA._toStore("OMB"):
1687         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
1688     if selfA._toStore("SigmaObs2"):
1689         TraceR = R.trace(Y.size)
1690         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
1691     if selfA._toStore("MahalanobisConsistency"):
1692         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
1693     if selfA._toStore("SimulationQuantiles"):
1694         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
1695     if selfA._toStore("SimulatedObservationAtBackground"):
1696         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
1697     if selfA._toStore("SimulatedObservationAtOptimum"):
1698         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
1699     #
1700     return 0
1701
1702 # ==============================================================================
1703 def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13",
1704     BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000):
1705     """
1706     Maximum Likelihood Ensemble Filter
1707     """
1708     if selfA._parameters["EstimationOf"] == "Parameters":
1709         selfA._parameters["StoreInternalVariables"] = True
1710     #
1711     # Opérateurs
1712     # ----------
1713     H = HO["Direct"].appliedControledFormTo
1714     #
1715     if selfA._parameters["EstimationOf"] == "State":
1716         M = EM["Direct"].appliedControledFormTo
1717     #
1718     if CM is not None and "Tangent" in CM and U is not None:
1719         Cm = CM["Tangent"].asMatrix(Xb)
1720     else:
1721         Cm = None
1722     #
1723     # Nombre de pas identique au nombre de pas d'observations
1724     # -------------------------------------------------------
1725     if hasattr(Y,"stepnumber"):
1726         duration = Y.stepnumber()
1727         __p = numpy.cumprod(Y.shape())[-1]
1728     else:
1729         duration = 2
1730         __p = numpy.array(Y).size
1731     #
1732     # Précalcul des inversions de B et R
1733     # ----------------------------------
1734     if selfA._parameters["StoreInternalVariables"] \
1735         or selfA._toStore("CostFunctionJ") \
1736         or selfA._toStore("CostFunctionJb") \
1737         or selfA._toStore("CostFunctionJo") \
1738         or selfA._toStore("CurrentOptimum") \
1739         or selfA._toStore("APosterioriCovariance"):
1740         BI = B.getI()
1741     RI = R.getI()
1742     #
1743     # Initialisation
1744     # --------------
1745     __n = Xb.size
1746     __m = selfA._parameters["NumberOfMembers"]
1747     if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
1748     else:                         Pn = B
1749     if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
1750     else:                         Rn = R
1751     Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
1752     #
1753     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
1754         selfA.StoredVariables["Analysis"].store( Xb )
1755         if selfA._toStore("APosterioriCovariance"):
1756             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
1757             covarianceXa = Pn
1758     #
1759     previousJMinimum = numpy.finfo(float).max
1760     #
1761     for step in range(duration-1):
1762         if hasattr(Y,"store"):
1763             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
1764         else:
1765             Ynpu = numpy.ravel( Y ).reshape((__p,1))
1766         #
1767         if U is not None:
1768             if hasattr(U,"store") and len(U)>1:
1769                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
1770             elif hasattr(U,"store") and len(U)==1:
1771                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
1772             else:
1773                 Un = numpy.asmatrix(numpy.ravel( U )).T
1774         else:
1775             Un = None
1776         #
1777         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
1778             Xn = CovarianceInflation( Xn,
1779                 selfA._parameters["InflationType"],
1780                 selfA._parameters["InflationFactor"],
1781                 )
1782         #
1783         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
1784             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
1785                 argsAsSerie = True,
1786                 returnSerieAsArrayMatrix = True )
1787             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
1788             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
1789                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
1790                 Xn_predicted = Xn_predicted + Cm * Un
1791         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
1792             # --- > Par principe, M = Id, Q = 0
1793             Xn_predicted = Xn
1794         #
1795         #--------------------------
1796         if VariantM == "MLEF13":
1797             Xfm = numpy.ravel(Xn_predicted.mean(axis=1, dtype=mfp).astype('float'))
1798             EaX = EnsembleOfAnomalies( Xn_predicted, Xfm, 1./math.sqrt(__m-1) )
1799             Ua  = numpy.identity(__m)
1800             __j = 0
1801             Deltaw = 1
1802             if not BnotT:
1803                 Ta  = numpy.identity(__m)
1804             vw  = numpy.zeros(__m)
1805             while numpy.linalg.norm(Deltaw) >= _e and __j <= _jmax:
1806                 vx1 = (Xfm + EaX @ vw).reshape((__n,1))
1807                 #
1808                 if BnotT:
1809                     E1 = vx1 + _epsilon * EaX
1810                 else:
1811                     E1 = vx1 + math.sqrt(__m-1) * EaX @ Ta
1812                 #
1813                 HE2 = H( [(E1[:,i,numpy.newaxis], Un) for i in range(__m)],
1814                     argsAsSerie = True,
1815                     returnSerieAsArrayMatrix = True )
1816                 vy2 = HE2.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
1817                 #
1818                 if BnotT:
1819                     EaY = (HE2 - vy2) / _epsilon
1820                 else:
1821                     EaY = ( (HE2 - vy2) @ numpy.linalg.inv(Ta) ) / math.sqrt(__m-1)
1822                 #
1823                 GradJ = numpy.ravel(vw[:,None] - EaY.transpose() @ (RI * ( Ynpu - vy2 )))
1824                 mH = numpy.identity(__m) + EaY.transpose() @ (RI * EaY)
1825                 Deltaw = - numpy.linalg.solve(mH,GradJ)
1826                 #
1827                 vw = vw + Deltaw
1828                 #
1829                 if not BnotT:
1830                     Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1831                 #
1832                 __j = __j + 1
1833             #
1834             if BnotT:
1835                 Ta = numpy.real(scipy.linalg.sqrtm(numpy.linalg.inv( mH )))
1836             #
1837             Xn = vx1 + math.sqrt(__m-1) * EaX @ Ta @ Ua
1838         #--------------------------
1839         else:
1840             raise ValueError("VariantM has to be chosen in the authorized methods list.")
1841         #
1842         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
1843             Xn = CovarianceInflation( Xn,
1844                 selfA._parameters["InflationType"],
1845                 selfA._parameters["InflationFactor"],
1846                 )
1847         #
1848         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
1849         #--------------------------
1850         #
1851         if selfA._parameters["StoreInternalVariables"] \
1852             or selfA._toStore("CostFunctionJ") \
1853             or selfA._toStore("CostFunctionJb") \
1854             or selfA._toStore("CostFunctionJo") \
1855             or selfA._toStore("APosterioriCovariance") \
1856             or selfA._toStore("InnovationAtCurrentAnalysis") \
1857             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
1858             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1859             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
1860             _Innovation = Ynpu - _HXa
1861         #
1862         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1863         # ---> avec analysis
1864         selfA.StoredVariables["Analysis"].store( Xa )
1865         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
1866             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
1867         if selfA._toStore("InnovationAtCurrentAnalysis"):
1868             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
1869         # ---> avec current state
1870         if selfA._parameters["StoreInternalVariables"] \
1871             or selfA._toStore("CurrentState"):
1872             selfA.StoredVariables["CurrentState"].store( Xn )
1873         if selfA._toStore("ForecastState"):
1874             selfA.StoredVariables["ForecastState"].store( EMX )
1875         if selfA._toStore("BMA"):
1876             selfA.StoredVariables["BMA"].store( EMX - Xa )
1877         if selfA._toStore("InnovationAtCurrentState"):
1878             selfA.StoredVariables["InnovationAtCurrentState"].store( - HE2 + Ynpu )
1879         if selfA._toStore("SimulatedObservationAtCurrentState") \
1880             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1881             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HE2 )
1882         # ---> autres
1883         if selfA._parameters["StoreInternalVariables"] \
1884             or selfA._toStore("CostFunctionJ") \
1885             or selfA._toStore("CostFunctionJb") \
1886             or selfA._toStore("CostFunctionJo") \
1887             or selfA._toStore("CurrentOptimum") \
1888             or selfA._toStore("APosterioriCovariance"):
1889             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
1890             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
1891             J   = Jb + Jo
1892             selfA.StoredVariables["CostFunctionJb"].store( Jb )
1893             selfA.StoredVariables["CostFunctionJo"].store( Jo )
1894             selfA.StoredVariables["CostFunctionJ" ].store( J )
1895             #
1896             if selfA._toStore("IndexOfOptimum") \
1897                 or selfA._toStore("CurrentOptimum") \
1898                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
1899                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
1900                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
1901                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1902                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
1903             if selfA._toStore("IndexOfOptimum"):
1904                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
1905             if selfA._toStore("CurrentOptimum"):
1906                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
1907             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
1908                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
1909             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
1910                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
1911             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
1912                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
1913             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
1914                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
1915         if selfA._toStore("APosterioriCovariance"):
1916             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
1917         if selfA._parameters["EstimationOf"] == "Parameters" \
1918             and J < previousJMinimum:
1919             previousJMinimum    = J
1920             XaMin               = Xa
1921             if selfA._toStore("APosterioriCovariance"):
1922                 covarianceXaMin = Pn
1923     #
1924     # Stockage final supplémentaire de l'optimum en estimation de paramètres
1925     # ----------------------------------------------------------------------
1926     if selfA._parameters["EstimationOf"] == "Parameters":
1927         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1928         selfA.StoredVariables["Analysis"].store( XaMin )
1929         if selfA._toStore("APosterioriCovariance"):
1930             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
1931         if selfA._toStore("BMA"):
1932             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
1933     #
1934     return 0
1935
1936 # ==============================================================================
1937 def mmqr(
1938         func     = None,
1939         x0       = None,
1940         fprime   = None,
1941         bounds   = None,
1942         quantile = 0.5,
1943         maxfun   = 15000,
1944         toler    = 1.e-06,
1945         y        = None,
1946         ):
1947     """
1948     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
1949     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
1950     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
1951     """
1952     #
1953     # Recuperation des donnees et informations initiales
1954     # --------------------------------------------------
1955     variables = numpy.ravel( x0 )
1956     mesures   = numpy.ravel( y )
1957     increment = sys.float_info[0]
1958     p         = variables.size
1959     n         = mesures.size
1960     quantile  = float(quantile)
1961     #
1962     # Calcul des parametres du MM
1963     # ---------------------------
1964     tn      = float(toler) / n
1965     e0      = -tn / math.log(tn)
1966     epsilon = (e0-tn)/(1+math.log(e0))
1967     #
1968     # Calculs d'initialisation
1969     # ------------------------
1970     residus  = mesures - numpy.ravel( func( variables ) )
1971     poids    = 1./(epsilon+numpy.abs(residus))
1972     veps     = 1. - 2. * quantile - residus * poids
1973     lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
1974     iteration = 0
1975     #
1976     # Recherche iterative
1977     # -------------------
1978     while (increment > toler) and (iteration < maxfun) :
1979         iteration += 1
1980         #
1981         Derivees  = numpy.array(fprime(variables))
1982         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
1983         DeriveesT = Derivees.transpose()
1984         M         =   numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
1985         SM        =   numpy.transpose(numpy.dot( DeriveesT , veps ))
1986         step      = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
1987         #
1988         variables = variables + step
1989         if bounds is not None:
1990             # Attention : boucle infinie à éviter si un intervalle est trop petit
1991             while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
1992                 step      = step/2.
1993                 variables = variables - step
1994         residus   = mesures - numpy.ravel( func(variables) )
1995         surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
1996         #
1997         while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
1998             step      = step/2.
1999             variables = variables - step
2000             residus   = mesures - numpy.ravel( func(variables) )
2001             surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
2002         #
2003         increment     = lastsurrogate-surrogate
2004         poids         = 1./(epsilon+numpy.abs(residus))
2005         veps          = 1. - 2. * quantile - residus * poids
2006         lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
2007     #
2008     # Mesure d'écart
2009     # --------------
2010     Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
2011     #
2012     return variables, Ecart, [n,p,iteration,increment,0]
2013
2014 # ==============================================================================
2015 def multi3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle):
2016     """
2017     3DVAR multi-pas et multi-méthodes
2018     """
2019     #
2020     # Initialisation
2021     # --------------
2022     Xn = numpy.ravel(Xb).reshape((-1,1))
2023     #
2024     if selfA._parameters["EstimationOf"] == "State":
2025         M = EM["Direct"].appliedTo
2026         #
2027         if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
2028             selfA.StoredVariables["Analysis"].store( Xn )
2029             if selfA._toStore("APosterioriCovariance"):
2030                 if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(Xn.size)
2031                 else:                         Pn = B
2032                 selfA.StoredVariables["APosterioriCovariance"].store( Pn )
2033             if selfA._toStore("ForecastState"):
2034                 selfA.StoredVariables["ForecastState"].store( Xn )
2035     #
2036     if hasattr(Y,"stepnumber"):
2037         duration = Y.stepnumber()
2038     else:
2039         duration = 2
2040     #
2041     # Multi-pas
2042     # ---------
2043     for step in range(duration-1):
2044         if hasattr(Y,"store"):
2045             Ynpu = numpy.ravel( Y[step+1] ).reshape((-1,1))
2046         else:
2047             Ynpu = numpy.ravel( Y ).reshape((-1,1))
2048         #
2049         if selfA._parameters["EstimationOf"] == "State": # Forecast
2050             Xn = selfA.StoredVariables["Analysis"][-1]
2051             Xn_predicted = M( Xn )
2052             if selfA._toStore("ForecastState"):
2053                 selfA.StoredVariables["ForecastState"].store( Xn_predicted )
2054         elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
2055             # --- > Par principe, M = Id, Q = 0
2056             Xn_predicted = Xn
2057         Xn_predicted = numpy.ravel(Xn_predicted).reshape((-1,1))
2058         #
2059         oneCycle(selfA, Xn_predicted, Ynpu, U, HO, None, None, R, B, None)
2060     #
2061     return 0
2062
2063 # ==============================================================================
2064 def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
2065     """
2066     3DVAR PSAS
2067     """
2068     #
2069     # Initialisations
2070     # ---------------
2071     #
2072     # Opérateurs
2073     Hm = HO["Direct"].appliedTo
2074     #
2075     # Utilisation éventuelle d'un vecteur H(Xb) précalculé
2076     if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
2077         HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
2078     else:
2079         HXb = Hm( Xb )
2080     HXb = numpy.asmatrix(numpy.ravel( HXb )).T
2081     if Y.size != HXb.size:
2082         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))
2083     if max(Y.shape) != max(HXb.shape):
2084         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))
2085     #
2086     if selfA._toStore("JacobianMatrixAtBackground"):
2087         HtMb = HO["Tangent"].asMatrix(ValueForMethodForm = Xb)
2088         HtMb = HtMb.reshape(Y.size,Xb.size) # ADAO & check shape
2089         selfA.StoredVariables["JacobianMatrixAtBackground"].store( HtMb )
2090     #
2091     Ht = HO["Tangent"].asMatrix(Xb)
2092     BHT = B * Ht.T
2093     HBHTpR = R + Ht * BHT
2094     Innovation = Y - HXb
2095     #
2096     # Point de démarrage de l'optimisation
2097     Xini = numpy.zeros(Xb.shape)
2098     #
2099     # Définition de la fonction-coût
2100     # ------------------------------
2101     def CostFunction(w):
2102         _W = numpy.asmatrix(numpy.ravel( w )).T
2103         if selfA._parameters["StoreInternalVariables"] or \
2104             selfA._toStore("CurrentState") or \
2105             selfA._toStore("CurrentOptimum"):
2106             selfA.StoredVariables["CurrentState"].store( Xb + BHT * _W )
2107         if selfA._toStore("SimulatedObservationAtCurrentState") or \
2108             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2109             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( Hm( Xb + BHT * _W ) )
2110         if selfA._toStore("InnovationAtCurrentState"):
2111             selfA.StoredVariables["InnovationAtCurrentState"].store( Innovation )
2112         #
2113         Jb  = float( 0.5 * _W.T * HBHTpR * _W )
2114         Jo  = float( - _W.T * Innovation )
2115         J   = Jb + Jo
2116         #
2117         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
2118         selfA.StoredVariables["CostFunctionJb"].store( Jb )
2119         selfA.StoredVariables["CostFunctionJo"].store( Jo )
2120         selfA.StoredVariables["CostFunctionJ" ].store( J )
2121         if selfA._toStore("IndexOfOptimum") or \
2122             selfA._toStore("CurrentOptimum") or \
2123             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
2124             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
2125             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
2126             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2127             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2128         if selfA._toStore("IndexOfOptimum"):
2129             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2130         if selfA._toStore("CurrentOptimum"):
2131             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
2132         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2133             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
2134         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2135             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2136         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2137             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2138         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2139             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2140         return J
2141     #
2142     def GradientOfCostFunction(w):
2143         _W = numpy.asmatrix(numpy.ravel( w )).T
2144         GradJb  = HBHTpR * _W
2145         GradJo  = - Innovation
2146         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
2147         return GradJ
2148     #
2149     # Minimisation de la fonctionnelle
2150     # --------------------------------
2151     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
2152     #
2153     if selfA._parameters["Minimizer"] == "LBFGSB":
2154         if "0.19" <= scipy.version.version <= "1.1.0":
2155             import lbfgsbhlt as optimiseur
2156         else:
2157             import scipy.optimize as optimiseur
2158         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
2159             func        = CostFunction,
2160             x0          = Xini,
2161             fprime      = GradientOfCostFunction,
2162             args        = (),
2163             bounds      = selfA._parameters["Bounds"],
2164             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
2165             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
2166             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2167             iprint      = selfA._parameters["optiprint"],
2168             )
2169         nfeval = Informations['funcalls']
2170         rc     = Informations['warnflag']
2171     elif selfA._parameters["Minimizer"] == "TNC":
2172         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
2173             func        = CostFunction,
2174             x0          = Xini,
2175             fprime      = GradientOfCostFunction,
2176             args        = (),
2177             bounds      = selfA._parameters["Bounds"],
2178             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
2179             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2180             ftol        = selfA._parameters["CostDecrementTolerance"],
2181             messages    = selfA._parameters["optmessages"],
2182             )
2183     elif selfA._parameters["Minimizer"] == "CG":
2184         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
2185             f           = CostFunction,
2186             x0          = Xini,
2187             fprime      = GradientOfCostFunction,
2188             args        = (),
2189             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2190             gtol        = selfA._parameters["GradientNormTolerance"],
2191             disp        = selfA._parameters["optdisp"],
2192             full_output = True,
2193             )
2194     elif selfA._parameters["Minimizer"] == "NCG":
2195         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
2196             f           = CostFunction,
2197             x0          = Xini,
2198             fprime      = GradientOfCostFunction,
2199             args        = (),
2200             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2201             avextol     = selfA._parameters["CostDecrementTolerance"],
2202             disp        = selfA._parameters["optdisp"],
2203             full_output = True,
2204             )
2205     elif selfA._parameters["Minimizer"] == "BFGS":
2206         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
2207             f           = CostFunction,
2208             x0          = Xini,
2209             fprime      = GradientOfCostFunction,
2210             args        = (),
2211             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2212             gtol        = selfA._parameters["GradientNormTolerance"],
2213             disp        = selfA._parameters["optdisp"],
2214             full_output = True,
2215             )
2216     else:
2217         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
2218     #
2219     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2220     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
2221     #
2222     # Correction pour pallier a un bug de TNC sur le retour du Minimum
2223     # ----------------------------------------------------------------
2224     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
2225         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
2226         Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
2227     else:
2228         Minimum = Xb + BHT * numpy.asmatrix(numpy.ravel( Minimum )).T
2229     #
2230     # Obtention de l'analyse
2231     # ----------------------
2232     Xa = Minimum
2233     #
2234     selfA.StoredVariables["Analysis"].store( Xa )
2235     #
2236     if selfA._toStore("OMA") or \
2237         selfA._toStore("SigmaObs2") or \
2238         selfA._toStore("SimulationQuantiles") or \
2239         selfA._toStore("SimulatedObservationAtOptimum"):
2240         if selfA._toStore("SimulatedObservationAtCurrentState"):
2241             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
2242         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2243             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
2244         else:
2245             HXa = Hm( Xa )
2246     #
2247     # Calcul de la covariance d'analyse
2248     # ---------------------------------
2249     if selfA._toStore("APosterioriCovariance") or \
2250         selfA._toStore("SimulationQuantiles") or \
2251         selfA._toStore("JacobianMatrixAtOptimum") or \
2252         selfA._toStore("KalmanGainAtOptimum"):
2253         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
2254         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
2255     if selfA._toStore("APosterioriCovariance") or \
2256         selfA._toStore("SimulationQuantiles") or \
2257         selfA._toStore("KalmanGainAtOptimum"):
2258         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
2259         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
2260     if selfA._toStore("APosterioriCovariance") or \
2261         selfA._toStore("SimulationQuantiles"):
2262         BI = B.getI()
2263         RI = R.getI()
2264         HessienneI = []
2265         nb = Xa.size
2266         for i in range(nb):
2267             _ee    = numpy.matrix(numpy.zeros(nb)).T
2268             _ee[i] = 1.
2269             _HtEE  = numpy.dot(HtM,_ee)
2270             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
2271             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
2272         HessienneI = numpy.matrix( HessienneI )
2273         A = HessienneI.I
2274         if min(A.shape) != max(A.shape):
2275             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)))
2276         if (numpy.diag(A) < 0).any():
2277             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,))
2278         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
2279             try:
2280                 L = numpy.linalg.cholesky( A )
2281             except:
2282                 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,))
2283     if selfA._toStore("APosterioriCovariance"):
2284         selfA.StoredVariables["APosterioriCovariance"].store( A )
2285     if selfA._toStore("JacobianMatrixAtOptimum"):
2286         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
2287     if selfA._toStore("KalmanGainAtOptimum"):
2288         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
2289         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
2290         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
2291     #
2292     # Calculs et/ou stockages supplémentaires
2293     # ---------------------------------------
2294     if selfA._toStore("Innovation") or \
2295         selfA._toStore("SigmaObs2") or \
2296         selfA._toStore("MahalanobisConsistency") or \
2297         selfA._toStore("OMB"):
2298         d  = Y - HXb
2299     if selfA._toStore("Innovation"):
2300         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
2301     if selfA._toStore("BMA"):
2302         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
2303     if selfA._toStore("OMA"):
2304         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
2305     if selfA._toStore("OMB"):
2306         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
2307     if selfA._toStore("SigmaObs2"):
2308         TraceR = R.trace(Y.size)
2309         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
2310     if selfA._toStore("MahalanobisConsistency"):
2311         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
2312     if selfA._toStore("SimulationQuantiles"):
2313         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
2314     if selfA._toStore("SimulatedObservationAtBackground"):
2315         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
2316     if selfA._toStore("SimulatedObservationAtOptimum"):
2317         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
2318     #
2319     return 0
2320
2321 # ==============================================================================
2322 def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"):
2323     """
2324     Stochastic EnKF
2325     """
2326     if selfA._parameters["EstimationOf"] == "Parameters":
2327         selfA._parameters["StoreInternalVariables"] = True
2328     #
2329     # Opérateurs
2330     H = HO["Direct"].appliedControledFormTo
2331     #
2332     if selfA._parameters["EstimationOf"] == "State":
2333         M = EM["Direct"].appliedControledFormTo
2334     #
2335     if CM is not None and "Tangent" in CM and U is not None:
2336         Cm = CM["Tangent"].asMatrix(Xb)
2337     else:
2338         Cm = None
2339     #
2340     # Durée d'observation et tailles
2341     if hasattr(Y,"stepnumber"):
2342         duration = Y.stepnumber()
2343         __p = numpy.cumprod(Y.shape())[-1]
2344     else:
2345         duration = 2
2346         __p = numpy.array(Y).size
2347     #
2348     # Précalcul des inversions de B et R
2349     if selfA._parameters["StoreInternalVariables"] \
2350         or selfA._toStore("CostFunctionJ") \
2351         or selfA._toStore("CostFunctionJb") \
2352         or selfA._toStore("CostFunctionJo") \
2353         or selfA._toStore("CurrentOptimum") \
2354         or selfA._toStore("APosterioriCovariance"):
2355         BI = B.getI()
2356         RI = R.getI()
2357     #
2358     __n = Xb.size
2359     __m = selfA._parameters["NumberOfMembers"]
2360     #
2361     if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
2362     else:                         Pn = B
2363     if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
2364     else:                         Rn = R
2365     Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m )
2366     #
2367     if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
2368         selfA.StoredVariables["Analysis"].store( Xb )
2369         if selfA._toStore("APosterioriCovariance"):
2370             selfA.StoredVariables["APosterioriCovariance"].store( Pn )
2371             covarianceXa = Pn
2372     #
2373     previousJMinimum = numpy.finfo(float).max
2374     #
2375     for step in range(duration-1):
2376         if hasattr(Y,"store"):
2377             Ynpu = numpy.ravel( Y[step+1] ).reshape((__p,1))
2378         else:
2379             Ynpu = numpy.ravel( Y ).reshape((__p,1))
2380         #
2381         if U is not None:
2382             if hasattr(U,"store") and len(U)>1:
2383                 Un = numpy.asmatrix(numpy.ravel( U[step] )).T
2384             elif hasattr(U,"store") and len(U)==1:
2385                 Un = numpy.asmatrix(numpy.ravel( U[0] )).T
2386             else:
2387                 Un = numpy.asmatrix(numpy.ravel( U )).T
2388         else:
2389             Un = None
2390         #
2391         if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
2392             Xn = CovarianceInflation( Xn,
2393                 selfA._parameters["InflationType"],
2394                 selfA._parameters["InflationFactor"],
2395                 )
2396         #
2397         if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
2398             EMX = M( [(Xn[:,i], Un) for i in range(__m)],
2399                 argsAsSerie = True,
2400                 returnSerieAsArrayMatrix = True )
2401             Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q )
2402             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
2403                 argsAsSerie = True,
2404                 returnSerieAsArrayMatrix = True )
2405             if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
2406                 Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
2407                 Xn_predicted = Xn_predicted + Cm * Un
2408         elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
2409             # --- > Par principe, M = Id, Q = 0
2410             Xn_predicted = Xn
2411             HX_predicted = H( [(Xn_predicted[:,i], Un) for i in range(__m)],
2412                 argsAsSerie = True,
2413                 returnSerieAsArrayMatrix = True )
2414         #
2415         # Mean of forecast and observation of forecast
2416         Xfm  = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
2417         Hfm  = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
2418         #
2419         #--------------------------
2420         if VariantM == "KalmanFilterFormula05":
2421             PfHT, HPfHT = 0., 0.
2422             for i in range(__m):
2423                 Exfi = Xn_predicted[:,i].reshape((__n,1)) - Xfm
2424                 Eyfi = HX_predicted[:,i].reshape((__p,1)) - Hfm
2425                 PfHT  += Exfi * Eyfi.T
2426                 HPfHT += Eyfi * Eyfi.T
2427             PfHT  = (1./(__m-1)) * PfHT
2428             HPfHT = (1./(__m-1)) * HPfHT
2429             Kn     = PfHT * ( R + HPfHT ).I
2430             del PfHT, HPfHT
2431             #
2432             for i in range(__m):
2433                 ri = numpy.random.multivariate_normal(numpy.zeros(__p), Rn)
2434                 Xn[:,i] = numpy.ravel(Xn_predicted[:,i]) + Kn @ (numpy.ravel(Ynpu) + ri - HX_predicted[:,i])
2435         #--------------------------
2436         elif VariantM == "KalmanFilterFormula16":
2437             EpY   = EnsembleOfCenteredPerturbations(Ynpu, Rn, __m)
2438             EpYm  = EpY.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1))
2439             #
2440             EaX   = EnsembleOfAnomalies( Xn_predicted ) / math.sqrt(__m-1)
2441             EaY = (HX_predicted - Hfm - EpY + EpYm) / math.sqrt(__m-1)
2442             #
2443             Kn = EaX @ EaY.T @ numpy.linalg.inv( EaY @ EaY.T)
2444             #
2445             for i in range(__m):
2446                 Xn[:,i] = numpy.ravel(Xn_predicted[:,i]) + Kn @ (numpy.ravel(EpY[:,i]) - HX_predicted[:,i])
2447         #--------------------------
2448         else:
2449             raise ValueError("VariantM has to be chosen in the authorized methods list.")
2450         #
2451         if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
2452             Xn = CovarianceInflation( Xn,
2453                 selfA._parameters["InflationType"],
2454                 selfA._parameters["InflationFactor"],
2455                 )
2456         #
2457         Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1))
2458         #--------------------------
2459         #
2460         if selfA._parameters["StoreInternalVariables"] \
2461             or selfA._toStore("CostFunctionJ") \
2462             or selfA._toStore("CostFunctionJb") \
2463             or selfA._toStore("CostFunctionJo") \
2464             or selfA._toStore("APosterioriCovariance") \
2465             or selfA._toStore("InnovationAtCurrentAnalysis") \
2466             or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
2467             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2468             _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
2469             _Innovation = Ynpu - _HXa
2470         #
2471         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2472         # ---> avec analysis
2473         selfA.StoredVariables["Analysis"].store( Xa )
2474         if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
2475             selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
2476         if selfA._toStore("InnovationAtCurrentAnalysis"):
2477             selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
2478         # ---> avec current state
2479         if selfA._parameters["StoreInternalVariables"] \
2480             or selfA._toStore("CurrentState"):
2481             selfA.StoredVariables["CurrentState"].store( Xn )
2482         if selfA._toStore("ForecastState"):
2483             selfA.StoredVariables["ForecastState"].store( EMX )
2484         if selfA._toStore("BMA"):
2485             selfA.StoredVariables["BMA"].store( EMX - Xa )
2486         if selfA._toStore("InnovationAtCurrentState"):
2487             selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
2488         if selfA._toStore("SimulatedObservationAtCurrentState") \
2489             or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2490             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
2491         # ---> autres
2492         if selfA._parameters["StoreInternalVariables"] \
2493             or selfA._toStore("CostFunctionJ") \
2494             or selfA._toStore("CostFunctionJb") \
2495             or selfA._toStore("CostFunctionJo") \
2496             or selfA._toStore("CurrentOptimum") \
2497             or selfA._toStore("APosterioriCovariance"):
2498             Jb  = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
2499             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
2500             J   = Jb + Jo
2501             selfA.StoredVariables["CostFunctionJb"].store( Jb )
2502             selfA.StoredVariables["CostFunctionJo"].store( Jo )
2503             selfA.StoredVariables["CostFunctionJ" ].store( J )
2504             #
2505             if selfA._toStore("IndexOfOptimum") \
2506                 or selfA._toStore("CurrentOptimum") \
2507                 or selfA._toStore("CostFunctionJAtCurrentOptimum") \
2508                 or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
2509                 or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
2510                 or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2511                 IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2512             if selfA._toStore("IndexOfOptimum"):
2513                 selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2514             if selfA._toStore("CurrentOptimum"):
2515                 selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
2516             if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2517                 selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
2518             if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2519                 selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2520             if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2521                 selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2522             if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2523                 selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2524         if selfA._toStore("APosterioriCovariance"):
2525             selfA.StoredVariables["APosterioriCovariance"].store( EnsembleErrorCovariance(Xn) )
2526         if selfA._parameters["EstimationOf"] == "Parameters" \
2527             and J < previousJMinimum:
2528             previousJMinimum    = J
2529             XaMin               = Xa
2530             if selfA._toStore("APosterioriCovariance"):
2531                 covarianceXaMin = Pn
2532     #
2533     # Stockage final supplémentaire de l'optimum en estimation de paramètres
2534     # ----------------------------------------------------------------------
2535     if selfA._parameters["EstimationOf"] == "Parameters":
2536         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
2537         selfA.StoredVariables["Analysis"].store( XaMin )
2538         if selfA._toStore("APosterioriCovariance"):
2539             selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
2540         if selfA._toStore("BMA"):
2541             selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
2542     #
2543     return 0
2544
2545 # ==============================================================================
2546 def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
2547     """
2548     3DVAR
2549     """
2550     #
2551     # Initialisations
2552     # ---------------
2553     #
2554     # Opérateurs
2555     Hm = HO["Direct"].appliedTo
2556     Ha = HO["Adjoint"].appliedInXTo
2557     #
2558     # Utilisation éventuelle d'un vecteur H(Xb) précalculé
2559     if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
2560         HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
2561     else:
2562         HXb = Hm( Xb )
2563     HXb = numpy.asmatrix(numpy.ravel( HXb )).T
2564     if Y.size != HXb.size:
2565         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))
2566     if max(Y.shape) != max(HXb.shape):
2567         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))
2568     #
2569     if selfA._toStore("JacobianMatrixAtBackground"):
2570         HtMb = HO["Tangent"].asMatrix(ValueForMethodForm = Xb)
2571         HtMb = HtMb.reshape(Y.size,Xb.size) # ADAO & check shape
2572         selfA.StoredVariables["JacobianMatrixAtBackground"].store( HtMb )
2573     #
2574     # Précalcul des inversions de B et R
2575     BI = B.getI()
2576     RI = R.getI()
2577     #
2578     # Point de démarrage de l'optimisation
2579     Xini = selfA._parameters["InitializationPoint"]
2580     #
2581     # Définition de la fonction-coût
2582     # ------------------------------
2583     def CostFunction(x):
2584         _X  = numpy.asmatrix(numpy.ravel( x )).T
2585         if selfA._parameters["StoreInternalVariables"] or \
2586             selfA._toStore("CurrentState") or \
2587             selfA._toStore("CurrentOptimum"):
2588             selfA.StoredVariables["CurrentState"].store( _X )
2589         _HX = Hm( _X )
2590         _HX = numpy.asmatrix(numpy.ravel( _HX )).T
2591         _Innovation = Y - _HX
2592         if selfA._toStore("SimulatedObservationAtCurrentState") or \
2593             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2594             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
2595         if selfA._toStore("InnovationAtCurrentState"):
2596             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
2597         #
2598         Jb  = float( 0.5 * (_X - Xb).T * BI * (_X - Xb) )
2599         Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
2600         J   = Jb + Jo
2601         #
2602         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
2603         selfA.StoredVariables["CostFunctionJb"].store( Jb )
2604         selfA.StoredVariables["CostFunctionJo"].store( Jo )
2605         selfA.StoredVariables["CostFunctionJ" ].store( J )
2606         if selfA._toStore("IndexOfOptimum") or \
2607             selfA._toStore("CurrentOptimum") or \
2608             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
2609             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
2610             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
2611             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2612             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2613         if selfA._toStore("IndexOfOptimum"):
2614             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2615         if selfA._toStore("CurrentOptimum"):
2616             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
2617         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2618             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
2619         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2620             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2621         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2622             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2623         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2624             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2625         return J
2626     #
2627     def GradientOfCostFunction(x):
2628         _X      = numpy.asmatrix(numpy.ravel( x )).T
2629         _HX     = Hm( _X )
2630         _HX     = numpy.asmatrix(numpy.ravel( _HX )).T
2631         GradJb  = BI * (_X - Xb)
2632         GradJo  = - Ha( (_X, RI * (Y - _HX)) )
2633         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
2634         return GradJ
2635     #
2636     # Minimisation de la fonctionnelle
2637     # --------------------------------
2638     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
2639     #
2640     if selfA._parameters["Minimizer"] == "LBFGSB":
2641         if "0.19" <= scipy.version.version <= "1.1.0":
2642             import lbfgsbhlt as optimiseur
2643         else:
2644             import scipy.optimize as optimiseur
2645         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
2646             func        = CostFunction,
2647             x0          = Xini,
2648             fprime      = GradientOfCostFunction,
2649             args        = (),
2650             bounds      = selfA._parameters["Bounds"],
2651             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
2652             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
2653             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2654             iprint      = selfA._parameters["optiprint"],
2655             )
2656         nfeval = Informations['funcalls']
2657         rc     = Informations['warnflag']
2658     elif selfA._parameters["Minimizer"] == "TNC":
2659         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
2660             func        = CostFunction,
2661             x0          = Xini,
2662             fprime      = GradientOfCostFunction,
2663             args        = (),
2664             bounds      = selfA._parameters["Bounds"],
2665             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
2666             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2667             ftol        = selfA._parameters["CostDecrementTolerance"],
2668             messages    = selfA._parameters["optmessages"],
2669             )
2670     elif selfA._parameters["Minimizer"] == "CG":
2671         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
2672             f           = CostFunction,
2673             x0          = Xini,
2674             fprime      = GradientOfCostFunction,
2675             args        = (),
2676             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2677             gtol        = selfA._parameters["GradientNormTolerance"],
2678             disp        = selfA._parameters["optdisp"],
2679             full_output = True,
2680             )
2681     elif selfA._parameters["Minimizer"] == "NCG":
2682         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
2683             f           = CostFunction,
2684             x0          = Xini,
2685             fprime      = GradientOfCostFunction,
2686             args        = (),
2687             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2688             avextol     = selfA._parameters["CostDecrementTolerance"],
2689             disp        = selfA._parameters["optdisp"],
2690             full_output = True,
2691             )
2692     elif selfA._parameters["Minimizer"] == "BFGS":
2693         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
2694             f           = CostFunction,
2695             x0          = Xini,
2696             fprime      = GradientOfCostFunction,
2697             args        = (),
2698             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2699             gtol        = selfA._parameters["GradientNormTolerance"],
2700             disp        = selfA._parameters["optdisp"],
2701             full_output = True,
2702             )
2703     else:
2704         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
2705     #
2706     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2707     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
2708     #
2709     # Correction pour pallier a un bug de TNC sur le retour du Minimum
2710     # ----------------------------------------------------------------
2711     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
2712         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
2713     #
2714     # Obtention de l'analyse
2715     # ----------------------
2716     Xa = numpy.asmatrix(numpy.ravel( Minimum )).T
2717     #
2718     selfA.StoredVariables["Analysis"].store( Xa )
2719     #
2720     if selfA._toStore("OMA") or \
2721         selfA._toStore("SigmaObs2") or \
2722         selfA._toStore("SimulationQuantiles") or \
2723         selfA._toStore("SimulatedObservationAtOptimum"):
2724         if selfA._toStore("SimulatedObservationAtCurrentState"):
2725             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
2726         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
2727             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
2728         else:
2729             HXa = Hm( Xa )
2730     #
2731     # Calcul de la covariance d'analyse
2732     # ---------------------------------
2733     if selfA._toStore("APosterioriCovariance") or \
2734         selfA._toStore("SimulationQuantiles") or \
2735         selfA._toStore("JacobianMatrixAtOptimum") or \
2736         selfA._toStore("KalmanGainAtOptimum"):
2737         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
2738         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
2739     if selfA._toStore("APosterioriCovariance") or \
2740         selfA._toStore("SimulationQuantiles") or \
2741         selfA._toStore("KalmanGainAtOptimum"):
2742         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
2743         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
2744     if selfA._toStore("APosterioriCovariance") or \
2745         selfA._toStore("SimulationQuantiles"):
2746         HessienneI = []
2747         nb = Xa.size
2748         for i in range(nb):
2749             _ee    = numpy.matrix(numpy.zeros(nb)).T
2750             _ee[i] = 1.
2751             _HtEE  = numpy.dot(HtM,_ee)
2752             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
2753             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
2754         HessienneI = numpy.matrix( HessienneI )
2755         A = HessienneI.I
2756         if min(A.shape) != max(A.shape):
2757             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)))
2758         if (numpy.diag(A) < 0).any():
2759             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,))
2760         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
2761             try:
2762                 L = numpy.linalg.cholesky( A )
2763             except:
2764                 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,))
2765     if selfA._toStore("APosterioriCovariance"):
2766         selfA.StoredVariables["APosterioriCovariance"].store( A )
2767     if selfA._toStore("JacobianMatrixAtOptimum"):
2768         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
2769     if selfA._toStore("KalmanGainAtOptimum"):
2770         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
2771         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
2772         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
2773     #
2774     # Calculs et/ou stockages supplémentaires
2775     # ---------------------------------------
2776     if selfA._toStore("Innovation") or \
2777         selfA._toStore("SigmaObs2") or \
2778         selfA._toStore("MahalanobisConsistency") or \
2779         selfA._toStore("OMB"):
2780         d  = Y - HXb
2781     if selfA._toStore("Innovation"):
2782         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
2783     if selfA._toStore("BMA"):
2784         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
2785     if selfA._toStore("OMA"):
2786         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
2787     if selfA._toStore("OMB"):
2788         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
2789     if selfA._toStore("SigmaObs2"):
2790         TraceR = R.trace(Y.size)
2791         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
2792     if selfA._toStore("MahalanobisConsistency"):
2793         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
2794     if selfA._toStore("SimulationQuantiles"):
2795         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
2796     if selfA._toStore("SimulatedObservationAtBackground"):
2797         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
2798     if selfA._toStore("SimulatedObservationAtOptimum"):
2799         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
2800     #
2801     return 0
2802
2803 # ==============================================================================
2804 def std4dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
2805     """
2806     4DVAR
2807     """
2808     #
2809     # Initialisations
2810     # ---------------
2811     #
2812     # Opérateurs
2813     Hm = HO["Direct"].appliedControledFormTo
2814     Mm = EM["Direct"].appliedControledFormTo
2815     #
2816     if CM is not None and "Tangent" in CM and U is not None:
2817         Cm = CM["Tangent"].asMatrix(Xb)
2818     else:
2819         Cm = None
2820     #
2821     def Un(_step):
2822         if U is not None:
2823             if hasattr(U,"store") and 1<=_step<len(U) :
2824                 _Un = numpy.asmatrix(numpy.ravel( U[_step] )).T
2825             elif hasattr(U,"store") and len(U)==1:
2826                 _Un = numpy.asmatrix(numpy.ravel( U[0] )).T
2827             else:
2828                 _Un = numpy.asmatrix(numpy.ravel( U )).T
2829         else:
2830             _Un = None
2831         return _Un
2832     def CmUn(_xn,_un):
2833         if Cm is not None and _un is not None: # Attention : si Cm est aussi dans M, doublon !
2834             _Cm   = Cm.reshape(_xn.size,_un.size) # ADAO & check shape
2835             _CmUn = _Cm * _un
2836         else:
2837             _CmUn = 0.
2838         return _CmUn
2839     #
2840     # Remarque : les observations sont exploitées à partir du pas de temps
2841     # numéro 1, et sont utilisées dans Yo comme rangées selon ces indices.
2842     # Donc le pas 0 n'est pas utilisé puisque la première étape commence
2843     # avec l'observation du pas 1.
2844     #
2845     # Nombre de pas identique au nombre de pas d'observations
2846     if hasattr(Y,"stepnumber"):
2847         duration = Y.stepnumber()
2848     else:
2849         duration = 2
2850     #
2851     # Précalcul des inversions de B et R
2852     BI = B.getI()
2853     RI = R.getI()
2854     #
2855     # Point de démarrage de l'optimisation
2856     Xini = selfA._parameters["InitializationPoint"]
2857     #
2858     # Définition de la fonction-coût
2859     # ------------------------------
2860     selfA.DirectCalculation = [None,] # Le pas 0 n'est pas observé
2861     selfA.DirectInnovation  = [None,] # Le pas 0 n'est pas observé
2862     def CostFunction(x):
2863         _X  = numpy.asmatrix(numpy.ravel( x )).T
2864         if selfA._parameters["StoreInternalVariables"] or \
2865             selfA._toStore("CurrentState") or \
2866             selfA._toStore("CurrentOptimum"):
2867             selfA.StoredVariables["CurrentState"].store( _X )
2868         Jb  = float( 0.5 * (_X - Xb).T * BI * (_X - Xb) )
2869         selfA.DirectCalculation = [None,]
2870         selfA.DirectInnovation  = [None,]
2871         Jo  = 0.
2872         _Xn = _X
2873         for step in range(0,duration-1):
2874             if hasattr(Y,"store"):
2875                 _Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
2876             else:
2877                 _Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
2878             _Un = Un(step)
2879             #
2880             # Etape d'évolution
2881             if selfA._parameters["EstimationOf"] == "State":
2882                 _Xn = Mm( (_Xn, _Un) ) + CmUn(_Xn, _Un)
2883             elif selfA._parameters["EstimationOf"] == "Parameters":
2884                 pass
2885             #
2886             if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection":
2887                 _Xn = numpy.max(numpy.hstack((_Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,0])),axis=1)
2888                 _Xn = numpy.min(numpy.hstack((_Xn,numpy.asmatrix(selfA._parameters["Bounds"])[:,1])),axis=1)
2889             #
2890             # Etape de différence aux observations
2891             if selfA._parameters["EstimationOf"] == "State":
2892                 _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, None) ) )).T
2893             elif selfA._parameters["EstimationOf"] == "Parameters":
2894                 _YmHMX = _Ynpu - numpy.asmatrix(numpy.ravel( Hm( (_Xn, _Un) ) )).T - CmUn(_Xn, _Un)
2895             #
2896             # Stockage de l'état
2897             selfA.DirectCalculation.append( _Xn )
2898             selfA.DirectInnovation.append( _YmHMX )
2899             #
2900             # Ajout dans la fonctionnelle d'observation
2901             Jo = Jo + 0.5 * float( _YmHMX.T * RI * _YmHMX )
2902         J = Jb + Jo
2903         #
2904         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
2905         selfA.StoredVariables["CostFunctionJb"].store( Jb )
2906         selfA.StoredVariables["CostFunctionJo"].store( Jo )
2907         selfA.StoredVariables["CostFunctionJ" ].store( J )
2908         if selfA._toStore("IndexOfOptimum") or \
2909             selfA._toStore("CurrentOptimum") or \
2910             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
2911             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
2912             selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2913             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
2914         if selfA._toStore("IndexOfOptimum"):
2915             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
2916         if selfA._toStore("CurrentOptimum"):
2917             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
2918         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
2919             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
2920         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
2921             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
2922         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
2923             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
2924         return J
2925     #
2926     def GradientOfCostFunction(x):
2927         _X      = numpy.asmatrix(numpy.ravel( x )).T
2928         GradJb  = BI * (_X - Xb)
2929         GradJo  = 0.
2930         for step in range(duration-1,0,-1):
2931             # Étape de récupération du dernier stockage de l'évolution
2932             _Xn = selfA.DirectCalculation.pop()
2933             # Étape de récupération du dernier stockage de l'innovation
2934             _YmHMX = selfA.DirectInnovation.pop()
2935             # Calcul des adjoints
2936             Ha = HO["Adjoint"].asMatrix(ValueForMethodForm = _Xn)
2937             Ha = Ha.reshape(_Xn.size,_YmHMX.size) # ADAO & check shape
2938             Ma = EM["Adjoint"].asMatrix(ValueForMethodForm = _Xn)
2939             Ma = Ma.reshape(_Xn.size,_Xn.size) # ADAO & check shape
2940             # Calcul du gradient par état adjoint
2941             GradJo = GradJo + Ha * RI * _YmHMX # Équivaut pour Ha linéaire à : Ha( (_Xn, RI * _YmHMX) )
2942             GradJo = Ma * GradJo               # Équivaut pour Ma linéaire à : Ma( (_Xn, GradJo) )
2943         GradJ = numpy.ravel( GradJb ) - numpy.ravel( GradJo )
2944         return GradJ
2945     #
2946     # Minimisation de la fonctionnelle
2947     # --------------------------------
2948     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
2949     #
2950     if selfA._parameters["Minimizer"] == "LBFGSB":
2951         if "0.19" <= scipy.version.version <= "1.1.0":
2952             import lbfgsbhlt as optimiseur
2953         else:
2954             import scipy.optimize as optimiseur
2955         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
2956             func        = CostFunction,
2957             x0          = Xini,
2958             fprime      = GradientOfCostFunction,
2959             args        = (),
2960             bounds      = selfA._parameters["Bounds"],
2961             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
2962             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
2963             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2964             iprint      = selfA._parameters["optiprint"],
2965             )
2966         nfeval = Informations['funcalls']
2967         rc     = Informations['warnflag']
2968     elif selfA._parameters["Minimizer"] == "TNC":
2969         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
2970             func        = CostFunction,
2971             x0          = Xini,
2972             fprime      = GradientOfCostFunction,
2973             args        = (),
2974             bounds      = selfA._parameters["Bounds"],
2975             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
2976             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
2977             ftol        = selfA._parameters["CostDecrementTolerance"],
2978             messages    = selfA._parameters["optmessages"],
2979             )
2980     elif selfA._parameters["Minimizer"] == "CG":
2981         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
2982             f           = CostFunction,
2983             x0          = Xini,
2984             fprime      = GradientOfCostFunction,
2985             args        = (),
2986             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2987             gtol        = selfA._parameters["GradientNormTolerance"],
2988             disp        = selfA._parameters["optdisp"],
2989             full_output = True,
2990             )
2991     elif selfA._parameters["Minimizer"] == "NCG":
2992         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
2993             f           = CostFunction,
2994             x0          = Xini,
2995             fprime      = GradientOfCostFunction,
2996             args        = (),
2997             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
2998             avextol     = selfA._parameters["CostDecrementTolerance"],
2999             disp        = selfA._parameters["optdisp"],
3000             full_output = True,
3001             )
3002     elif selfA._parameters["Minimizer"] == "BFGS":
3003         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
3004             f           = CostFunction,
3005             x0          = Xini,
3006             fprime      = GradientOfCostFunction,
3007             args        = (),
3008             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3009             gtol        = selfA._parameters["GradientNormTolerance"],
3010             disp        = selfA._parameters["optdisp"],
3011             full_output = True,
3012             )
3013     else:
3014         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
3015     #
3016     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3017     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
3018     #
3019     # Correction pour pallier a un bug de TNC sur le retour du Minimum
3020     # ----------------------------------------------------------------
3021     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
3022         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
3023     #
3024     # Obtention de l'analyse
3025     # ----------------------
3026     Xa = numpy.asmatrix(numpy.ravel( Minimum )).T
3027     #
3028     selfA.StoredVariables["Analysis"].store( Xa )
3029     #
3030     # Calculs et/ou stockages supplémentaires
3031     # ---------------------------------------
3032     if selfA._toStore("BMA"):
3033         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
3034     #
3035     return 0
3036
3037 # ==============================================================================
3038 def van3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
3039     """
3040     3DVAR variational analysis with no inversion of B
3041     """
3042     #
3043     # Initialisations
3044     # ---------------
3045     #
3046     # Opérateurs
3047     Hm = HO["Direct"].appliedTo
3048     Ha = HO["Adjoint"].appliedInXTo
3049     #
3050     # Précalcul des inversions de B et R
3051     BT = B.getT()
3052     RI = R.getI()
3053     #
3054     # Point de démarrage de l'optimisation
3055     Xini = numpy.zeros(Xb.shape)
3056     #
3057     # Définition de la fonction-coût
3058     # ------------------------------
3059     def CostFunction(v):
3060         _V = numpy.asmatrix(numpy.ravel( v )).T
3061         _X = Xb + B * _V
3062         if selfA._parameters["StoreInternalVariables"] or \
3063             selfA._toStore("CurrentState") or \
3064             selfA._toStore("CurrentOptimum"):
3065             selfA.StoredVariables["CurrentState"].store( _X )
3066         _HX = Hm( _X )
3067         _HX = numpy.asmatrix(numpy.ravel( _HX )).T
3068         _Innovation = Y - _HX
3069         if selfA._toStore("SimulatedObservationAtCurrentState") or \
3070             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3071             selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
3072         if selfA._toStore("InnovationAtCurrentState"):
3073             selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
3074         #
3075         Jb  = float( 0.5 * _V.T * BT * _V )
3076         Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
3077         J   = Jb + Jo
3078         #
3079         selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) )
3080         selfA.StoredVariables["CostFunctionJb"].store( Jb )
3081         selfA.StoredVariables["CostFunctionJo"].store( Jo )
3082         selfA.StoredVariables["CostFunctionJ" ].store( J )
3083         if selfA._toStore("IndexOfOptimum") or \
3084             selfA._toStore("CurrentOptimum") or \
3085             selfA._toStore("CostFunctionJAtCurrentOptimum") or \
3086             selfA._toStore("CostFunctionJbAtCurrentOptimum") or \
3087             selfA._toStore("CostFunctionJoAtCurrentOptimum") or \
3088             selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3089             IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3090         if selfA._toStore("IndexOfOptimum"):
3091             selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
3092         if selfA._toStore("CurrentOptimum"):
3093             selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["CurrentState"][IndexMin] )
3094         if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3095             selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
3096         if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
3097             selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
3098         if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
3099             selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
3100         if selfA._toStore("CostFunctionJAtCurrentOptimum"):
3101             selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
3102         return J
3103     #
3104     def GradientOfCostFunction(v):
3105         _V = numpy.asmatrix(numpy.ravel( v )).T
3106         _X = Xb + B * _V
3107         _HX     = Hm( _X )
3108         _HX     = numpy.asmatrix(numpy.ravel( _HX )).T
3109         GradJb  = BT * _V
3110         GradJo  = - Ha( (_X, RI * (Y - _HX)) )
3111         GradJ   = numpy.ravel( GradJb ) + numpy.ravel( GradJo )
3112         return GradJ
3113     #
3114     # Minimisation de la fonctionnelle
3115     # --------------------------------
3116     nbPreviousSteps = selfA.StoredVariables["CostFunctionJ"].stepnumber()
3117     #
3118     if selfA._parameters["Minimizer"] == "LBFGSB":
3119         if "0.19" <= scipy.version.version <= "1.1.0":
3120             import lbfgsbhlt as optimiseur
3121         else:
3122             import scipy.optimize as optimiseur
3123         Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
3124             func        = CostFunction,
3125             x0          = Xini,
3126             fprime      = GradientOfCostFunction,
3127             args        = (),
3128             bounds      = selfA._parameters["Bounds"],
3129             maxfun      = selfA._parameters["MaximumNumberOfSteps"]-1,
3130             factr       = selfA._parameters["CostDecrementTolerance"]*1.e14,
3131             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3132             iprint      = selfA._parameters["optiprint"],
3133             )
3134         nfeval = Informations['funcalls']
3135         rc     = Informations['warnflag']
3136     elif selfA._parameters["Minimizer"] == "TNC":
3137         Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
3138             func        = CostFunction,
3139             x0          = Xini,
3140             fprime      = GradientOfCostFunction,
3141             args        = (),
3142             bounds      = selfA._parameters["Bounds"],
3143             maxfun      = selfA._parameters["MaximumNumberOfSteps"],
3144             pgtol       = selfA._parameters["ProjectedGradientTolerance"],
3145             ftol        = selfA._parameters["CostDecrementTolerance"],
3146             messages    = selfA._parameters["optmessages"],
3147             )
3148     elif selfA._parameters["Minimizer"] == "CG":
3149         Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
3150             f           = CostFunction,
3151             x0          = Xini,
3152             fprime      = GradientOfCostFunction,
3153             args        = (),
3154             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3155             gtol        = selfA._parameters["GradientNormTolerance"],
3156             disp        = selfA._parameters["optdisp"],
3157             full_output = True,
3158             )
3159     elif selfA._parameters["Minimizer"] == "NCG":
3160         Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
3161             f           = CostFunction,
3162             x0          = Xini,
3163             fprime      = GradientOfCostFunction,
3164             args        = (),
3165             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3166             avextol     = selfA._parameters["CostDecrementTolerance"],
3167             disp        = selfA._parameters["optdisp"],
3168             full_output = True,
3169             )
3170     elif selfA._parameters["Minimizer"] == "BFGS":
3171         Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
3172             f           = CostFunction,
3173             x0          = Xini,
3174             fprime      = GradientOfCostFunction,
3175             args        = (),
3176             maxiter     = selfA._parameters["MaximumNumberOfSteps"],
3177             gtol        = selfA._parameters["GradientNormTolerance"],
3178             disp        = selfA._parameters["optdisp"],
3179             full_output = True,
3180             )
3181     else:
3182         raise ValueError("Error in Minimizer name: %s"%selfA._parameters["Minimizer"])
3183     #
3184     IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
3185     MinJ     = selfA.StoredVariables["CostFunctionJ"][IndexMin]
3186     #
3187     # Correction pour pallier a un bug de TNC sur le retour du Minimum
3188     # ----------------------------------------------------------------
3189     if selfA._parameters["StoreInternalVariables"] or selfA._toStore("CurrentState"):
3190         Minimum = selfA.StoredVariables["CurrentState"][IndexMin]
3191         Minimum = numpy.asmatrix(numpy.ravel( Minimum )).T
3192     else:
3193         Minimum = Xb + B * numpy.asmatrix(numpy.ravel( Minimum )).T
3194     #
3195     # Obtention de l'analyse
3196     # ----------------------
3197     Xa = Minimum
3198     #
3199     selfA.StoredVariables["Analysis"].store( Xa )
3200     #
3201     if selfA._toStore("OMA") or \
3202         selfA._toStore("SigmaObs2") or \
3203         selfA._toStore("SimulationQuantiles") or \
3204         selfA._toStore("SimulatedObservationAtOptimum"):
3205         if selfA._toStore("SimulatedObservationAtCurrentState"):
3206             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
3207         elif selfA._toStore("SimulatedObservationAtCurrentOptimum"):
3208             HXa = selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
3209         else:
3210             HXa = Hm( Xa )
3211     #
3212     # Calcul de la covariance d'analyse
3213     # ---------------------------------
3214     if selfA._toStore("APosterioriCovariance") or \
3215         selfA._toStore("SimulationQuantiles") or \
3216         selfA._toStore("JacobianMatrixAtOptimum") or \
3217         selfA._toStore("KalmanGainAtOptimum"):
3218         HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
3219         HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
3220     if selfA._toStore("APosterioriCovariance") or \
3221         selfA._toStore("SimulationQuantiles") or \
3222         selfA._toStore("KalmanGainAtOptimum"):
3223         HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
3224         HaM = HaM.reshape(Xa.size,Y.size) # ADAO & check shape
3225     if selfA._toStore("APosterioriCovariance") or \
3226         selfA._toStore("SimulationQuantiles"):
3227         BI = B.getI()
3228         HessienneI = []
3229         nb = Xa.size
3230         for i in range(nb):
3231             _ee    = numpy.matrix(numpy.zeros(nb)).T
3232             _ee[i] = 1.
3233             _HtEE  = numpy.dot(HtM,_ee)
3234             _HtEE  = numpy.asmatrix(numpy.ravel( _HtEE )).T
3235             HessienneI.append( numpy.ravel( BI*_ee + HaM * (RI * _HtEE) ) )
3236         HessienneI = numpy.matrix( HessienneI )
3237         A = HessienneI.I
3238         if min(A.shape) != max(A.shape):
3239             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)))
3240         if (numpy.diag(A) < 0).any():
3241             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,))
3242         if logging.getLogger().level < logging.WARNING: # La verification n'a lieu qu'en debug
3243             try:
3244                 L = numpy.linalg.cholesky( A )
3245             except:
3246                 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,))
3247     if selfA._toStore("APosterioriCovariance"):
3248         selfA.StoredVariables["APosterioriCovariance"].store( A )
3249     if selfA._toStore("JacobianMatrixAtOptimum"):
3250         selfA.StoredVariables["JacobianMatrixAtOptimum"].store( HtM )
3251     if selfA._toStore("KalmanGainAtOptimum"):
3252         if   (Y.size <= Xb.size): KG  = B * HaM * (R + numpy.dot(HtM, B * HaM)).I
3253         elif (Y.size >  Xb.size): KG = (BI + numpy.dot(HaM, RI * HtM)).I * HaM * RI
3254         selfA.StoredVariables["KalmanGainAtOptimum"].store( KG )
3255     #
3256     # Calculs et/ou stockages supplémentaires
3257     # ---------------------------------------
3258     if selfA._toStore("Innovation") or \
3259         selfA._toStore("SigmaObs2") or \
3260         selfA._toStore("MahalanobisConsistency") or \
3261         selfA._toStore("OMB"):
3262         d  = Y - HXb
3263     if selfA._toStore("Innovation"):
3264         selfA.StoredVariables["Innovation"].store( numpy.ravel(d) )
3265     if selfA._toStore("BMA"):
3266         selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
3267     if selfA._toStore("OMA"):
3268         selfA.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
3269     if selfA._toStore("OMB"):
3270         selfA.StoredVariables["OMB"].store( numpy.ravel(d) )
3271     if selfA._toStore("SigmaObs2"):
3272         TraceR = R.trace(Y.size)
3273         selfA.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
3274     if selfA._toStore("MahalanobisConsistency"):
3275         selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
3276     if selfA._toStore("SimulationQuantiles"):
3277         QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM)
3278     if selfA._toStore("SimulatedObservationAtBackground"):
3279         selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
3280     if selfA._toStore("SimulatedObservationAtOptimum"):
3281         selfA.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
3282     #
3283     return 0
3284
3285 # ==============================================================================
3286 if __name__ == "__main__":
3287     print('\n AUTODIAGNOSTIC\n')