src/daComposant/daCore/NumericObjects.py

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