src/daComposant/daCore/NumericObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2024 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Définit les objets numériques génériques.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, copy, types, sys, logging, math, numpy, scipy, itertools, warnings
  29 import scipy.linalg  # Py3.6
  30 from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm
  31 from daCore.PlatformInfo import PlatformInfo, vt, vfloat
  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.ravel( X ).reshape((-1, 1))
  41     __sys_path_tmp = sys.path
  42     sys.path.insert(0, funcrepr["__userFunction__path"])
  43     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  44     __fonction = getattr(__module, funcrepr["__userFunction__name"])
  45     sys.path = __sys_path_tmp
  46     del __sys_path_tmp
  47     if isinstance(xArgs, dict):
  48         __HX  = __fonction( __X, **xArgs )
  49     else:
  50         __HX  = __fonction( __X )
  51     return numpy.ravel( __HX )
  52
  53 # ==============================================================================
  54 class FDApproximation(object):
  55     """
  56     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  57     création d'un objet, en fournissant une fonction "Function", on obtient un
  58     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  59     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  60     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  61     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  62     centrées si le booléen "centeredDF" est vrai.
  63     """
  64     __slots__ = (
  65         "__name", "__extraArgs", "__mpEnabled", "__mpWorkers", "__mfEnabled",
  66         "__rmEnabled", "__avoidRC", "__tolerBP", "__centeredDF", "__lengthRJ",
  67         "__listJPCP", "__listJPCI", "__listJPCR", "__listJPPN", "__listJPIN",
  68         "__userOperator", "__userFunction", "__increment", "__pool", "__dX",
  69         "__userFunction__name", "__userFunction__modl", "__userFunction__path",
  70     )
  71
  72     def __init__(self,
  73                  name                  = "FDApproximation",
  74                  Function              = None,
  75                  centeredDF            = False,
  76                  increment             = 0.01,
  77                  dX                    = None,
  78                  extraArguments        = None,
  79                  reducingMemoryUse     = False,
  80                  avoidingRedundancy    = True,
  81                  toleranceInRedundancy = 1.e-18,
  82                  lengthOfRedundancy    = -1,
  83                  mpEnabled             = False,
  84                  mpWorkers             = None,
  85                  mfEnabled             = False ):
  86         #
  87         self.__name = str(name)
  88         self.__extraArgs = extraArguments
  89         #
  90         if mpEnabled:
  91             try:
  92                 import multiprocessing  # noqa: F401
  93                 self.__mpEnabled = True
  94             except ImportError:
  95                 self.__mpEnabled = False
  96         else:
  97             self.__mpEnabled = False
  98         self.__mpWorkers = mpWorkers
  99         if self.__mpWorkers is not None and self.__mpWorkers < 1:
 100             self.__mpWorkers = None
 101         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled, self.__mpWorkers))
 102         #
 103         self.__mfEnabled = bool(mfEnabled)
 104         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
 105         #
 106         self.__rmEnabled = bool(reducingMemoryUse)
 107         logging.debug("FDA Calculs avec réduction mémoire : %s"%(self.__rmEnabled,))
 108         #
 109         if avoidingRedundancy:
 110             self.__avoidRC = True
 111             self.__tolerBP = float(toleranceInRedundancy)
 112             self.__lengthRJ = int(lengthOfRedundancy)
 113             self.__listJPCP = []  # Jacobian Previous Calculated Points
 114             self.__listJPCI = []  # Jacobian Previous Calculated Increment
 115             self.__listJPCR = []  # Jacobian Previous Calculated Results
 116             self.__listJPPN = []  # Jacobian Previous Calculated Point Norms
 117             self.__listJPIN = []  # Jacobian Previous Calculated Increment Norms
 118         else:
 119             self.__avoidRC = False
 120         logging.debug("FDA Calculs avec réduction des doublons : %s"%self.__avoidRC)
 121         if self.__avoidRC:
 122             logging.debug("FDA Tolérance de détermination des doublons : %.2e"%self.__tolerBP)
 123         #
 124         if self.__mpEnabled:
 125             if isinstance(Function, types.FunctionType):
 126                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 127                 self.__userFunction__name = Function.__name__
 128                 try:
 129                     mod = os.path.join(Function.__globals__['filepath'], Function.__globals__['filename'])
 130                 except Exception:
 131                     mod = os.path.abspath(Function.__globals__['__file__'])
 132                 if not os.path.isfile(mod):
 133                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 134                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc', '').replace('.pyo', '').replace('.py', '')
 135                 self.__userFunction__path = os.path.dirname(mod)
 136                 del mod
 137                 self.__userOperator = Operator(
 138                     name                 = self.__name,
 139                     fromMethod           = Function,
 140                     avoidingRedundancy   = self.__avoidRC,
 141                     inputAsMultiFunction = self.__mfEnabled,
 142                     extraArguments       = self.__extraArgs )
 143                 self.__userFunction = self.__userOperator.appliedTo  # Pour le calcul Direct
 144             elif isinstance(Function, types.MethodType):
 145                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 146                 self.__userFunction__name = Function.__name__
 147                 try:
 148                     mod = os.path.join(Function.__globals__['filepath'], Function.__globals__['filename'])
 149                 except Exception:
 150                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 151                 if not os.path.isfile(mod):
 152                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 153                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc', '').replace('.pyo', '').replace('.py', '')
 154                 self.__userFunction__path = os.path.dirname(mod)
 155                 del mod
 156                 self.__userOperator = Operator(
 157                     name                 = self.__name,
 158                     fromMethod           = Function,
 159                     avoidingRedundancy   = self.__avoidRC,
 160                     inputAsMultiFunction = self.__mfEnabled,
 161                     extraArguments       = self.__extraArgs )
 162                 self.__userFunction = self.__userOperator.appliedTo  # Pour le calcul Direct
 163             else:
 164                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 165         else:
 166             self.__userOperator = Operator(
 167                 name                 = self.__name,
 168                 fromMethod           = Function,
 169                 avoidingRedundancy   = self.__avoidRC,
 170                 inputAsMultiFunction = self.__mfEnabled,
 171                 extraArguments       = self.__extraArgs )
 172             self.__userFunction = self.__userOperator.appliedTo
 173         #
 174         self.__centeredDF = bool(centeredDF)
 175         if abs(float(increment)) > 1.e-15:
 176             self.__increment  = float(increment)
 177         else:
 178             self.__increment  = 0.01
 179         if dX is None:
 180             self.__dX     = None
 181         else:
 182             self.__dX     = numpy.ravel( dX )
 183
 184     # ---------------------------------------------------------
 185     def __doublon__(self, __e, __l, __n, __v=None):
 186         __ac, __iac = False, -1
 187         for i in range(len(__l) - 1, -1, -1):
 188             if numpy.linalg.norm(__e - __l[i]) < self.__tolerBP * __n[i]:
 189                 __ac, __iac = True, i
 190                 if __v is not None:
 191                     logging.debug("FDA Cas%s déjà calculé, récupération du doublon %i"%(__v, __iac))
 192                 break
 193         return __ac, __iac
 194
 195     # ---------------------------------------------------------
 196     def __listdotwith__(self, __LMatrix, __dotWith = None, __dotTWith = None):
 197         "Produit incrémental d'une matrice liste de colonnes avec un vecteur"
 198         if not isinstance(__LMatrix, (list, tuple)):
 199             raise TypeError("Columnwise list matrix has not the proper type: %s"%type(__LMatrix))
 200         if __dotWith is not None:
 201             __Idwx = numpy.ravel( __dotWith )
 202             assert len(__LMatrix) == __Idwx.size, "Incorrect size of elements"
 203             __Produit = numpy.zeros(__LMatrix[0].size)
 204             for i, col in enumerate(__LMatrix):
 205                 __Produit += float(__Idwx[i]) * col
 206             return __Produit
 207         elif __dotTWith is not None:
 208             _Idwy = numpy.ravel( __dotTWith ).T
 209             assert __LMatrix[0].size == _Idwy.size, "Incorrect size of elements"
 210             __Produit = numpy.zeros(len(__LMatrix))
 211             for i, col in enumerate(__LMatrix):
 212                 __Produit[i] = vfloat( _Idwy @ col)
 213             return __Produit
 214         else:
 215             __Produit = None
 216         return __Produit
 217
 218     # ---------------------------------------------------------
 219     def DirectOperator(self, X, **extraArgs ):
 220         """
 221         Calcul du direct à l'aide de la fonction fournie.
 222
 223         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 224         ne doivent pas être données ici à la fonction utilisateur.
 225         """
 226         logging.debug("FDA Calcul DirectOperator (explicite)")
 227         if self.__mfEnabled:
 228             _HX = self.__userFunction( X, argsAsSerie = True )
 229         else:
 230             _HX = numpy.ravel(self.__userFunction( numpy.ravel(X) ))
 231         #
 232         return _HX
 233
 234     # ---------------------------------------------------------
 235     def TangentMatrix(self, X, dotWith = None, dotTWith = None ):
 236         """
 237         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 238         c'est-à-dire le gradient de H en X. On utilise des différences finies
 239         directionnelles autour du point X. X est un numpy.ndarray.
 240
 241         Différences finies centrées (approximation d'ordre 2):
 242         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 243            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 244            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 245            H( X_moins_dXi )
 246         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 247            le pas 2*dXi
 248         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 249
 250         Différences finies non centrées (approximation d'ordre 1):
 251         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 252            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 253            HX_plus_dXi = H( X_plus_dXi )
 254         2/ On calcule la valeur centrale HX = H(X)
 255         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 256            le pas dXi
 257         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 258
 259         """
 260         logging.debug("FDA Début du calcul de la Jacobienne")
 261         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 262         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 263         #
 264         if X is None or len(X) == 0:
 265             raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
 266         #
 267         _X = numpy.ravel( X )
 268         #
 269         if self.__dX is None:
 270             _dX  = self.__increment * _X
 271         else:
 272             _dX = numpy.ravel( self.__dX )
 273         assert len(_X) == len(_dX), "Inconsistent dX increment length with respect to the X one"
 274         assert _X.size == _dX.size, "Inconsistent dX increment size with respect to the X one"
 275         #
 276         if (_dX == 0.).any():
 277             moyenne = _dX.mean()
 278             if moyenne == 0.:
 279                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 280             else:
 281                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 282         #
 283         __alreadyCalculated  = False
 284         if self.__avoidRC:
 285             __bidon, __alreadyCalculatedP = self.__doublon__( _X, self.__listJPCP, self.__listJPPN, None)
 286             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 287             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 288                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 289                 logging.debug("FDA Cas J déjà calculé, récupération du doublon %i"%__i)
 290         #
 291         if __alreadyCalculated:
 292             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 293             _Jacobienne = self.__listJPCR[__i]
 294             logging.debug("FDA Fin du calcul de la Jacobienne")
 295             if dotWith is not None:
 296                 return numpy.dot(  _Jacobienne, numpy.ravel( dotWith ))
 297             elif dotTWith is not None:
 298                 return numpy.dot(_Jacobienne.T, numpy.ravel( dotTWith ))
 299         else:
 300             logging.debug("FDA   Calcul Jacobienne (explicite)")
 301             if self.__centeredDF:
 302                 #
 303                 if self.__mpEnabled and not self.__mfEnabled:
 304                     funcrepr = {
 305                         "__userFunction__path": self.__userFunction__path,
 306                         "__userFunction__modl": self.__userFunction__modl,
 307                         "__userFunction__name": self.__userFunction__name,
 308                     }
 309                     _jobs = []
 310                     for i in range( len(_dX) ):
 311                         _dXi            = _dX[i]
 312                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 313                         _X_plus_dXi[i]  = _X[i] + _dXi
 314                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 315                         _X_moins_dXi[i] = _X[i] - _dXi
 316                         #
 317                         _jobs.append( ( _X_plus_dXi, self.__extraArgs, funcrepr) )
 318                         _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
 319                     #
 320                     import multiprocessing
 321                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 322                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 323                     self.__pool.close()
 324                     self.__pool.join()
 325                     #
 326                     _Jacobienne  = []
 327                     for i in range( len(_dX) ):
 328                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2 * i] - _HX_plusmoins_dX[2 * i + 1] ) / (2. * _dX[i]) )
 329                     #
 330                 elif self.__mfEnabled:
 331                     _xserie = []
 332                     for i in range( len(_dX) ):
 333                         _dXi            = _dX[i]
 334                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 335                         _X_plus_dXi[i]  = _X[i] + _dXi
 336                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 337                         _X_moins_dXi[i] = _X[i] - _dXi
 338                         #
 339                         _xserie.append( _X_plus_dXi )
 340                         _xserie.append( _X_moins_dXi )
 341                     #
 342                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 343                     #
 344                     _Jacobienne  = []
 345                     for i in range( len(_dX) ):
 346                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2 * i] - _HX_plusmoins_dX[2 * i + 1] ) / (2. * _dX[i]) )
 347                     #
 348                 else:
 349                     _Jacobienne  = []
 350                     for i in range( _dX.size ):
 351                         _dXi            = _dX[i]
 352                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 353                         _X_plus_dXi[i]  = _X[i] + _dXi
 354                         _X_moins_dXi    = numpy.array( _X, dtype=float )
 355                         _X_moins_dXi[i] = _X[i] - _dXi
 356                         #
 357                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 358                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 359                         #
 360                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2. * _dXi) )
 361                 #
 362             else:
 363                 #
 364                 if self.__mpEnabled and not self.__mfEnabled:
 365                     funcrepr = {
 366                         "__userFunction__path": self.__userFunction__path,
 367                         "__userFunction__modl": self.__userFunction__modl,
 368                         "__userFunction__name": self.__userFunction__name,
 369                     }
 370                     _jobs = []
 371                     _jobs.append( (_X, self.__extraArgs, funcrepr) )
 372                     for i in range( len(_dX) ):
 373                         _X_plus_dXi    = numpy.array( _X, dtype=float )
 374                         _X_plus_dXi[i] = _X[i] + _dX[i]
 375                         #
 376                         _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
 377                     #
 378                     import multiprocessing
 379                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 380                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 381                     self.__pool.close()
 382                     self.__pool.join()
 383                     #
 384                     _HX = _HX_plus_dX.pop(0)
 385                     #
 386                     _Jacobienne = []
 387                     for i in range( len(_dX) ):
 388                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 389                     #
 390                 elif self.__mfEnabled:
 391                     _xserie = []
 392                     _xserie.append( _X )
 393                     for i in range( len(_dX) ):
 394                         _X_plus_dXi    = numpy.array( _X, dtype=float )
 395                         _X_plus_dXi[i] = _X[i] + _dX[i]
 396                         #
 397                         _xserie.append( _X_plus_dXi )
 398                     #
 399                     _HX_plus_dX = self.DirectOperator( _xserie )
 400                     #
 401                     _HX = _HX_plus_dX.pop(0)
 402                     #
 403                     _Jacobienne = []
 404                     for i in range( len(_dX) ):
 405                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 406                     #
 407                 else:
 408                     _Jacobienne  = []
 409                     _HX = self.DirectOperator( _X )
 410                     for i in range( _dX.size ):
 411                         _dXi            = _dX[i]
 412                         _X_plus_dXi     = numpy.array( _X, dtype=float )
 413                         _X_plus_dXi[i]  = _X[i] + _dXi
 414                         #
 415                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 416                         #
 417                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 418             #
 419             if (dotWith is not None) or (dotTWith is not None):
 420                 __Produit = self.__listdotwith__(_Jacobienne, dotWith, dotTWith)
 421             else:
 422                 __Produit = None
 423             if __Produit is None or self.__avoidRC:
 424                 _Jacobienne = numpy.transpose( numpy.vstack( _Jacobienne ) )
 425                 if self.__avoidRC:
 426                     if self.__lengthRJ < 0:
 427                         self.__lengthRJ = 2 * _X.size
 428                     while len(self.__listJPCP) > self.__lengthRJ:
 429                         self.__listJPCP.pop(0)
 430                         self.__listJPCI.pop(0)
 431                         self.__listJPCR.pop(0)
 432                         self.__listJPPN.pop(0)
 433                         self.__listJPIN.pop(0)
 434                     self.__listJPCP.append( copy.copy(_X) )
 435                     self.__listJPCI.append( copy.copy(_dX) )
 436                     self.__listJPCR.append( copy.copy(_Jacobienne) )
 437                     self.__listJPPN.append( numpy.linalg.norm(_X) )
 438                     self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 439             logging.debug("FDA Fin du calcul de la Jacobienne")
 440             if __Produit is not None:
 441                 return __Produit
 442         #
 443         return _Jacobienne
 444
 445     # ---------------------------------------------------------
 446     def TangentOperator(self, paire, **extraArgs ):
 447         """
 448         Calcul du tangent à l'aide de la Jacobienne.
 449
 450         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 451         ne doivent pas être données ici à la fonction utilisateur.
 452         """
 453         if self.__mfEnabled:
 454             assert len(paire) == 1, "Incorrect length of arguments"
 455             _paire = paire[0]
 456             assert len(_paire) == 2, "Incorrect number of arguments"
 457         else:
 458             assert len(paire) == 2, "Incorrect number of arguments"
 459             _paire = paire
 460         X, dX = _paire
 461         if dX is None or len(dX) == 0:
 462             #
 463             # Calcul de la forme matricielle si le second argument est None
 464             # -------------------------------------------------------------
 465             _Jacobienne = self.TangentMatrix( X )
 466             if self.__mfEnabled:
 467                 return [_Jacobienne,]
 468             else:
 469                 return _Jacobienne
 470         else:
 471             #
 472             # Calcul de la valeur linéarisée de H en X appliqué à dX
 473             # ------------------------------------------------------
 474             _HtX = self.TangentMatrix( X, dotWith = dX )
 475             if self.__mfEnabled:
 476                 return [_HtX,]
 477             else:
 478                 return _HtX
 479
 480     # ---------------------------------------------------------
 481     def AdjointOperator(self, paire, **extraArgs ):
 482         """
 483         Calcul de l'adjoint à l'aide de la Jacobienne.
 484
 485         NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
 486         ne doivent pas être données ici à la fonction utilisateur.
 487         """
 488         if self.__mfEnabled:
 489             assert len(paire) == 1, "Incorrect length of arguments"
 490             _paire = paire[0]
 491             assert len(_paire) == 2, "Incorrect number of arguments"
 492         else:
 493             assert len(paire) == 2, "Incorrect number of arguments"
 494             _paire = paire
 495         X, Y = _paire
 496         if Y is None or len(Y) == 0:
 497             #
 498             # Calcul de la forme matricielle si le second argument est None
 499             # -------------------------------------------------------------
 500             _JacobienneT = self.TangentMatrix( X ).T
 501             if self.__mfEnabled:
 502                 return [_JacobienneT,]
 503             else:
 504                 return _JacobienneT
 505         else:
 506             #
 507             # Calcul de la valeur de l'adjoint en X appliqué à Y
 508             # --------------------------------------------------
 509             _HaY = self.TangentMatrix( X, dotTWith = Y )
 510             if self.__mfEnabled:
 511                 return [_HaY,]
 512             else:
 513                 return _HaY
 514
 515 # ==============================================================================
 516 def SetInitialDirection( __Direction = [], __Amplitude = 1., __Position = None ):
 517     "Établit ou élabore une direction avec une amplitude"
 518     #
 519     if len(__Direction) == 0 and __Position is None:
 520         raise ValueError("If initial direction is void, current position has to be given")
 521     if abs(float(__Amplitude)) < mpr:
 522         raise ValueError("Amplitude of perturbation can not be zero")
 523     #
 524     if len(__Direction) > 0:
 525         __dX0 = numpy.asarray(__Direction)
 526     else:
 527         __dX0 = []
 528         __X0 = numpy.ravel(numpy.asarray(__Position))
 529         __mX0 = numpy.mean( __X0, dtype=mfp )
 530         if abs(__mX0) < 2 * mpr:
 531             __mX0 = 1.  # Évite le problème de position nulle
 532         for v in __X0:
 533             if abs(v) > 1.e-8:
 534                 __dX0.append( numpy.random.normal(0., abs(v)) )
 535             else:
 536                 __dX0.append( numpy.random.normal(0., __mX0) )
 537     #
 538     __dX0 = numpy.asarray(__dX0, float)  # Évite le problème d'array de taille 1
 539     __dX0 = numpy.ravel( __dX0 )         # Redresse les vecteurs
 540     __dX0 = float(__Amplitude) * __dX0
 541     #
 542     return __dX0
 543
 544 # ==============================================================================
 545 def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
 546     "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
 547     #
 548     __bgCenter = numpy.ravel(__bgCenter)[:, None]
 549     if __nbMembers < 1:
 550         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
 551     #
 552     if __bgCovariance is None:
 553         _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 554     else:
 555         _Z = numpy.random.multivariate_normal(numpy.zeros(__bgCenter.size), __bgCovariance, size=__nbMembers).T
 556         _Perturbations = numpy.tile( __bgCenter, __nbMembers) + _Z
 557     #
 558     return _Perturbations
 559
 560 # ==============================================================================
 561 def EnsembleOfBackgroundPerturbations(
 562         __bgCenter,
 563         __bgCovariance,
 564         __nbMembers,
 565         __withSVD = True ):
 566     "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
 567     def __CenteredRandomAnomalies(Zr, N):
 568         """
 569         Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
 570         notes manuscrites de MB et conforme au code de PS avec eps = -1
 571         """
 572         eps = -1
 573         Q = numpy.identity(N - 1) - numpy.ones((N - 1, N - 1)) / numpy.sqrt(N) / (numpy.sqrt(N) - eps)
 574         Q = numpy.concatenate((Q, [eps * numpy.ones(N - 1) / numpy.sqrt(N)]), axis=0)
 575         R, _ = numpy.linalg.qr(numpy.random.normal(size = (N - 1, N - 1)))
 576         Q = numpy.dot(Q, R)
 577         Zr = numpy.dot(Q, Zr)
 578         return Zr.T
 579     #
 580     __bgCenter = numpy.ravel(__bgCenter).reshape((-1, 1))
 581     if __nbMembers < 1:
 582         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
 583     if __bgCovariance is None:
 584         _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 585     else:
 586         if __withSVD:
 587             _U, _s, _V = numpy.linalg.svd(__bgCovariance, full_matrices=False)
 588             _nbctl = __bgCenter.size
 589             if __nbMembers > _nbctl:
 590                 _Z = numpy.concatenate((numpy.dot(
 591                     numpy.diag(numpy.sqrt(_s[:_nbctl])), _V[:_nbctl]),
 592                     numpy.random.multivariate_normal(numpy.zeros(_nbctl), __bgCovariance, __nbMembers - 1 - _nbctl)), axis = 0)
 593             else:
 594                 _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:__nbMembers - 1])), _V[:__nbMembers - 1])
 595             _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
 596             _Perturbations = __bgCenter + _Zca
 597         else:
 598             if max(abs(__bgCovariance.flatten())) > 0:
 599                 _nbctl = __bgCenter.size
 600                 _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl), __bgCovariance, __nbMembers - 1)
 601                 _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
 602                 _Perturbations = __bgCenter + _Zca
 603             else:
 604                 _Perturbations = numpy.tile( __bgCenter, __nbMembers)
 605     #
 606     return _Perturbations
 607
 608 # ==============================================================================
 609 def EnsembleMean( __Ensemble ):
 610     "Renvoie la moyenne empirique d'un ensemble"
 611     return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1, 1))
 612
 613 # ==============================================================================
 614 def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1. ):
 615     "Renvoie les anomalies centrées à partir d'un ensemble"
 616     if __OptMean is None:
 617         __Em = EnsembleMean( __Ensemble )
 618     else:
 619         __Em = numpy.ravel( __OptMean ).reshape((-1, 1))
 620     #
 621     return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
 622
 623 # ==============================================================================
 624 def EnsembleErrorCovariance( __Ensemble, __Quick = False ):
 625     "Renvoie l'estimation empirique de la covariance d'ensemble"
 626     if __Quick:
 627         # Covariance rapide mais rarement définie positive
 628         __Covariance = numpy.cov( __Ensemble )
 629     else:
 630         # Résultat souvent identique à numpy.cov, mais plus robuste
 631         __n, __m = numpy.asarray( __Ensemble ).shape
 632         __Anomalies = EnsembleOfAnomalies( __Ensemble )
 633         # Estimation empirique
 634         __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m - 1)
 635         # Assure la symétrie
 636         __Covariance = ( __Covariance + __Covariance.T ) * 0.5
 637         # Assure la positivité
 638         __epsilon    = mpr * numpy.trace( __Covariance )
 639         __Covariance = __Covariance + __epsilon * numpy.identity(__n)
 640     #
 641     return __Covariance
 642
 643 # ==============================================================================
 644 def SingularValuesEstimation( __Ensemble, __Using = "SVDVALS"):
 645     "Renvoie les valeurs singulières de l'ensemble et leur carré"
 646     if __Using == "SVDVALS":  # Recommandé
 647         __sv   = scipy.linalg.svdvals( __Ensemble )
 648         __svsq = __sv**2
 649     elif __Using == "SVD":
 650         _, __sv, _ = numpy.linalg.svd( __Ensemble )
 651         __svsq = __sv**2
 652     elif __Using == "EIG":  # Lent
 653         __eva, __eve = numpy.linalg.eig( __Ensemble @ __Ensemble.T )
 654         __svsq = numpy.sort(numpy.abs(numpy.real( __eva )))[::-1]
 655         __sv   = numpy.sqrt( __svsq )
 656     elif __Using == "EIGH":
 657         __eva, __eve = numpy.linalg.eigh( __Ensemble @ __Ensemble.T )
 658         __svsq = numpy.sort(numpy.abs(numpy.real( __eva )))[::-1]
 659         __sv   = numpy.sqrt( __svsq )
 660     elif __Using == "EIGVALS":
 661         __eva  = numpy.linalg.eigvals( __Ensemble @ __Ensemble.T )
 662         __svsq = numpy.sort(numpy.abs(numpy.real( __eva )))[::-1]
 663         __sv   = numpy.sqrt( __svsq )
 664     elif __Using == "EIGVALSH":
 665         __eva = numpy.linalg.eigvalsh( __Ensemble @ __Ensemble.T )
 666         __svsq = numpy.sort(numpy.abs(numpy.real( __eva )))[::-1]
 667         __sv   = numpy.sqrt( __svsq )
 668     else:
 669         raise ValueError("Error in requested variant name: %s"%__Using)
 670     #
 671     __tisv = __svsq / __svsq.sum()
 672     __qisv = 1. - __svsq.cumsum() / __svsq.sum()
 673     # Différence à 1.e-16 : __qisv = 1. - __tisv.cumsum()
 674     #
 675     return __sv, __svsq, __tisv, __qisv
 676
 677 # ==============================================================================
 678 def MaxL2NormByColumn(__Ensemble, __LcCsts = False, __IncludedPoints = []):
 679     "Maximum des normes L2 calculées par colonne"
 680     if __LcCsts and len(__IncludedPoints) > 0:
 681         normes = numpy.linalg.norm(
 682             numpy.take(__Ensemble, __IncludedPoints, axis=0, mode='clip'),
 683             axis = 0,
 684         )
 685     else:
 686         normes = numpy.linalg.norm( __Ensemble, axis = 0)
 687     nmax = numpy.max(normes)
 688     imax = numpy.argmax(normes)
 689     return nmax, imax, normes
 690
 691 def MaxLinfNormByColumn(__Ensemble, __LcCsts = False, __IncludedPoints = []):
 692     "Maximum des normes Linf calculées par colonne"
 693     if __LcCsts and len(__IncludedPoints) > 0:
 694         normes = numpy.linalg.norm(
 695             numpy.take(__Ensemble, __IncludedPoints, axis=0, mode='clip'),
 696             axis = 0, ord=numpy.inf,
 697         )
 698     else:
 699         normes = numpy.linalg.norm( __Ensemble, axis = 0, ord=numpy.inf)
 700     nmax = numpy.max(normes)
 701     imax = numpy.argmax(normes)
 702     return nmax, imax, normes
 703
 704 def InterpolationErrorByColumn(
 705         __Ensemble = None, __Basis = None, __Points = None, __M = 2,  # Usage 1
 706         __Differences = None,                                         # Usage 2
 707         __ErrorNorm = None,                                           # Commun
 708         __LcCsts = False, __IncludedPoints = [],                      # Commun
 709         __CDM = False,  # ComputeMaxDifference                        # Commun
 710         __RMU = False,  # ReduceMemoryUse                             # Commun
 711         __FTL = False,  # ForceTril                                   # Commun
 712         ):   # noqa: E123
 713     "Analyse des normes d'erreurs d'interpolation calculées par colonne"
 714     if __ErrorNorm == "L2":
 715         NormByColumn = MaxL2NormByColumn
 716     else:
 717         NormByColumn = MaxLinfNormByColumn
 718     #
 719     if __Differences is None and not __RMU:  # Usage 1
 720         if __FTL:
 721             rBasis = numpy.tril( __Basis[__Points, :] )
 722         else:
 723             rBasis = __Basis[__Points, :]
 724         rEnsemble = __Ensemble[__Points, :]
 725         #
 726         if __M > 1:
 727             rBasis_inv = numpy.linalg.inv(rBasis)
 728             Interpolator = numpy.dot(__Basis, numpy.dot(rBasis_inv, rEnsemble))
 729         else:
 730             rBasis_inv = 1. / rBasis
 731             Interpolator = numpy.outer(__Basis, numpy.outer(rBasis_inv, rEnsemble))
 732         #
 733         differences = __Ensemble - Interpolator
 734         #
 735         error, nbr, _ = NormByColumn(differences, __LcCsts, __IncludedPoints)
 736         #
 737         if __CDM:
 738             maxDifference = differences[:, nbr]
 739         #
 740     elif __Differences is None and __RMU:  # Usage 1
 741         if __FTL:
 742             rBasis = numpy.tril( __Basis[__Points, :] )
 743         else:
 744             rBasis = __Basis[__Points, :]
 745         rEnsemble = __Ensemble[__Points, :]
 746         #
 747         if __M > 1:
 748             rBasis_inv = numpy.linalg.inv(rBasis)
 749             rCoordinates = numpy.dot(rBasis_inv, rEnsemble)
 750         else:
 751             rBasis_inv = 1. / rBasis
 752             rCoordinates = numpy.outer(rBasis_inv, rEnsemble)
 753         #
 754         error = 0.
 755         nbr = -1
 756         for iCol in range(__Ensemble.shape[1]):
 757             if __M > 1:
 758                 iDifference = __Ensemble[:, iCol] - numpy.dot(__Basis, rCoordinates[:, iCol])
 759             else:
 760                 iDifference = __Ensemble[:, iCol] - numpy.ravel(numpy.outer(__Basis, rCoordinates[:, iCol]))
 761             #
 762             normDifference, _, _ = NormByColumn(iDifference, __LcCsts, __IncludedPoints)
 763             #
 764             if normDifference > error:
 765                 error         = normDifference
 766                 nbr           = iCol
 767         #
 768         if __CDM:
 769             maxDifference = __Ensemble[:, nbr] - numpy.dot(__Basis, rCoordinates[:, nbr])
 770         #
 771     else:  # Usage 2
 772         differences = __Differences
 773         #
 774         error, nbr, _ = NormByColumn(differences, __LcCsts, __IncludedPoints)
 775         #
 776         if __CDM:
 777             # faire cette variable intermédiaire coûte cher
 778             maxDifference = differences[:, nbr]
 779     #
 780     if __CDM:
 781         return error, nbr, maxDifference
 782     else:
 783         return error, nbr
 784
 785 # ==============================================================================
 786 def EnsemblePerturbationWithGivenCovariance(
 787         __Ensemble,
 788         __Covariance,
 789         __Seed = None ):
 790     "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
 791     if hasattr(__Covariance, "assparsematrix"):
 792         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix()) / abs(__Ensemble).mean() < mpr).all():
 793             # Traitement d'une covariance nulle ou presque
 794             return __Ensemble
 795         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
 796             # Traitement d'une covariance nulle ou presque
 797             return __Ensemble
 798     else:
 799         if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance) / abs(__Ensemble).mean() < mpr).all():
 800             # Traitement d'une covariance nulle ou presque
 801             return __Ensemble
 802         if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
 803             # Traitement d'une covariance nulle ou presque
 804             return __Ensemble
 805     #
 806     __n, __m = __Ensemble.shape
 807     if __Seed is not None:
 808         numpy.random.seed(__Seed)
 809     #
 810     if hasattr(__Covariance, "isscalar") and __Covariance.isscalar():
 811         # Traitement d'une covariance multiple de l'identité
 812         __zero = 0.
 813         __std  = numpy.sqrt(__Covariance.assparsematrix())
 814         __Ensemble += numpy.random.normal(__zero, __std, size=(__m, __n)).T
 815     #
 816     elif hasattr(__Covariance, "isvector") and __Covariance.isvector():
 817         # Traitement d'une covariance diagonale avec variances non identiques
 818         __zero = numpy.zeros(__n)
 819         __std  = numpy.sqrt(__Covariance.assparsematrix())
 820         __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
 821     #
 822     elif hasattr(__Covariance, "ismatrix") and __Covariance.ismatrix():
 823         # Traitement d'une covariance pleine
 824         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
 825     #
 826     elif isinstance(__Covariance, numpy.ndarray):
 827         # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
 828         __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
 829     #
 830     else:
 831         raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
 832     #
 833     return __Ensemble
 834
 835 # ==============================================================================
 836 def CovarianceInflation(
 837         __InputCovOrEns,
 838         __InflationType   = None,
 839         __InflationFactor = None,
 840         __BackgroundCov   = None ):
 841     """
 842     Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
 843
 844     Synthèse : Hunt 2007, section 2.3.5
 845     """
 846     if __InflationFactor is None:
 847         return __InputCovOrEns
 848     else:
 849         __InflationFactor = float(__InflationFactor)
 850     #
 851     __InputCovOrEns = numpy.asarray(__InputCovOrEns)
 852     if __InputCovOrEns.size == 0:
 853         return __InputCovOrEns
 854     #
 855     if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
 856         if __InflationFactor < 1.:
 857             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 858         if __InflationFactor < 1. + mpr:  # No inflation = 1
 859             return __InputCovOrEns
 860         __OutputCovOrEns = __InflationFactor**2 * __InputCovOrEns
 861     #
 862     elif __InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
 863         if __InflationFactor < 1.:
 864             raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
 865         if __InflationFactor < 1. + mpr:  # No inflation = 1
 866             return __InputCovOrEns
 867         __InputCovOrEnsMean = __InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
 868         __OutputCovOrEns = __InputCovOrEnsMean[:, numpy.newaxis] \
 869             + __InflationFactor * (__InputCovOrEns - __InputCovOrEnsMean[:, numpy.newaxis])
 870     #
 871     elif __InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
 872         if __InflationFactor < 0.:
 873             raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
 874         if __InflationFactor < mpr:  # No inflation = 0
 875             return __InputCovOrEns
 876         __n, __m = __InputCovOrEns.shape
 877         if __n != __m:
 878             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 879         __tr = __InputCovOrEns.trace() / __n
 880         if __InflationFactor > __tr:
 881             raise ValueError("Inflation factor for additive inflation has to be small over %.0e."%__tr)
 882         __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * numpy.identity(__n)
 883     #
 884     elif __InflationType == "HybridOnBackgroundCovariance":
 885         if __InflationFactor < 0.:
 886             raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
 887         if __InflationFactor < mpr:  # No inflation = 0
 888             return __InputCovOrEns
 889         __n, __m = __InputCovOrEns.shape
 890         if __n != __m:
 891             raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
 892         if __BackgroundCov is None:
 893             raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
 894         if __InputCovOrEns.shape != __BackgroundCov.shape:
 895             raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
 896         __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov
 897     #
 898     elif __InflationType == "Relaxation":
 899         raise NotImplementedError("Relaxation inflation type not implemented")
 900     #
 901     else:
 902         raise ValueError("Error in inflation type, '%s' is not a valid keyword."%__InflationType)
 903     #
 904     return __OutputCovOrEns
 905
 906 # ==============================================================================
 907 def HessienneEstimation( __selfA, __nb, __HaM, __HtM, __BI, __RI ):
 908     "Estimation de la Hessienne"
 909     #
 910     __HessienneI = []
 911     for i in range(int(__nb)):
 912         __ee    = numpy.zeros((__nb, 1))
 913         __ee[i] = 1.
 914         __HtEE  = numpy.dot(__HtM, __ee).reshape((-1, 1))
 915         __HessienneI.append( numpy.ravel( __BI * __ee + __HaM * (__RI * __HtEE) ) )
 916     #
 917     __A = numpy.linalg.inv(numpy.array( __HessienneI ))
 918     __A = (__A + __A.T) * 0.5  # Symétrie
 919     __A = __A + mpr * numpy.trace( __A ) * numpy.identity(__nb)  # Positivité
 920     #
 921     if min(__A.shape) != max(__A.shape):
 922         raise ValueError(
 923             "The %s a posteriori covariance matrix A"%(__selfA._name,) + \
 924             " is of shape %s, despites it has to be a"%(str(__A.shape),) + \
 925             " squared matrix. There is an error in the observation operator," + \
 926             " please check it.")
 927     if (numpy.diag(__A) < 0).any():
 928         raise ValueError(
 929             "The %s a posteriori covariance matrix A"%(__selfA._name,) + \
 930             " has at least one negative value on its diagonal. There is an" + \
 931             " error in the observation operator, please check it.")
 932     if logging.getLogger().level < logging.WARNING:  # La vérification n'a lieu qu'en debug
 933         try:
 934             numpy.linalg.cholesky( __A )
 935             logging.debug("%s La matrice de covariance a posteriori A est bien symétrique définie positive."%(__selfA._name,))
 936         except Exception:
 937             raise ValueError(
 938                 "The %s a posteriori covariance matrix A"%(__selfA._name,) + \
 939                 " is not symmetric positive-definite. Please check your a" + \
 940                 " priori covariances and your observation operator.")
 941     #
 942     return __A
 943
 944 # ==============================================================================
 945 def QuantilesEstimations( selfA, A, Xa, HXa = None, Hm = None, HtM = None ):
 946     "Estimation des quantiles a posteriori à partir de A>0 (selfA est modifié)"
 947     nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
 948     #
 949     # Traitement des bornes
 950     if "StateBoundsForQuantiles" in selfA._parameters:
 951         LBounds = selfA._parameters["StateBoundsForQuantiles"]  # Prioritaire
 952     elif "Bounds" in selfA._parameters:
 953         LBounds = selfA._parameters["Bounds"]  # Défaut raisonnable
 954     else:
 955         LBounds = None
 956     if LBounds is not None:
 957         LBounds = ForceNumericBounds( LBounds )
 958     __Xa = numpy.ravel(Xa)
 959     #
 960     # Échantillonnage des états
 961     YfQ  = None
 962     EXr  = None
 963     for i in range(nbsamples):
 964         if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None and HXa is not None:
 965             dXr = (numpy.random.multivariate_normal(__Xa, A) - __Xa).reshape((-1, 1))
 966             if LBounds is not None:  # "EstimateProjection" par défaut
 967                 dXr = numpy.max(numpy.hstack((dXr, LBounds[:, 0].reshape((-1, 1))) - __Xa.reshape((-1, 1))), axis=1)
 968                 dXr = numpy.min(numpy.hstack((dXr, LBounds[:, 1].reshape((-1, 1))) - __Xa.reshape((-1, 1))), axis=1)
 969             dYr = HtM @ dXr
 970             Yr = HXa.reshape((-1, 1)) + dYr
 971             if selfA._toStore("SampledStateForQuantiles"):
 972                 Xr = __Xa + numpy.ravel(dXr)
 973         elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
 974             Xr = numpy.random.multivariate_normal(__Xa, A)
 975             if LBounds is not None:  # "EstimateProjection" par défaut
 976                 Xr = numpy.max(numpy.hstack((Xr.reshape((-1, 1)), LBounds[:, 0].reshape((-1, 1)))), axis=1)
 977                 Xr = numpy.min(numpy.hstack((Xr.reshape((-1, 1)), LBounds[:, 1].reshape((-1, 1)))), axis=1)
 978             Yr = numpy.asarray(Hm( Xr ))
 979         else:
 980             raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
 981         #
 982         if YfQ is None:
 983             YfQ = Yr.reshape((-1, 1))
 984             if selfA._toStore("SampledStateForQuantiles"):
 985                 EXr = Xr.reshape((-1, 1))
 986         else:
 987             YfQ = numpy.hstack((YfQ, Yr.reshape((-1, 1))))
 988             if selfA._toStore("SampledStateForQuantiles"):
 989                 EXr = numpy.hstack((EXr, Xr.reshape((-1, 1))))
 990     #
 991     # Extraction des quantiles
 992     YfQ.sort(axis=-1)
 993     YQ = None
 994     for quantile in selfA._parameters["Quantiles"]:
 995         if not (0. <= float(quantile) <= 1.):
 996             continue
 997         indice = int(nbsamples * float(quantile) - 1. / nbsamples)
 998         if YQ is None:
 999             YQ = YfQ[:, indice].reshape((-1, 1))
1000         else:
1001             YQ = numpy.hstack((YQ, YfQ[:, indice].reshape((-1, 1))))
1002     if YQ is not None:  # Liste non vide de quantiles
1003         selfA.StoredVariables["SimulationQuantiles"].store( YQ )
1004     if selfA._toStore("SampledStateForQuantiles"):
1005         selfA.StoredVariables["SampledStateForQuantiles"].store( EXr )
1006     #
1007     return 0
1008
1009 # ==============================================================================
1010 def ForceNumericBounds( __Bounds, __infNumbers = True ):
1011     "Force les bornes à être des valeurs numériques, sauf si globalement None"
1012     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
1013     if __Bounds is None:
1014         return None
1015     #
1016     # Converti toutes les bornes individuelles None à +/- l'infini chiffré
1017     __Bounds = numpy.asarray( __Bounds, dtype=float ).reshape((-1, 2))
1018     if len(__Bounds.shape) != 2 or __Bounds.shape[0] == 0 or __Bounds.shape[1] != 2:
1019         raise ValueError("Incorrectly shaped bounds data (effective shape is %s)"%(__Bounds.shape,))
1020     if __infNumbers:
1021         __Bounds[numpy.isnan(__Bounds[:, 0]), 0] = -float('inf')
1022         __Bounds[numpy.isnan(__Bounds[:, 1]), 1] = float('inf')
1023     else:
1024         __Bounds[numpy.isnan(__Bounds[:, 0]), 0] = -sys.float_info.max
1025         __Bounds[numpy.isnan(__Bounds[:, 1]), 1] = sys.float_info.max
1026     return __Bounds
1027
1028 # ==============================================================================
1029 def RecentredBounds( __Bounds, __Center, __Scale = None ):
1030     "Recentre les bornes autour de 0, sauf si globalement None"
1031     # Conserve une valeur par défaut à None s'il n'y a pas de bornes
1032     if __Bounds is None:
1033         return None
1034     #
1035     if __Scale is None:
1036         # Recentre les valeurs numériques de bornes
1037         return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1, 1))
1038     else:
1039         # Recentre les valeurs numériques de bornes et change l'échelle par une matrice
1040         return __Scale @ (ForceNumericBounds( __Bounds, False ) - numpy.ravel( __Center ).reshape((-1, 1)))
1041
1042 # ==============================================================================
1043 def ApplyBounds( __Vector, __Bounds, __newClip = True ):
1044     "Applique des bornes numériques à un point"
1045     # Conserve une valeur par défaut s'il n'y a pas de bornes
1046     if __Bounds is None:
1047         return __Vector
1048     #
1049     if not isinstance(__Vector, numpy.ndarray):  # Is an array
1050         raise ValueError("Incorrect array definition of vector data")
1051     if not isinstance(__Bounds, numpy.ndarray):  # Is an array
1052         raise ValueError("Incorrect array definition of bounds data")
1053     if 2 * __Vector.size != __Bounds.size:  # Is a 2 column array of vector length
1054         raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size, __Vector.size))
1055     if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
1056         raise ValueError("Incorrectly shaped bounds data")
1057     #
1058     if __newClip:
1059         __Vector = __Vector.clip(
1060             __Bounds[:, 0].reshape(__Vector.shape),
1061             __Bounds[:, 1].reshape(__Vector.shape),
1062         )
1063     else:
1064         __Vector = numpy.max(numpy.hstack((__Vector.reshape((-1, 1)), numpy.asmatrix(__Bounds)[:, 0])), axis=1)
1065         __Vector = numpy.min(numpy.hstack((__Vector.reshape((-1, 1)), numpy.asmatrix(__Bounds)[:, 1])), axis=1)
1066         __Vector = numpy.asarray(__Vector)
1067     #
1068     return __Vector
1069
1070 # ==============================================================================
1071 def VariablesAndIncrementsBounds( __Bounds, __BoxBounds, __Xini, __Name, __Multiplier = 1. ):
1072     __Bounds    = ForceNumericBounds( __Bounds )
1073     __BoxBounds = ForceNumericBounds( __BoxBounds )
1074     if __Bounds is None and __BoxBounds is None:
1075         raise ValueError("Algorithm %s requires bounds on all variables (by Bounds), or on all variable increments (by BoxBounds), or both, to be explicitly given."%(__Name,))
1076     elif __Bounds is None and __BoxBounds is not None:
1077         __Bounds    = __BoxBounds
1078         logging.debug("%s Definition of parameter bounds from current parameter increment bounds"%(__Name,))
1079     elif __Bounds is not None and __BoxBounds is None:
1080         __BoxBounds = __Multiplier * (__Bounds - __Xini.reshape((-1, 1)))  # "M * [Xmin,Xmax]-Xini"
1081         logging.debug("%s Definition of parameter increment bounds from current parameter bounds"%(__Name,))
1082     return __Bounds, __BoxBounds
1083
1084 # ==============================================================================
1085 def Apply3DVarRecentringOnEnsemble( __EnXn, __EnXf, __Ynpu, __HO, __R, __B, __SuppPars ):
1086     "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
1087     __Betaf = __SuppPars["HybridCovarianceEquilibrium"]
1088     #
1089     Xf = EnsembleMean( __EnXf )
1090     Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
1091     Pf = (1 - __Betaf) * __B.asfullmatrix(Xf.size) + __Betaf * Pf
1092     #
1093     selfB = PartialAlgorithm("3DVAR")
1094     selfB._parameters["Minimizer"] = "LBFGSB"
1095     selfB._parameters["MaximumNumberOfIterations"] = __SuppPars["HybridMaximumNumberOfIterations"]
1096     selfB._parameters["CostDecrementTolerance"] = __SuppPars["HybridCostDecrementTolerance"]
1097     selfB._parameters["ProjectedGradientTolerance"] = -1
1098     selfB._parameters["GradientNormTolerance"] = 1.e-05
1099     selfB._parameters["StoreInternalVariables"] = False
1100     selfB._parameters["optiprint"] = -1
1101     selfB._parameters["optdisp"] = 0
1102     selfB._parameters["Bounds"] = None
1103     selfB._parameters["InitializationPoint"] = Xf
1104     from daAlgorithms.Atoms import std3dvar
1105     std3dvar.std3dvar(selfB, Xf, __Ynpu, None, __HO, None, __R, Pf)
1106     Xa = selfB.get("Analysis")[-1].reshape((-1, 1))
1107     del selfB
1108     #
1109     return Xa + EnsembleOfAnomalies( __EnXn )
1110
1111 # ==============================================================================
1112 def GenerateRandomPointInHyperSphere( __Center, __Radius ):
1113     "Génère un point aléatoire uniformément à l'intérieur d'une hyper-sphère"
1114     __Dimension  = numpy.asarray( __Center ).size
1115     __GaussDelta = numpy.random.normal( 0, 1, size=__Center.shape )
1116     __VectorNorm = numpy.linalg.norm( __GaussDelta )
1117     __PointOnHS  = __Radius * (__GaussDelta / __VectorNorm)
1118     __MoveInHS   = math.exp( math.log(numpy.random.uniform()) / __Dimension)  # rand()**1/n
1119     __PointInHS  = __MoveInHS * __PointOnHS
1120     return __Center + __PointInHS
1121
1122 # ==============================================================================
1123 class GenerateWeightsAndSigmaPoints(object):
1124     "Génère les points sigma et les poids associés"
1125
1126     def __init__(self,
1127                  Nn=0, EO="State", VariantM="UKF",
1128                  Alpha=None, Beta=2., Kappa=0.):
1129         self.Nn = int(Nn)
1130         self.Alpha = numpy.longdouble( Alpha )
1131         self.Beta  = numpy.longdouble( Beta )
1132         if abs(Kappa) < 2 * mpr:
1133             if EO == "Parameters" and VariantM == "UKF":
1134                 self.Kappa = 3 - self.Nn
1135             else:  # EO == "State":
1136                 self.Kappa = 0
1137         else:
1138             self.Kappa = Kappa
1139         self.Kappa  = numpy.longdouble( self.Kappa )
1140         self.Lambda = self.Alpha**2 * ( self.Nn + self.Kappa ) - self.Nn
1141         self.Gamma  = self.Alpha * numpy.sqrt( self.Nn + self.Kappa )
1142         # Rq.: Gamma = sqrt(n+Lambda) = Alpha*sqrt(n+Kappa)
1143         assert 0. < self.Alpha <= 1., "Alpha has to be between 0 strictly and 1 included"
1144         #
1145         if VariantM == "UKF":
1146             self.Wm, self.Wc, self.SC = self.__UKF2000()
1147         elif VariantM == "S3F":
1148             self.Wm, self.Wc, self.SC = self.__S3F2022()
1149         elif VariantM == "MSS":
1150             self.Wm, self.Wc, self.SC = self.__MSS2011()
1151         elif VariantM == "5OS":
1152             self.Wm, self.Wc, self.SC = self.__5OS2002()
1153         else:
1154             raise ValueError("Variant \"%s\" is not a valid one."%VariantM)
1155
1156     def __UKF2000(self):
1157         "Standard Set, Julier et al. 2000 (aka Canonical UKF)"
1158         # Rq.: W^{(m)}_{i=/=0} = 1. / (2.*(n + Lambda))
1159         Winn = 1. / (2. * ( self.Nn + self.Kappa ) * self.Alpha**2)
1160         Ww = []
1161         Ww.append( 0. )
1162         for point in range(2 * self.Nn):
1163             Ww.append( Winn )
1164         # Rq.: LsLpL = Lambda / (n + Lambda)
1165         LsLpL = 1. - self.Nn / (self.Alpha**2 * ( self.Nn + self.Kappa ))
1166         Wm = numpy.array( Ww )
1167         Wm[0] = LsLpL
1168         Wc = numpy.array( Ww )
1169         Wc[0] = LsLpL + (1. - self.Alpha**2 + self.Beta)
1170         # OK: assert abs(Wm.sum()-1.) < self.Nn*mpr, "UKF ill-conditioned %s >= %s"%(abs(Wm.sum()-1.), self.Nn*mpr)
1171         #
1172         SC = numpy.zeros((self.Nn, len(Ww)))
1173         for ligne in range(self.Nn):
1174             it = ligne + 1
1175             SC[ligne, it          ] = self.Gamma
1176             SC[ligne, self.Nn + it] = -self.Gamma
1177         #
1178         return Wm, Wc, SC
1179
1180     def __S3F2022(self):
1181         "Scaled Spherical Simplex Set, Papakonstantinou et al. 2022"
1182         # Rq.: W^{(m)}_{i=/=0} = (n + Kappa) / ((n + Lambda) * (n + 1 + Kappa))
1183         Winn = 1. / ((self.Nn + 1. + self.Kappa) * self.Alpha**2)
1184         Ww = []
1185         Ww.append( 0. )
1186         for point in range(self.Nn + 1):
1187             Ww.append( Winn )
1188         # Rq.: LsLpL = Lambda / (n + Lambda)
1189         LsLpL = 1. - self.Nn / (self.Alpha**2 * ( self.Nn + self.Kappa ))
1190         Wm = numpy.array( Ww )
1191         Wm[0] = LsLpL
1192         Wc = numpy.array( Ww )
1193         Wc[0] = LsLpL + (1. - self.Alpha**2 + self.Beta)
1194         # OK: assert abs(Wm.sum()-1.) < self.Nn*mpr, "S3F ill-conditioned %s >= %s"%(abs(Wm.sum()-1.), self.Nn*mpr)
1195         #
1196         SC = numpy.zeros((self.Nn, len(Ww)))
1197         for ligne in range(self.Nn):
1198             it = ligne + 1
1199             q_t = it / math.sqrt( it * (it + 1) * Winn )
1200             SC[ligne, 1:it + 1] = -q_t / it
1201             SC[ligne, it + 1  ] = q_t
1202         #
1203         return Wm, Wc, SC
1204
1205     def __MSS2011(self):
1206         "Minimum Set, Menegaz et al. 2011"
1207         rho2 = (1 - self.Alpha) / self.Nn
1208         Cc = numpy.real(scipy.linalg.sqrtm( numpy.identity(self.Nn) - rho2 ))
1209         Ww = self.Alpha * rho2 * scipy.linalg.inv(Cc) @ numpy.ones(self.Nn) @ scipy.linalg.inv(Cc.T)
1210         Wm = Wc = numpy.concatenate((Ww, [self.Alpha]))
1211         # OK: assert abs(Wm.sum()-1.) < self.Nn*mpr, "MSS ill-conditioned %s >= %s"%(abs(Wm.sum()-1.), self.Nn*mpr)
1212         #
1213         # inv(sqrt(W)) = diag(inv(sqrt(W)))
1214         SC1an = Cc @ numpy.diag(1. / numpy.sqrt( Ww ))
1215         SCnpu = (- numpy.sqrt(rho2) / numpy.sqrt(self.Alpha)) * numpy.ones(self.Nn).reshape((-1, 1))
1216         SC = numpy.concatenate((SC1an, SCnpu), axis=1)
1217         #
1218         return Wm, Wc, SC
1219
1220     def __5OS2002(self):
1221         "Fifth Order Set, Lerner 2002"
1222         Ww = []
1223         for point in range(2 * self.Nn):
1224             Ww.append( (4. - self.Nn) / 18. )
1225         for point in range(2 * self.Nn, 2 * self.Nn**2):
1226             Ww.append( 1. / 36. )
1227         Ww.append( (self.Nn**2 - 7 * self.Nn) / 18. + 1.)
1228         Wm = Wc = numpy.array( Ww )
1229         # OK: assert abs(Wm.sum()-1.) < self.Nn*mpr, "5OS ill-conditioned %s >= %s"%(abs(Wm.sum()-1.), self.Nn*mpr)
1230         #
1231         xi1n  = numpy.diag( math.sqrt(3) * numpy.ones( self.Nn ) )
1232         xi2n  = numpy.diag( -math.sqrt(3) * numpy.ones( self.Nn ) )
1233         #
1234         xi3n1 = numpy.zeros((int((self.Nn - 1) * self.Nn / 2), self.Nn), dtype=float)
1235         xi3n2 = numpy.zeros((int((self.Nn - 1) * self.Nn / 2), self.Nn), dtype=float)
1236         xi4n1 = numpy.zeros((int((self.Nn - 1) * self.Nn / 2), self.Nn), dtype=float)
1237         xi4n2 = numpy.zeros((int((self.Nn - 1) * self.Nn / 2), self.Nn), dtype=float)
1238         ia = 0
1239         for i1 in range(self.Nn - 1):
1240             for i2 in range(i1 + 1, self.Nn):
1241                 xi3n1[ia, i1] = xi3n2[ia, i2] = math.sqrt(3)
1242                 xi3n2[ia, i1] = xi3n1[ia, i2] = -math.sqrt(3)
1243                 # --------------------------------
1244                 xi4n1[ia, i1] = xi4n1[ia, i2] = math.sqrt(3)
1245                 xi4n2[ia, i1] = xi4n2[ia, i2] = -math.sqrt(3)
1246                 ia += 1
1247         SC = numpy.concatenate((xi1n, xi2n, xi3n1, xi3n2, xi4n1, xi4n2, numpy.zeros((1, self.Nn)))).T
1248         #
1249         return Wm, Wc, SC
1250
1251     def nbOfPoints(self):
1252         assert self.Nn      == self.SC.shape[0], "Size mismatch %i =/= %i"%(self.Nn, self.SC.shape[0])
1253         assert self.Wm.size == self.SC.shape[1], "Size mismatch %i =/= %i"%(self.Wm.size, self.SC.shape[1])
1254         assert self.Wm.size == self.Wc.size, "Size mismatch %i =/= %i"%(self.Wm.size, self.Wc.size)
1255         return self.Wm.size
1256
1257     def get(self):
1258         return self.Wm, self.Wc, self.SC
1259
1260     def __repr__(self):
1261         "x.__repr__() <==> repr(x)"
1262         msg  = ""
1263         msg += "    Alpha   = %s\n"%self.Alpha
1264         msg += "    Beta    = %s\n"%self.Beta
1265         msg += "    Kappa   = %s\n"%self.Kappa
1266         msg += "    Lambda  = %s\n"%self.Lambda
1267         msg += "    Gamma   = %s\n"%self.Gamma
1268         msg += "    Wm      = %s\n"%self.Wm
1269         msg += "    Wc      = %s\n"%self.Wc
1270         msg += "    sum(Wm) = %s\n"%numpy.sum(self.Wm)
1271         msg += "    SC      =\n%s\n"%self.SC
1272         return msg
1273
1274 # ==============================================================================
1275 def GetNeighborhoodTopology( __ntype, __ipop ):
1276     "Renvoi une topologie de connexion pour une population de points"
1277     if __ntype in ["FullyConnectedNeighborhood", "FullyConnectedNeighbourhood", "gbest"]:
1278         __topology = [__ipop for __i in __ipop]
1279     elif __ntype in ["RingNeighborhoodWithRadius1", "RingNeighbourhoodWithRadius1", "lbest"]:
1280         __cpop = list(__ipop[-1:]) + list(__ipop) + list(__ipop[:1])
1281         __topology = [__cpop[__n:__n + 3] for __n in range(len(__ipop))]
1282     elif __ntype in ["RingNeighborhoodWithRadius2", "RingNeighbourhoodWithRadius2"]:
1283         __cpop = list(__ipop[-2:]) + list(__ipop) + list(__ipop[:2])
1284         __topology = [__cpop[__n:__n + 5] for __n in range(len(__ipop))]
1285     elif __ntype in ["AdaptativeRandomWith3Neighbors", "AdaptativeRandomWith3Neighbours", "abest"]:
1286         __cpop = 3 * list(__ipop)
1287         __topology = [[__i] + list(numpy.random.choice(__cpop, 3)) for __i in __ipop]
1288     elif __ntype in ["AdaptativeRandomWith5Neighbors", "AdaptativeRandomWith5Neighbours"]:
1289         __cpop = 5 * list(__ipop)
1290         __topology = [[__i] + list(numpy.random.choice(__cpop, 5)) for __i in __ipop]
1291     else:
1292         raise ValueError("Swarm topology type unavailable because name \"%s\" is unknown."%__ntype)
1293     return __topology
1294
1295 # ==============================================================================
1296 def FindIndexesFromNames( __NameOfLocations = None, __ExcludeLocations = None, ForceArray = False ):
1297     "Exprime les indices des noms exclus, en ignorant les absents"
1298     if __ExcludeLocations is None:
1299         __ExcludeIndexes = ()
1300     elif isinstance(__ExcludeLocations, (list, numpy.ndarray, tuple)) and len(__ExcludeLocations) == 0:
1301         __ExcludeIndexes = ()
1302     # ----------
1303     elif __NameOfLocations is None:
1304         try:
1305             __ExcludeIndexes = numpy.asarray(__ExcludeLocations, dtype=int)
1306         except ValueError as e:
1307             if "invalid literal for int() with base 10:" in str(e):
1308                 raise ValueError("to exclude named locations, initial location name list can not be void and has to have the same length as one state")
1309             else:
1310                 raise ValueError(str(e))
1311     elif isinstance(__NameOfLocations, (list, numpy.ndarray, tuple)) and len(__NameOfLocations) == 0:
1312         try:
1313             __ExcludeIndexes = numpy.asarray(__ExcludeLocations, dtype=int)
1314         except ValueError as e:
1315             if "invalid literal for int() with base 10:" in str(e):
1316                 raise ValueError("to exclude named locations, initial location name list can not be void and has to have the same length as one state")
1317             else:
1318                 raise ValueError(str(e))
1319     # ----------
1320     else:
1321         try:
1322             __ExcludeIndexes = numpy.asarray(__ExcludeLocations, dtype=int)
1323         except ValueError as e:
1324             if "invalid literal for int() with base 10:" in str(e):
1325                 if len(__NameOfLocations) < 1.e6 + 1 and len(__ExcludeLocations) > 1500:
1326                     __Heuristic = True
1327                 else:
1328                     __Heuristic = False
1329                 if ForceArray or __Heuristic:
1330                     # Recherche par array permettant des noms invalides, peu efficace
1331                     __NameToIndex = dict(numpy.array((
1332                         __NameOfLocations,
1333                         range(len(__NameOfLocations))
1334                     )).T)
1335                     __ExcludeIndexes = numpy.asarray([__NameToIndex.get(k, -1) for k in __ExcludeLocations], dtype=int)
1336                     #
1337                 else:
1338                     # Recherche par liste permettant des noms invalides, très efficace
1339                     def __NameToIndex_get( cle, default = -1 ):
1340                         if cle in __NameOfLocations:
1341                             return __NameOfLocations.index(cle)
1342                         else:
1343                             return default
1344                     __ExcludeIndexes = numpy.asarray([__NameToIndex_get(k, -1) for k in __ExcludeLocations], dtype=int)
1345                     #
1346                     # Recherche par liste interdisant des noms invalides, mais encore un peu plus efficace
1347                     # __ExcludeIndexes = numpy.asarray([__NameOfLocations.index(k) for k in __ExcludeLocations], dtype=int)
1348                     #
1349                 # Ignore les noms absents
1350                 __ExcludeIndexes = numpy.compress(__ExcludeIndexes > -1, __ExcludeIndexes)
1351                 if len(__ExcludeIndexes) == 0:
1352                     __ExcludeIndexes = ()
1353             else:
1354                 raise ValueError(str(e))
1355     # ----------
1356     return __ExcludeIndexes
1357
1358 # ==============================================================================
1359 def BuildComplexSampleList(
1360         __SampleAsnUplet,
1361         __SampleAsExplicitHyperCube,
1362         __SampleAsMinMaxStepHyperCube,
1363         __SampleAsMinMaxLatinHyperCube,
1364         __SampleAsMinMaxSobolSequence,
1365         __SampleAsIndependantRandomVariables,
1366         __X0,
1367         __Seed = None ):
1368     # ---------------------------
1369     if len(__SampleAsnUplet) > 0:
1370         sampleList = __SampleAsnUplet
1371         for i, Xx in enumerate(sampleList):
1372             if numpy.ravel(Xx).size != __X0.size:
1373                 raise ValueError("The size %i of the %ith state X in the sample and %i of the checking point Xb are different, they have to be identical."%(numpy.ravel(Xx).size, i + 1, __X0.size))
1374     # ---------------------------
1375     elif len(__SampleAsExplicitHyperCube) > 0:
1376         sampleList = itertools.product(*list(__SampleAsExplicitHyperCube))
1377     # ---------------------------
1378     elif len(__SampleAsMinMaxStepHyperCube) > 0:
1379         coordinatesList = []
1380         for i, dim in enumerate(__SampleAsMinMaxStepHyperCube):
1381             if len(dim) != 3:
1382                 raise ValueError("For dimension %i, the variable definition \"%s\" is incorrect, it should be [min,max,step]."%(i, dim))
1383             else:
1384                 coordinatesList.append(numpy.linspace(dim[0], dim[1], 1 + int((float(dim[1]) - float(dim[0])) / float(dim[2]))))
1385         sampleList = itertools.product(*coordinatesList)
1386     # ---------------------------
1387     elif len(__SampleAsMinMaxLatinHyperCube) > 0:
1388         if vt(scipy.version.version) <= vt("1.7.0"):
1389             __msg = "In order to use Latin Hypercube sampling, you must at least use Scipy version 1.7.0 (and you are presently using Scipy %s). A void sample is then generated."%scipy.version.version
1390             warnings.warn(__msg, FutureWarning, stacklevel=50)
1391             sampleList = []
1392         else:
1393             __spDesc = list(__SampleAsMinMaxLatinHyperCube)
1394             __nbDime, __nbSamp  = map(int, __spDesc.pop())  # Réduction du dernier
1395             __sample = scipy.stats.qmc.LatinHypercube(
1396                 d = len(__spDesc),
1397                 seed = numpy.random.default_rng(__Seed),
1398             )
1399             __sample = __sample.random(n = __nbSamp)
1400             __bounds = numpy.array(__spDesc)[:, 0:2]
1401             __l_bounds = __bounds[:, 0]
1402             __u_bounds = __bounds[:, 1]
1403             sampleList = scipy.stats.qmc.scale(__sample, __l_bounds, __u_bounds)
1404     # ---------------------------
1405     elif len(__SampleAsMinMaxSobolSequence) > 0:
1406         if vt(scipy.version.version) <= vt("1.7.0"):
1407             __msg = "In order to use Latin Hypercube sampling, you must at least use Scipy version 1.7.0 (and you are presently using Scipy %s). A void sample is then generated."%scipy.version.version
1408             warnings.warn(__msg, FutureWarning, stacklevel=50)
1409             sampleList = []
1410         else:
1411             __spDesc = list(__SampleAsMinMaxSobolSequence)
1412             __nbDime, __nbSamp  = map(int, __spDesc.pop())  # Réduction du dernier
1413             if __nbDime != len(__spDesc):
1414                 warnings.warn("Declared space dimension (%i) is not equal to number of bounds (%i), the last one will be used."%(__nbDime, len(__spDesc)), FutureWarning, stacklevel=50)
1415             __sample = scipy.stats.qmc.Sobol(
1416                 d = len(__spDesc),
1417                 seed = numpy.random.default_rng(__Seed),
1418             )
1419             __sample = __sample.random_base2(m = int(math.log2(__nbSamp)) + 1)
1420             __bounds = numpy.array(__spDesc)[:, 0:2]
1421             __l_bounds = __bounds[:, 0]
1422             __u_bounds = __bounds[:, 1]
1423             sampleList = scipy.stats.qmc.scale(__sample, __l_bounds, __u_bounds)
1424     # ---------------------------
1425     elif len(__SampleAsIndependantRandomVariables) > 0:
1426         coordinatesList = []
1427         for i, dim in enumerate(__SampleAsIndependantRandomVariables):
1428             if len(dim) != 3:
1429                 raise ValueError("For dimension %i, the variable definition \"%s\" is incorrect, it should be ('distribution',(parameters),length) with distribution in ['normal'(mean,std),'lognormal'(mean,sigma),'uniform'(low,high),'weibull'(shape)]."%(i, dim))
1430             elif not ( str(dim[0]) in ['normal', 'lognormal', 'uniform', 'weibull'] \
1431                        and hasattr(numpy.random, str(dim[0])) ):
1432                 raise ValueError("For dimension %i, the distribution name \"%s\" is not allowed, please choose in ['normal'(mean,std),'lognormal'(mean,sigma),'uniform'(low,high),'weibull'(shape)]"%(i, str(dim[0])))
1433             else:
1434                 distribution = getattr(numpy.random, str(dim[0]), 'normal')
1435                 coordinatesList.append(distribution(*dim[1], size=max(1, int(dim[2]))))
1436         sampleList = itertools.product(*coordinatesList)
1437     else:
1438         sampleList = iter([__X0,])
1439     # ----------
1440     return sampleList
1441
1442 # ==============================================================================
1443 def multiXOsteps(
1444         selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle,
1445         __CovForecast = False ):
1446     """
1447     Prévision multi-pas avec une correction par pas (multi-méthodes)
1448     """
1449     #
1450     # Initialisation
1451     # --------------
1452     if selfA._parameters["EstimationOf"] == "State":
1453         if len(selfA.StoredVariables["Analysis"]) == 0 or not selfA._parameters["nextStep"]:
1454             Xn = numpy.asarray(Xb)
1455             if __CovForecast:
1456                 Pn = B
1457             selfA.StoredVariables["Analysis"].store( Xn )
1458             if selfA._toStore("APosterioriCovariance"):
1459                 if hasattr(B, "asfullmatrix"):
1460                     selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
1461                 else:
1462                     selfA.StoredVariables["APosterioriCovariance"].store( B )
1463             selfA._setInternalState("seed", numpy.random.get_state())
1464         elif selfA._parameters["nextStep"]:
1465             Xn = selfA._getInternalState("Xn")
1466             if __CovForecast:
1467                 Pn = selfA._getInternalState("Pn")
1468     else:
1469         Xn = numpy.asarray(Xb)
1470         if __CovForecast:
1471             Pn = B
1472     #
1473     if hasattr(Y, "stepnumber"):
1474         duration = Y.stepnumber()
1475     else:
1476         duration = 2
1477     #
1478     # Multi-steps
1479     # -----------
1480     for step in range(duration - 1):
1481         selfA.StoredVariables["CurrentStepNumber"].store( len(selfA.StoredVariables["Analysis"]) )
1482         #
1483         if hasattr(Y, "store"):
1484             Ynpu = numpy.asarray( Y[step + 1] ).reshape((-1, 1))
1485         else:
1486             Ynpu = numpy.asarray( Y ).reshape((-1, 1))
1487         #
1488         if U is not None:
1489             if hasattr(U, "store") and len(U) > 1:
1490                 Un = numpy.asarray( U[step] ).reshape((-1, 1))
1491             elif hasattr(U, "store") and len(U) == 1:
1492                 Un = numpy.asarray( U[0] ).reshape((-1, 1))
1493             else:
1494                 Un = numpy.asarray( U ).reshape((-1, 1))
1495         else:
1496             Un = None
1497         #
1498         # Predict (Time Update)
1499         # ---------------------
1500         if selfA._parameters["EstimationOf"] == "State":
1501             if __CovForecast:
1502                 Mt = EM["Tangent"].asMatrix(Xn)
1503                 Mt = Mt.reshape(Xn.size, Xn.size)  # ADAO & check shape
1504                 Ma = EM["Adjoint"].asMatrix(Xn)
1505                 Ma = Ma.reshape(Xn.size, Xn.size)  # ADAO & check shape
1506                 Pn_predicted = Q + Mt @ (Pn @ Ma)
1507             Mm = EM["Direct"].appliedControledFormTo
1508             Xn_predicted = Mm( (Xn, Un) ).reshape((-1, 1))
1509             if CM is not None and "Tangent" in CM and Un is not None:  # Attention : si Cm est aussi dans M, doublon !
1510                 Cm = CM["Tangent"].asMatrix(Xn_predicted)
1511                 Cm = Cm.reshape(Xn.size, Un.size)  # ADAO & check shape
1512                 Xn_predicted = Xn_predicted + (Cm @ Un).reshape((-1, 1))
1513         elif selfA._parameters["EstimationOf"] == "Parameters":  # No forecast
1514             # --- > Par principe, M = Id, Q = 0
1515             Xn_predicted = Xn
1516             if __CovForecast:
1517                 Pn_predicted = Pn
1518         Xn_predicted = numpy.asarray(Xn_predicted).reshape((-1, 1))
1519         if selfA._toStore("ForecastState"):
1520             selfA.StoredVariables["ForecastState"].store( Xn_predicted )
1521         if __CovForecast:
1522             if hasattr(Pn_predicted, "asfullmatrix"):
1523                 Pn_predicted = Pn_predicted.asfullmatrix(Xn.size)
1524             else:
1525                 Pn_predicted = numpy.asarray(Pn_predicted).reshape((Xn.size, Xn.size))
1526             if selfA._toStore("ForecastCovariance"):
1527                 selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
1528         #
1529         # Correct (Measurement Update)
1530         # ----------------------------
1531         if __CovForecast:
1532             oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, Pn_predicted, True)
1533         else:
1534             oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, B, True)
1535         #
1536         # --------------------------
1537         Xn = selfA._getInternalState("Xn")
1538         if __CovForecast:
1539             Pn = selfA._getInternalState("Pn")
1540     #
1541     return 0
1542
1543 # ==============================================================================
1544 if __name__ == "__main__":
1545     print("\n AUTODIAGNOSTIC\n")