src/daComposant/daCore/NumericObjects.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2020 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Définit les objets numériques génériques.
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import os, time, copy, types, sys, logging
  29 import math, numpy, scipy
  30 from daCore.BasicObjects import Operator
  31 from daCore.PlatformInfo import PlatformInfo
  32 mpr = PlatformInfo().MachinePrecision()
  33 mfp = PlatformInfo().MaximumPrecision()
  34 # logging.getLogger().setLevel(logging.DEBUG)
  35
  36 # ==============================================================================
  37 def ExecuteFunction( paire ):
  38     assert len(paire) == 2, "Incorrect number of arguments"
  39     X, funcrepr = paire
  40     __X = numpy.asmatrix(numpy.ravel( X )).T
  41     __sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
  42     __module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
  43     __fonction = getattr(__module,funcrepr["__userFunction__name"])
  44     sys.path = __sys_path_tmp ; del __sys_path_tmp
  45     __HX  = __fonction( __X )
  46     return numpy.ravel( __HX )
  47
  48 # ==============================================================================
  49 class FDApproximation(object):
  50     """
  51     Cette classe sert d'interface pour définir les opérateurs approximés. A la
  52     création d'un objet, en fournissant une fonction "Function", on obtient un
  53     objet qui dispose de 3 méthodes "DirectOperator", "TangentOperator" et
  54     "AdjointOperator". On contrôle l'approximation DF avec l'incrément
  55     multiplicatif "increment" valant par défaut 1%, ou avec l'incrément fixe
  56     "dX" qui sera multiplié par "increment" (donc en %), et on effectue de DF
  57     centrées si le booléen "centeredDF" est vrai.
  58     """
  59     def __init__(self,
  60             name                  = "FDApproximation",
  61             Function              = None,
  62             centeredDF            = False,
  63             increment             = 0.01,
  64             dX                    = None,
  65             avoidingRedundancy    = True,
  66             toleranceInRedundancy = 1.e-18,
  67             lenghtOfRedundancy    = -1,
  68             mpEnabled             = False,
  69             mpWorkers             = None,
  70             mfEnabled             = False,
  71             ):
  72         self.__name = str(name)
  73         if mpEnabled:
  74             try:
  75                 import multiprocessing
  76                 self.__mpEnabled = True
  77             except ImportError:
  78                 self.__mpEnabled = False
  79         else:
  80             self.__mpEnabled = False
  81         self.__mpWorkers = mpWorkers
  82         if self.__mpWorkers is not None and self.__mpWorkers < 1:
  83             self.__mpWorkers = None
  84         logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
  85         #
  86         if mfEnabled:
  87             self.__mfEnabled = True
  88         else:
  89             self.__mfEnabled = False
  90         logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
  91         #
  92         if avoidingRedundancy:
  93             self.__avoidRC = True
  94             self.__tolerBP = float(toleranceInRedundancy)
  95             self.__lenghtRJ = int(lenghtOfRedundancy)
  96             self.__listJPCP = [] # Jacobian Previous Calculated Points
  97             self.__listJPCI = [] # Jacobian Previous Calculated Increment
  98             self.__listJPCR = [] # Jacobian Previous Calculated Results
  99             self.__listJPPN = [] # Jacobian Previous Calculated Point Norms
 100             self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
 101         else:
 102             self.__avoidRC = False
 103         #
 104         if self.__mpEnabled:
 105             if isinstance(Function,types.FunctionType):
 106                 logging.debug("FDA Calculs en multiprocessing : FunctionType")
 107                 self.__userFunction__name = Function.__name__
 108                 try:
 109                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 110                 except:
 111                     mod = os.path.abspath(Function.__globals__['__file__'])
 112                 if not os.path.isfile(mod):
 113                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 114                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 115                 self.__userFunction__path = os.path.dirname(mod)
 116                 del mod
 117                 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 118                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 119             elif isinstance(Function,types.MethodType):
 120                 logging.debug("FDA Calculs en multiprocessing : MethodType")
 121                 self.__userFunction__name = Function.__name__
 122                 try:
 123                     mod = os.path.join(Function.__globals__['filepath'],Function.__globals__['filename'])
 124                 except:
 125                     mod = os.path.abspath(Function.__func__.__globals__['__file__'])
 126                 if not os.path.isfile(mod):
 127                     raise ImportError("No user defined function or method found with the name %s"%(mod,))
 128                 self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
 129                 self.__userFunction__path = os.path.dirname(mod)
 130                 del mod
 131                 self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 132                 self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
 133             else:
 134                 raise TypeError("User defined function or method has to be provided for finite differences approximation.")
 135         else:
 136             self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
 137             self.__userFunction = self.__userOperator.appliedTo
 138         #
 139         self.__centeredDF = bool(centeredDF)
 140         if abs(float(increment)) > 1.e-15:
 141             self.__increment  = float(increment)
 142         else:
 143             self.__increment  = 0.01
 144         if dX is None:
 145             self.__dX     = None
 146         else:
 147             self.__dX     = numpy.asmatrix(numpy.ravel( dX )).T
 148         logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
 149         if self.__avoidRC:
 150             logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
 151
 152     # ---------------------------------------------------------
 153     def __doublon__(self, e, l, n, v=None):
 154         __ac, __iac = False, -1
 155         for i in range(len(l)-1,-1,-1):
 156             if numpy.linalg.norm(e - l[i]) < self.__tolerBP * n[i]:
 157                 __ac, __iac = True, i
 158                 if v is not None: logging.debug("FDA Cas%s déja calculé, récupération du doublon %i"%(v,__iac))
 159                 break
 160         return __ac, __iac
 161
 162     # ---------------------------------------------------------
 163     def DirectOperator(self, X ):
 164         """
 165         Calcul du direct à l'aide de la fonction fournie.
 166         """
 167         logging.debug("FDA Calcul DirectOperator (explicite)")
 168         if self.__mfEnabled:
 169             _HX = self.__userFunction( X, argsAsSerie = True )
 170         else:
 171             _X = numpy.asmatrix(numpy.ravel( X )).T
 172             _HX = numpy.ravel(self.__userFunction( _X ))
 173         #
 174         return _HX
 175
 176     # ---------------------------------------------------------
 177     def TangentMatrix(self, X ):
 178         """
 179         Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
 180         c'est-à-dire le gradient de H en X. On utilise des différences finies
 181         directionnelles autour du point X. X est un numpy.matrix.
 182
 183         Différences finies centrées (approximation d'ordre 2):
 184         1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
 185            dX[i] à la  composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
 186            on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
 187            H( X_moins_dXi )
 188         2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
 189            le pas 2*dXi
 190         3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 191
 192         Différences finies non centrées (approximation d'ordre 1):
 193         1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
 194            composante X[i] pour composer X_plus_dXi, et on calcule la réponse
 195            HX_plus_dXi = H( X_plus_dXi )
 196         2/ On calcule la valeur centrale HX = H(X)
 197         3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
 198            le pas dXi
 199         4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
 200
 201         """
 202         logging.debug("FDA Début du calcul de la Jacobienne")
 203         logging.debug("FDA   Incrément de............: %s*X"%float(self.__increment))
 204         logging.debug("FDA   Approximation centrée...: %s"%(self.__centeredDF))
 205         #
 206         if X is None or len(X)==0:
 207             raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
 208         #
 209         _X = numpy.asmatrix(numpy.ravel( X )).T
 210         #
 211         if self.__dX is None:
 212             _dX  = self.__increment * _X
 213         else:
 214             _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
 215         #
 216         if (_dX == 0.).any():
 217             moyenne = _dX.mean()
 218             if moyenne == 0.:
 219                 _dX = numpy.where( _dX == 0., float(self.__increment), _dX )
 220             else:
 221                 _dX = numpy.where( _dX == 0., moyenne, _dX )
 222         #
 223         __alreadyCalculated  = False
 224         if self.__avoidRC:
 225             __bidon, __alreadyCalculatedP = self.__doublon__(_X,  self.__listJPCP, self.__listJPPN, None)
 226             __bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
 227             if __alreadyCalculatedP == __alreadyCalculatedI > -1:
 228                 __alreadyCalculated, __i = True, __alreadyCalculatedP
 229                 logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
 230         #
 231         if __alreadyCalculated:
 232             logging.debug("FDA   Calcul Jacobienne (par récupération du doublon %i)"%__i)
 233             _Jacobienne = self.__listJPCR[__i]
 234         else:
 235             logging.debug("FDA   Calcul Jacobienne (explicite)")
 236             if self.__centeredDF:
 237                 #
 238                 if self.__mpEnabled and not self.__mfEnabled:
 239                     funcrepr = {
 240                         "__userFunction__path" : self.__userFunction__path,
 241                         "__userFunction__modl" : self.__userFunction__modl,
 242                         "__userFunction__name" : self.__userFunction__name,
 243                     }
 244                     _jobs = []
 245                     for i in range( len(_dX) ):
 246                         _dXi            = _dX[i]
 247                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 248                         _X_plus_dXi[i]  = _X[i] + _dXi
 249                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 250                         _X_moins_dXi[i] = _X[i] - _dXi
 251                         #
 252                         _jobs.append( (_X_plus_dXi,  funcrepr) )
 253                         _jobs.append( (_X_moins_dXi, funcrepr) )
 254                     #
 255                     import multiprocessing
 256                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 257                     _HX_plusmoins_dX = self.__pool.map( ExecuteFunction, _jobs )
 258                     self.__pool.close()
 259                     self.__pool.join()
 260                     #
 261                     _Jacobienne  = []
 262                     for i in range( len(_dX) ):
 263                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 264                     #
 265                 elif self.__mfEnabled:
 266                     _xserie = []
 267                     for i in range( len(_dX) ):
 268                         _dXi            = _dX[i]
 269                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 270                         _X_plus_dXi[i]  = _X[i] + _dXi
 271                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 272                         _X_moins_dXi[i] = _X[i] - _dXi
 273                         #
 274                         _xserie.append( _X_plus_dXi )
 275                         _xserie.append( _X_moins_dXi )
 276                     #
 277                     _HX_plusmoins_dX = self.DirectOperator( _xserie )
 278                      #
 279                     _Jacobienne  = []
 280                     for i in range( len(_dX) ):
 281                         _Jacobienne.append( numpy.ravel( _HX_plusmoins_dX[2*i] - _HX_plusmoins_dX[2*i+1] ) / (2.*_dX[i]) )
 282                     #
 283                 else:
 284                     _Jacobienne  = []
 285                     for i in range( _dX.size ):
 286                         _dXi            = _dX[i]
 287                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 288                         _X_plus_dXi[i]  = _X[i] + _dXi
 289                         _X_moins_dXi    = numpy.array( _X.A1, dtype=float )
 290                         _X_moins_dXi[i] = _X[i] - _dXi
 291                         #
 292                         _HX_plus_dXi    = self.DirectOperator( _X_plus_dXi )
 293                         _HX_moins_dXi   = self.DirectOperator( _X_moins_dXi )
 294                         #
 295                         _Jacobienne.append( numpy.ravel( _HX_plus_dXi - _HX_moins_dXi ) / (2.*_dXi) )
 296                 #
 297             else:
 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                     _jobs.append( (_X.A1, funcrepr) )
 307                     for i in range( len(_dX) ):
 308                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 309                         _X_plus_dXi[i] = _X[i] + _dX[i]
 310                         #
 311                         _jobs.append( (_X_plus_dXi, funcrepr) )
 312                     #
 313                     import multiprocessing
 314                     self.__pool = multiprocessing.Pool(self.__mpWorkers)
 315                     _HX_plus_dX = self.__pool.map( ExecuteFunction, _jobs )
 316                     self.__pool.close()
 317                     self.__pool.join()
 318                     #
 319                     _HX = _HX_plus_dX.pop(0)
 320                     #
 321                     _Jacobienne = []
 322                     for i in range( len(_dX) ):
 323                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 324                     #
 325                 elif self.__mfEnabled:
 326                     _xserie = []
 327                     _xserie.append( _X.A1 )
 328                     for i in range( len(_dX) ):
 329                         _X_plus_dXi    = numpy.array( _X.A1, dtype=float )
 330                         _X_plus_dXi[i] = _X[i] + _dX[i]
 331                         #
 332                         _xserie.append( _X_plus_dXi )
 333                     #
 334                     _HX_plus_dX = self.DirectOperator( _xserie )
 335                     #
 336                     _HX = _HX_plus_dX.pop(0)
 337                     #
 338                     _Jacobienne = []
 339                     for i in range( len(_dX) ):
 340                         _Jacobienne.append( numpy.ravel(( _HX_plus_dX[i] - _HX ) / _dX[i]) )
 341                    #
 342                 else:
 343                     _Jacobienne  = []
 344                     _HX = self.DirectOperator( _X )
 345                     for i in range( _dX.size ):
 346                         _dXi            = _dX[i]
 347                         _X_plus_dXi     = numpy.array( _X.A1, dtype=float )
 348                         _X_plus_dXi[i]  = _X[i] + _dXi
 349                         #
 350                         _HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
 351                         #
 352                         _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
 353                 #
 354             #
 355             _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
 356             if self.__avoidRC:
 357                 if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
 358                 while len(self.__listJPCP) > self.__lenghtRJ:
 359                     self.__listJPCP.pop(0)
 360                     self.__listJPCI.pop(0)
 361                     self.__listJPCR.pop(0)
 362                     self.__listJPPN.pop(0)
 363                     self.__listJPIN.pop(0)
 364                 self.__listJPCP.append( copy.copy(_X) )
 365                 self.__listJPCI.append( copy.copy(_dX) )
 366                 self.__listJPCR.append( copy.copy(_Jacobienne) )
 367                 self.__listJPPN.append( numpy.linalg.norm(_X) )
 368                 self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
 369         #
 370         logging.debug("FDA Fin du calcul de la Jacobienne")
 371         #
 372         return _Jacobienne
 373
 374     # ---------------------------------------------------------
 375     def TangentOperator(self, paire ):
 376         """
 377         Calcul du tangent à l'aide de la Jacobienne.
 378         """
 379         if self.__mfEnabled:
 380             assert len(paire) == 1, "Incorrect lenght of arguments"
 381             _paire = paire[0]
 382             assert len(_paire) == 2, "Incorrect number of arguments"
 383         else:
 384             assert len(paire) == 2, "Incorrect number of arguments"
 385             _paire = paire
 386         X, dX = _paire
 387         _Jacobienne = self.TangentMatrix( X )
 388         if dX is None or len(dX) == 0:
 389             #
 390             # Calcul de la forme matricielle si le second argument est None
 391             # -------------------------------------------------------------
 392             if self.__mfEnabled: return [_Jacobienne,]
 393             else:                return _Jacobienne
 394         else:
 395             #
 396             # Calcul de la valeur linéarisée de H en X appliqué à dX
 397             # ------------------------------------------------------
 398             _dX = numpy.asmatrix(numpy.ravel( dX )).T
 399             _HtX = numpy.dot(_Jacobienne, _dX)
 400             if self.__mfEnabled: return [_HtX.A1,]
 401             else:                return _HtX.A1
 402
 403     # ---------------------------------------------------------
 404     def AdjointOperator(self, paire ):
 405         """
 406         Calcul de l'adjoint à l'aide de la Jacobienne.
 407         """
 408         if self.__mfEnabled:
 409             assert len(paire) == 1, "Incorrect lenght of arguments"
 410             _paire = paire[0]
 411             assert len(_paire) == 2, "Incorrect number of arguments"
 412         else:
 413             assert len(paire) == 2, "Incorrect number of arguments"
 414             _paire = paire
 415         X, Y = _paire
 416         _JacobienneT = self.TangentMatrix( X ).T
 417         if Y is None or len(Y) == 0:
 418             #
 419             # Calcul de la forme matricielle si le second argument est None
 420             # -------------------------------------------------------------
 421             if self.__mfEnabled: return [_JacobienneT,]
 422             else:                return _JacobienneT
 423         else:
 424             #
 425             # Calcul de la valeur de l'adjoint en X appliqué à Y
 426             # --------------------------------------------------
 427             _Y = numpy.asmatrix(numpy.ravel( Y )).T
 428             _HaY = numpy.dot(_JacobienneT, _Y)
 429             if self.__mfEnabled: return [_HaY.A1,]
 430             else:                return _HaY.A1
 431
 432 # ==============================================================================
 433 def mmqr(
 434         func     = None,
 435         x0       = None,
 436         fprime   = None,
 437         bounds   = None,
 438         quantile = 0.5,
 439         maxfun   = 15000,
 440         toler    = 1.e-06,
 441         y        = None,
 442         ):
 443     """
 444     Implémentation informatique de l'algorithme MMQR, basée sur la publication :
 445     David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
 446     Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
 447     """
 448     #
 449     # Recuperation des donnees et informations initiales
 450     # --------------------------------------------------
 451     variables = numpy.ravel( x0 )
 452     mesures   = numpy.ravel( y )
 453     increment = sys.float_info[0]
 454     p         = variables.size
 455     n         = mesures.size
 456     quantile  = float(quantile)
 457     #
 458     # Calcul des parametres du MM
 459     # ---------------------------
 460     tn      = float(toler) / n
 461     e0      = -tn / math.log(tn)
 462     epsilon = (e0-tn)/(1+math.log(e0))
 463     #
 464     # Calculs d'initialisation
 465     # ------------------------
 466     residus  = mesures - numpy.ravel( func( variables ) )
 467     poids    = 1./(epsilon+numpy.abs(residus))
 468     veps     = 1. - 2. * quantile - residus * poids
 469     lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
 470     iteration = 0
 471     #
 472     # Recherche iterative
 473     # -------------------
 474     while (increment > toler) and (iteration < maxfun) :
 475         iteration += 1
 476         #
 477         Derivees  = numpy.array(fprime(variables))
 478         Derivees  = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
 479         DeriveesT = Derivees.transpose()
 480         M         =   numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
 481         SM        =   numpy.transpose(numpy.dot( DeriveesT , veps ))
 482         step      = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
 483         #
 484         variables = variables + step
 485         if bounds is not None:
 486             while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
 487                 step      = step/2.
 488                 variables = variables - step
 489         residus   = mesures - numpy.ravel( func(variables) )
 490         surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
 491         #
 492         while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
 493             step      = step/2.
 494             variables = variables - step
 495             residus   = mesures - numpy.ravel( func(variables) )
 496             surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
 497         #
 498         increment     = lastsurrogate-surrogate
 499         poids         = 1./(epsilon+numpy.abs(residus))
 500         veps          = 1. - 2. * quantile - residus * poids
 501         lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
 502     #
 503     # Mesure d'écart
 504     # --------------
 505     Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
 506     #
 507     return variables, Ecart, [n,p,iteration,increment,0]
 508
 509 # ==============================================================================
 510
 511 def _BackgroundEnsembleGeneration( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True):
 512     "Génération d'un ensemble d'ébauche de taille _nbmembers-1"
 513     # ~ numpy.random.seed(1234567)
 514     if _nbmembers < 1:
 515         raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),))
 516     if _withSVD:
 517         U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False)
 518         _nbctl = len(_bgcenter)
 519         if _nbmembers > _nbctl:
 520             _Z = numpy.concatenate((numpy.dot(
 521                 numpy.diag(numpy.sqrt(s[:_nbctl])), V[:_nbctl]),
 522                 numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1-_nbctl)), axis = 0)
 523         else:
 524             _Z = numpy.dot(numpy.diag(numpy.sqrt(s[:_nbmembers-1])), V[:_nbmembers-1])
 525         _Zca = _CenteredAnomalies(_Z, _nbmembers)
 526         BackgroundEnsemble = (_bgcenter + _Zca.T).T
 527     else:
 528         if max(abs(_bgcovariance.flatten())) > 0:
 529             _nbctl = len(_bgcenter)
 530             _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),_bgcovariance,_nbmembers-1)
 531             _Zca = _CenteredAnomalies(_Z, _nbmembers)
 532             BackgroundEnsemble = (_bgcenter + _Zca.T).T
 533         else:
 534             BackgroundEnsemble = numpy.tile([_bgcenter],(_nbmembers,1)).T
 535     return BackgroundEnsemble
 536
 537 def _CenteredAnomalies(Zr, N):
 538     """
 539     Génère une matrice d'anomalies centrées selon les notes manuscrites de MB
 540     et conforme au code de PS avec eps = -1
 541     """
 542     eps = -1
 543     Q = numpy.eye(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
 544     Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
 545     R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
 546     Q = numpy.dot(Q,R)
 547     Zr = numpy.dot(Q,Zr)
 548     return Zr.T
 549
 550 def _IEnKF_cycle_Lag_1_SDA_GN(
 551         E0         = None,
 552         yObs       = None,
 553         RIdemi     = None,
 554         Mnnpu      = None,
 555         Hn         = None,
 556         variant    = "IEnKF", # IEnKF or IEKF
 557         iMaximum   = 15000,
 558         sTolerance = mfp,
 559         jTolerance = mfp,
 560         epsilonE   = 1e-5,
 561         nbPS       = 0,  # nbPreviousSteps
 562         ):
 563     # 201206
 564     if logging.getLogger().level < logging.WARNING:
 565         assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
 566         assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
 567         assert variant in ("IEnKF", "IEKF"), "Variant has to be IEnKF or IEKF"
 568     #
 569     nbCtl, nbMbr = E0.shape
 570     nbObs = yObs.size
 571     #
 572     if logging.getLogger().level < logging.WARNING:
 573         assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
 574     #
 575     yo  = yObs.reshape((nbObs,1))
 576     IN  = numpy.identity(nbMbr)
 577     if variant == "IEnKF":
 578         T    = numpy.identity(nbMbr)
 579         Tinv = numpy.identity(nbMbr)
 580     x00 = numpy.mean(E0, axis = 1)
 581     Ah0 = E0 - x00
 582     Ap0 = numpy.linalg.pinv( Ah0.T.dot(Ah0) )
 583     if logging.getLogger().level < logging.WARNING:
 584         assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
 585         assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
 586         assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
 587     #
 588     def _convergence_condition(j, dx, JCurr, JPrev):
 589         if j > iMaximum:
 590             logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
 591             return True
 592         #---------
 593         if j == 1:
 594             _deltaOnJ = 1.
 595         else:
 596             _deltaOnJ = abs(JCurr - JPrev) / JPrev
 597         if _deltaOnJ <= jTolerance:
 598             logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
 599             return True
 600         #---------
 601         _deltaOnX = numpy.linalg.norm(dx)
 602         if _deltaOnX <= sTolerance:
 603             logging.debug("Convergence on norm of state correction, that is below the threshold of %.1e."%sTolerance)
 604             return True # En correction de l'état
 605         #---------
 606         return False
 607     #
 608     St = dict([(k,[]) for k in [
 609         "CurrentState", "CurrentEnsemble",
 610         "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
 611         ]])
 612     #
 613     j, convergence, JPrev = 1, False, numpy.nan
 614     x1 = x00
 615     while not convergence:
 616         logging.debug("Internal IEnKS step number %i"%j)
 617         St["CurrentState"].append( x1.squeeze() )
 618         if variant == "IEnKF": # Transform
 619             E1 = x1 + Ah0.dot(T)
 620         else: # IEKF
 621             E1 = x1 + epsilonE * Ah0
 622         St["CurrentEnsemble"].append( E1 )
 623         E2  = numpy.array([Mnnpu(_x) for _x in E1.T]).reshape((nbCtl, nbMbr)) # Evolution 1->2
 624         HEL = numpy.array([Hn(_x) for _x in E2.T]).T     # Observation à 2
 625         yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
 626         HA2 = HEL - yLm
 627         if variant == "IEnKF":
 628             HA2 = HA2.dot(Tinv)
 629         else:
 630             HA2 = HA2 / epsilonE
 631         RIdemidy = RIdemi.dot(yo - yLm)
 632         xs = RIdemidy / math.sqrt(nbMbr-1)
 633         ES = RIdemi.dot(HA2) / math.sqrt(nbMbr-1)
 634         G  = numpy.linalg.inv(IN + ES.T.dot(ES))
 635         xb = G.dot(ES.T.dot(xs))
 636         dx = Ah0.dot(xb) + Ah0.dot(G.dot(Ap0.dot(Ah0.T.dot(x00 - x1))))
 637         #
 638         Jb = float(dx.T.dot(dx))
 639         Jo = float(RIdemidy.T.dot(RIdemidy))
 640         J  = Jo + Jb
 641         logging.debug("Values for cost functions are: J = %.5e  Jo = %.5e  Jb = %.5e"%(J,Jo,Jb))
 642         St["CostFunctionJb"].append( Jb )
 643         St["CostFunctionJo"].append( Jo )
 644         St["CostFunctionJ"].append( J )
 645         #
 646         x1 = x1 + dx
 647         j = j + 1
 648         convergence = _convergence_condition(j, dx, J, JPrev)
 649         JPrev = J
 650         #
 651         if variant == "IEnKF":
 652             T    = numpy.real_if_close(scipy.linalg.sqrtm(G))
 653             Tinv = numpy.linalg.inv(T)
 654     #
 655     # Stocke le dernier pas
 656     x2 = numpy.mean( E2, axis = 1)
 657     if variant == "IEKF":
 658         A2 = E2 - x2
 659         A2 = A2.dot(numpy.linalg.cholesky(G)) / epsilonE
 660         E2 = x2 + A2
 661     St["CurrentState"].append( x2.squeeze() )
 662     St["CurrentEnsemble"].append( E2 )
 663     #
 664     IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
 665     xa = St["CurrentState"][IndexMin]
 666     Ea = St["CurrentEnsemble"][IndexMin]
 667     #
 668     return (xa, Ea, St)
 669
 670 def ienkf(
 671         xb         = None,          # Background (None si E0)
 672         E0         = None,          # Background ensemble (None si xb)
 673         yObs       = None,          # Observation (série)
 674         B          = None,          # B
 675         RIdemi     = None,          # R^(-1/2)
 676         Mnnpu      = None,          # Evolution operator
 677         Hn         = None,          # Observation operator
 678         variant    = "IEnKF",       # IEnKF or IEKF
 679         nMembers   = 5,             # Number of members
 680         sMaximum   = 0,             # Number of spinup steps
 681         cMaximum   = 15000,         # Number of steps or cycles
 682         iMaximum   = 15000,         # Number of iterations per cycle
 683         sTolerance = mfp,           # State correction tolerance
 684         jTolerance = mfp,           # Cost decrement tolerance
 685         epsilon    = 1e-5,
 686         inflation  = 1.,
 687         nbPS       = 0,             # Number of previous steps
 688         setSeed    = None,
 689         ):
 690     #
 691     # Initial
 692     if setSeed is not None: numpy.random.seed(setSeed)
 693     if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
 694     #
 695     # Spinup
 696     # ------
 697     #
 698     # Cycles
 699     # ------
 700     xa, Ea, Sa = [xb,], [E0,], [{}]
 701     for step in range(cMaximum):
 702         if hasattr(yObs,"store"):         Ynpu = numpy.ravel( yObs[step+1] )
 703         elif type(yObs) in [list, tuple]: Ynpu = numpy.ravel( yObs[step+1] )
 704         else:                             Ynpu = numpy.ravel( yObs )
 705         #
 706         (xa_c, Ea_c, Sa_c) = _IEnKF_cycle_Lag_1_SDA_GN(
 707             E0,
 708             Ynpu,
 709             RIdemi,
 710             Mnnpu,
 711             Hn,
 712             variant,
 713             iMaximum,
 714             sTolerance,
 715             jTolerance,
 716             epsilon,
 717             nbPS,
 718             )
 719         xa.append( xa_c )
 720         Ea.append( Ea_c )
 721         Sa.append( Sa_c )
 722         #
 723         # Inflation for next cycle
 724         E0 = xa_c + inflation * (Ea_c - xa_c)
 725     #
 726     return (xa, Ea, Sa)
 727
 728 def _IEnKS_cycle_Lag_L_SDA_GN(
 729         E0         = None,
 730         yObs       = None,
 731         RIdemi     = None,
 732         Mnnpu      = None,
 733         Hn         = None,
 734         method     = "Transform",
 735         iMaximum   = 15000,
 736         sTolerance = mfp,
 737         jTolerance = mfp,
 738         Lag        = 1,
 739         epsilon    = -1.,
 740         nbPS       = 0,
 741         ):
 742     # 201407 & 201905
 743     if logging.getLogger().level < logging.WARNING:
 744         assert len(E0.shape) == 2, "Ensemble E0 is not well formed: not of shape 2!"
 745         assert len(RIdemi.shape) == 2, "R^{-1/2} is not well formed: not of shape 2!"
 746         assert method in ("Transform", "Bundle"), "Method has to be Transform or Bundle"
 747     #
 748     nbCtl, nbMbr = E0.shape
 749     nbObs = yObs.size
 750     #
 751     if logging.getLogger().level < logging.WARNING:
 752         assert RIdemi.shape[0] == RIdemi.shape[1] == nbObs, "R^{-1} not of good size: not of size nbObs!"
 753     #
 754     yo  = yObs.reshape((nbObs,1))
 755     IN  = numpy.identity(nbMbr)
 756     if method == "Transform":
 757         T    = numpy.identity(nbMbr)
 758         Tinv = numpy.identity(nbMbr)
 759     x00 = numpy.mean(E0, axis = 1)
 760     Ah0 = E0 - x00
 761     Am0  = (1/math.sqrt(nbMbr - 1)) * Ah0
 762     w   = numpy.zeros((nbMbr,1))
 763     if logging.getLogger().level < logging.WARNING:
 764         assert len(Ah0.shape) == 2, "Ensemble A0 is not well formed, of shape 2!"
 765         assert Ah0.shape[0] == nbCtl and Ah0.shape[1] == nbMbr, "Ensemble A0 is not well shaped!"
 766         assert abs(max(numpy.mean(Ah0, axis = 1))) < nbMbr*mpr, "Ensemble A0 seems not to be centered!"
 767     #
 768     def _convergence_condition(j, dw, JCurr, JPrev):
 769         if j > iMaximum:
 770             logging.debug("Convergence on maximum number of iterations per cycle, that reach the limit of %i."%iMaximum)
 771             return True
 772         #---------
 773         if j == 1:
 774             _deltaOnJ = 1.
 775         else:
 776             _deltaOnJ = abs(JCurr - JPrev) / JPrev
 777         if _deltaOnJ <= jTolerance:
 778             logging.debug("Convergence on cost decrement tolerance, that is below the threshold of %.1e."%jTolerance)
 779             return True
 780         #---------
 781         _deltaOnW = numpy.sqrt(numpy.mean(dw.squeeze()**2))
 782         if _deltaOnW <= sTolerance:
 783             logging.debug("Convergence on norm of weights correction, that is below the threshold of %.1e."%sTolerance)
 784             return True # En correction des poids
 785         #---------
 786         return False
 787     #
 788     St = dict([(k,[]) for k in [
 789         "CurrentState", "CurrentEnsemble", "CurrentWeights",
 790         "CostFunctionJb", "CostFunctionJo", "CostFunctionJ",
 791         ]])
 792     #
 793     j, convergence, JPrev = 1, False, numpy.nan
 794     while not convergence:
 795         logging.debug("Internal IEnKS step number %i"%j)
 796         x0 = x00 + Am0.dot( w )
 797         St["CurrentState"].append( x0.squeeze() )
 798         if method == "Transform":
 799             E0 = x0 + Ah0.dot(T)
 800         else:
 801             E0 = x0 + epsilon * Am0
 802         St["CurrentEnsemble"].append( E0 )
 803         Ek = E0
 804         yHmean = numpy.mean(E0, axis = 1)
 805         for k in range(1, Lag+1):
 806             Ek  = numpy.array([Mnnpu(_x) for _x in Ek.T]).reshape((nbCtl, nbMbr)) # Evolution 0->L
 807             if method == "Transform":
 808                 yHmean = Mnnpu(yHmean)
 809         HEL = numpy.array([Hn(_x) for _x in Ek.T]).T     # Observation à L
 810         #
 811         if method == "Transform":
 812             yLm = Hn( yHmean ).reshape((nbObs,1))
 813             YL = RIdemi.dot( (HEL - numpy.mean( HEL, axis = 1).reshape((nbObs,1))).dot(Tinv) ) / math.sqrt(nbMbr-1)
 814         else:
 815             yLm = numpy.mean( HEL, axis = 1).reshape((nbObs,1))
 816             YL = RIdemi.dot(HEL - yLm) / epsilon
 817         dy = RIdemi.dot(yo - yLm)
 818         #
 819         Jb = float(w.T.dot(w))
 820         Jo = float(dy.T.dot(dy))
 821         J  = Jo + Jb
 822         logging.debug("Values for cost functions are: J = %.5e  Jo = %.5e  Jb = %.5e"%(J,Jo,Jb))
 823         St["CurrentWeights"].append( w.squeeze() )
 824         St["CostFunctionJb"].append( Jb )
 825         St["CostFunctionJo"].append( Jo )
 826         St["CostFunctionJ"].append( J )
 827         if method == "Transform":
 828             GradJ = w - YL.T.dot(dy)
 829             HTild = IN + YL.T.dot(YL)
 830         else:
 831             GradJ = (nbMbr - 1)*w - YL.T.dot(RIdemi.dot(dy))
 832             HTild = (nbMbr - 1)*IN + YL.T.dot(RIdemi.dot(YL))
 833         HTild = numpy.array(HTild, dtype=float)
 834         dw = numpy.linalg.solve( HTild, numpy.array(GradJ, dtype=float) )
 835         w = w - dw
 836         j = j + 1
 837         convergence = _convergence_condition(j, dw, J, JPrev)
 838         JPrev = J
 839         #
 840         if method == "Transform":
 841             (U, s, _) = numpy.linalg.svd(HTild, full_matrices=False) # Hess = U s V
 842             T    = U.dot(numpy.diag(numpy.sqrt(1./s)).dot(U.T))   # T = Hess^(-1/2)
 843             Tinv = U.dot(numpy.diag(numpy.sqrt(s)).dot(U.T))      # Tinv = T^(-1)
 844     #
 845     # Stocke le dernier pas
 846     St["CurrentState"].append( numpy.mean( Ek, axis = 1).squeeze() )
 847     St["CurrentEnsemble"].append( Ek )
 848     #
 849     IndexMin = numpy.argmin( St["CostFunctionJ"][nbPS:] ) + nbPS
 850     xa = St["CurrentState"][IndexMin]
 851     Ea = St["CurrentEnsemble"][IndexMin]
 852     #
 853     return (xa, Ea, St)
 854
 855 def ienks(
 856         xb         = None,          # Background
 857         yObs       = None,          # Observation (série)
 858         E0         = None,          # Background ensemble
 859         B          = None,          # B
 860         RIdemi     = None,          # R^(-1/2)
 861         Mnnpu      = None,          # Evolution operator
 862         Hn         = None,          # Observation operator
 863         method     = "Transform",   # Bundle ou Transform
 864         nMembers   = 5,             # Number of members
 865         cMaximum   = 15000,         # Number of steps or cycles
 866         iMaximum   = 15000,         # Number of iterations per cycle
 867         sTolerance = mfp,           # Weights correction tolerance
 868         jTolerance = mfp,           # Cost decrement tolerance
 869         Lag        = 1,             # Lenght of smoothing window
 870         epsilon    = -1.,
 871         inflation  = 1.,
 872         nbPS       = 0,             # Number of previous steps
 873         setSeed    = None,
 874         ):
 875     #
 876     # Initial
 877     if setSeed is not None: numpy.random.seed(setSeed)
 878     if E0 is None: E0 = _BackgroundEnsembleGeneration( xb, B, nMembers)
 879     #
 880     # Spinup
 881     # ------
 882     #
 883     # Cycles
 884     # ------
 885     xa, Ea, Sa = [], [], []
 886     for i in range(Lag): # Lag void results
 887         xa.append([])
 888         Ea.append([])
 889         Sa.append([])
 890     for i in range(Lag,cMaximum):
 891         (xa_c, Ea_c, Sa_c) = _IEnKS_cycle_Lag_L_SDA_GN(
 892             E0,
 893             yObs[i-Lag:i],
 894             RIdemi,
 895             Mnnpu,
 896             Hn,
 897             method,
 898             iMaximum,
 899             sTolerance,
 900             jTolerance,
 901             Lag,
 902             epsilon,
 903             nbPS,
 904             )
 905         xa.append( xa_c )
 906         Ea.append( Ea_c )
 907         Sa.append( Sa_c )
 908         #
 909         # Inflation for next cycle
 910         E0 = xa_c + inflation * (Ea_c - xa_c)
 911     #
 912     return (xa, Ea, Sa)
 913
 914 # ==============================================================================
 915 if __name__ == "__main__":
 916     print('\n AUTODIAGNOSTIC\n')