# -*- coding: utf-8 -*-
#
-# Copyright (C) 2008-2021 EDF R&D
+# Copyright (C) 2008-2022 EDF R&D
#
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
"""
__author__ = "Jean-Philippe ARGAUD"
-import os, time, copy, types, sys, logging
-import math, numpy, scipy, scipy.optimize
-from daCore.BasicObjects import Operator
+import os, copy, types, sys, logging, numpy
+from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm
from daCore.PlatformInfo import PlatformInfo
mpr = PlatformInfo().MachinePrecision()
mfp = PlatformInfo().MaximumPrecision()
# logging.getLogger().setLevel(logging.DEBUG)
# ==============================================================================
-def ExecuteFunction( paire ):
- assert len(paire) == 2, "Incorrect number of arguments"
- X, funcrepr = paire
- __X = numpy.asmatrix(numpy.ravel( X )).T
+def ExecuteFunction( triplet ):
+ assert len(triplet) == 3, "Incorrect number of arguments"
+ X, xArgs, funcrepr = triplet
+ __X = numpy.ravel( X ).reshape((-1,1))
__sys_path_tmp = sys.path ; sys.path.insert(0,funcrepr["__userFunction__path"])
__module = __import__(funcrepr["__userFunction__modl"], globals(), locals(), [])
__fonction = getattr(__module,funcrepr["__userFunction__name"])
sys.path = __sys_path_tmp ; del __sys_path_tmp
- __HX = __fonction( __X )
+ if isinstance(xArgs, dict):
+ __HX = __fonction( __X, **xArgs )
+ else:
+ __HX = __fonction( __X )
return numpy.ravel( __HX )
# ==============================================================================
centeredDF = False,
increment = 0.01,
dX = None,
+ extraArguments = None,
+ reducingMemoryUse = False,
avoidingRedundancy = True,
toleranceInRedundancy = 1.e-18,
lenghtOfRedundancy = -1,
mfEnabled = False,
):
self.__name = str(name)
+ self.__extraArgs = extraArguments
+ #
if mpEnabled:
try:
import multiprocessing
self.__mpWorkers = None
logging.debug("FDA Calculs en multiprocessing : %s (nombre de processus : %s)"%(self.__mpEnabled,self.__mpWorkers))
#
- if mfEnabled:
- self.__mfEnabled = True
- else:
- self.__mfEnabled = False
+ self.__mfEnabled = bool(mfEnabled)
logging.debug("FDA Calculs en multifonctions : %s"%(self.__mfEnabled,))
#
+ self.__rmEnabled = bool(reducingMemoryUse)
+ logging.debug("FDA Calculs avec réduction mémoire : %s"%(self.__rmEnabled,))
+ #
if avoidingRedundancy:
self.__avoidRC = True
self.__tolerBP = float(toleranceInRedundancy)
self.__listJPIN = [] # Jacobian Previous Calculated Increment Norms
else:
self.__avoidRC = False
+ logging.debug("FDA Calculs avec réduction des doublons : %s"%self.__avoidRC)
+ if self.__avoidRC:
+ logging.debug("FDA Tolérance de détermination des doublons : %.2e"%self.__tolerBP)
#
if self.__mpEnabled:
if isinstance(Function,types.FunctionType):
self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
self.__userFunction__path = os.path.dirname(mod)
del mod
- self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
+ self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
elif isinstance(Function,types.MethodType):
logging.debug("FDA Calculs en multiprocessing : MethodType")
self.__userFunction__modl = os.path.basename(mod).replace('.pyc','').replace('.pyo','').replace('.py','')
self.__userFunction__path = os.path.dirname(mod)
del mod
- self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
+ self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
self.__userFunction = self.__userOperator.appliedTo # Pour le calcul Direct
else:
raise TypeError("User defined function or method has to be provided for finite differences approximation.")
else:
- self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled )
+ self.__userOperator = Operator( name = self.__name, fromMethod = Function, avoidingRedundancy = self.__avoidRC, inputAsMultiFunction = self.__mfEnabled, extraArguments = self.__extraArgs )
self.__userFunction = self.__userOperator.appliedTo
#
self.__centeredDF = bool(centeredDF)
if dX is None:
self.__dX = None
else:
- self.__dX = numpy.asmatrix(numpy.ravel( dX )).T
- logging.debug("FDA Reduction des doublons de calcul : %s"%self.__avoidRC)
- if self.__avoidRC:
- logging.debug("FDA Tolerance de determination des doublons : %.2e"%self.__tolerBP)
+ self.__dX = numpy.ravel( dX )
# ---------------------------------------------------------
def __doublon__(self, e, l, n, v=None):
return __ac, __iac
# ---------------------------------------------------------
- def DirectOperator(self, X ):
+ def __listdotwith__(self, __LMatrix, __dotWith = None, __dotTWith = None):
+ "Produit incrémental d'une matrice liste de colonnes avec un vecteur"
+ if not isinstance(__LMatrix, (list,tuple)):
+ raise TypeError("Columnwise list matrix has not the proper type: %s"%type(__LMatrix))
+ if __dotWith is not None:
+ __Idwx = numpy.ravel( __dotWith )
+ assert len(__LMatrix) == __Idwx.size, "Incorrect size of elements"
+ __Produit = numpy.zeros(__LMatrix[0].size)
+ for i, col in enumerate(__LMatrix):
+ __Produit += float(__Idwx[i]) * col
+ return __Produit
+ elif __dotTWith is not None:
+ _Idwy = numpy.ravel( __dotTWith ).T
+ assert __LMatrix[0].size == _Idwy.size, "Incorrect size of elements"
+ __Produit = numpy.zeros(len(__LMatrix))
+ for i, col in enumerate(__LMatrix):
+ __Produit[i] = float( _Idwy @ col)
+ return __Produit
+ else:
+ __Produit = None
+ return __Produit
+
+ # ---------------------------------------------------------
+ def DirectOperator(self, X, **extraArgs ):
"""
Calcul du direct à l'aide de la fonction fournie.
+
+ NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
+ ne doivent pas être données ici à la fonction utilisateur.
"""
logging.debug("FDA Calcul DirectOperator (explicite)")
if self.__mfEnabled:
_HX = self.__userFunction( X, argsAsSerie = True )
else:
- _X = numpy.asmatrix(numpy.ravel( X )).T
- _HX = numpy.ravel(self.__userFunction( _X ))
+ _HX = numpy.ravel(self.__userFunction( numpy.ravel(X) ))
#
return _HX
# ---------------------------------------------------------
- def TangentMatrix(self, X ):
+ def TangentMatrix(self, X, dotWith = None, dotTWith = None ):
"""
Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
c'est-à-dire le gradient de H en X. On utilise des différences finies
- directionnelles autour du point X. X est un numpy.matrix.
+ directionnelles autour du point X. X est un numpy.ndarray.
Différences finies centrées (approximation d'ordre 2):
1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
if X is None or len(X)==0:
raise ValueError("Nominal point X for approximate derivatives can not be None or void (given X: %s)."%(str(X),))
#
- _X = numpy.asmatrix(numpy.ravel( X )).T
+ _X = numpy.ravel( X )
#
if self.__dX is None:
_dX = self.__increment * _X
else:
- _dX = numpy.asmatrix(numpy.ravel( self.__dX )).T
+ _dX = numpy.ravel( self.__dX )
+ assert len(_X) == len(_dX), "Inconsistent dX increment length with respect to the X one"
+ assert _X.size == _dX.size, "Inconsistent dX increment size with respect to the X one"
#
if (_dX == 0.).any():
moyenne = _dX.mean()
__bidon, __alreadyCalculatedI = self.__doublon__(_dX, self.__listJPCI, self.__listJPIN, None)
if __alreadyCalculatedP == __alreadyCalculatedI > -1:
__alreadyCalculated, __i = True, __alreadyCalculatedP
- logging.debug("FDA Cas J déja calculé, récupération du doublon %i"%__i)
+ logging.debug("FDA Cas J déjà calculé, récupération du doublon %i"%__i)
#
if __alreadyCalculated:
logging.debug("FDA Calcul Jacobienne (par récupération du doublon %i)"%__i)
_Jacobienne = self.__listJPCR[__i]
+ logging.debug("FDA Fin du calcul de la Jacobienne")
+ if dotWith is not None:
+ return numpy.dot(_Jacobienne, numpy.ravel( dotWith ))
+ elif dotTWith is not None:
+ return numpy.dot(_Jacobienne.T, numpy.ravel( dotTWith ))
else:
logging.debug("FDA Calcul Jacobienne (explicite)")
if self.__centeredDF:
_jobs = []
for i in range( len(_dX) ):
_dXi = _dX[i]
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dXi
- _X_moins_dXi = numpy.array( _X.A1, dtype=float )
+ _X_moins_dXi = numpy.array( _X, dtype=float )
_X_moins_dXi[i] = _X[i] - _dXi
#
- _jobs.append( (_X_plus_dXi, funcrepr) )
- _jobs.append( (_X_moins_dXi, funcrepr) )
+ _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
+ _jobs.append( (_X_moins_dXi, self.__extraArgs, funcrepr) )
#
import multiprocessing
self.__pool = multiprocessing.Pool(self.__mpWorkers)
_xserie = []
for i in range( len(_dX) ):
_dXi = _dX[i]
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dXi
- _X_moins_dXi = numpy.array( _X.A1, dtype=float )
+ _X_moins_dXi = numpy.array( _X, dtype=float )
_X_moins_dXi[i] = _X[i] - _dXi
#
_xserie.append( _X_plus_dXi )
_Jacobienne = []
for i in range( _dX.size ):
_dXi = _dX[i]
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dXi
- _X_moins_dXi = numpy.array( _X.A1, dtype=float )
+ _X_moins_dXi = numpy.array( _X, dtype=float )
_X_moins_dXi[i] = _X[i] - _dXi
#
_HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
"__userFunction__name" : self.__userFunction__name,
}
_jobs = []
- _jobs.append( (_X.A1, funcrepr) )
+ _jobs.append( (_X, self.__extraArgs, funcrepr) )
for i in range( len(_dX) ):
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dX[i]
#
- _jobs.append( (_X_plus_dXi, funcrepr) )
+ _jobs.append( (_X_plus_dXi, self.__extraArgs, funcrepr) )
#
import multiprocessing
self.__pool = multiprocessing.Pool(self.__mpWorkers)
#
elif self.__mfEnabled:
_xserie = []
- _xserie.append( _X.A1 )
+ _xserie.append( _X )
for i in range( len(_dX) ):
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dX[i]
#
_xserie.append( _X_plus_dXi )
_HX = self.DirectOperator( _X )
for i in range( _dX.size ):
_dXi = _dX[i]
- _X_plus_dXi = numpy.array( _X.A1, dtype=float )
+ _X_plus_dXi = numpy.array( _X, dtype=float )
_X_plus_dXi[i] = _X[i] + _dXi
#
_HX_plus_dXi = self.DirectOperator( _X_plus_dXi )
#
_Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) )
- #
#
- _Jacobienne = numpy.asmatrix( numpy.vstack( _Jacobienne ) ).T
- if self.__avoidRC:
- if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
- while len(self.__listJPCP) > self.__lenghtRJ:
- self.__listJPCP.pop(0)
- self.__listJPCI.pop(0)
- self.__listJPCR.pop(0)
- self.__listJPPN.pop(0)
- self.__listJPIN.pop(0)
- self.__listJPCP.append( copy.copy(_X) )
- self.__listJPCI.append( copy.copy(_dX) )
- self.__listJPCR.append( copy.copy(_Jacobienne) )
- self.__listJPPN.append( numpy.linalg.norm(_X) )
- self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
- #
- logging.debug("FDA Fin du calcul de la Jacobienne")
+ if (dotWith is not None) or (dotTWith is not None):
+ __Produit = self.__listdotwith__(_Jacobienne, dotWith, dotTWith)
+ else:
+ __Produit = None
+ if __Produit is None or self.__avoidRC:
+ _Jacobienne = numpy.transpose( numpy.vstack( _Jacobienne ) )
+ if self.__avoidRC:
+ if self.__lenghtRJ < 0: self.__lenghtRJ = 2 * _X.size
+ while len(self.__listJPCP) > self.__lenghtRJ:
+ self.__listJPCP.pop(0)
+ self.__listJPCI.pop(0)
+ self.__listJPCR.pop(0)
+ self.__listJPPN.pop(0)
+ self.__listJPIN.pop(0)
+ self.__listJPCP.append( copy.copy(_X) )
+ self.__listJPCI.append( copy.copy(_dX) )
+ self.__listJPCR.append( copy.copy(_Jacobienne) )
+ self.__listJPPN.append( numpy.linalg.norm(_X) )
+ self.__listJPIN.append( numpy.linalg.norm(_Jacobienne) )
+ logging.debug("FDA Fin du calcul de la Jacobienne")
+ if __Produit is not None:
+ return __Produit
#
return _Jacobienne
# ---------------------------------------------------------
- def TangentOperator(self, paire ):
+ def TangentOperator(self, paire, **extraArgs ):
"""
Calcul du tangent à l'aide de la Jacobienne.
+
+ NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
+ ne doivent pas être données ici à la fonction utilisateur.
"""
if self.__mfEnabled:
- assert len(paire) == 1, "Incorrect lenght of arguments"
+ assert len(paire) == 1, "Incorrect length of arguments"
_paire = paire[0]
assert len(_paire) == 2, "Incorrect number of arguments"
else:
assert len(paire) == 2, "Incorrect number of arguments"
_paire = paire
X, dX = _paire
- _Jacobienne = self.TangentMatrix( X )
if dX is None or len(dX) == 0:
#
# Calcul de la forme matricielle si le second argument est None
# -------------------------------------------------------------
+ _Jacobienne = self.TangentMatrix( X )
if self.__mfEnabled: return [_Jacobienne,]
else: return _Jacobienne
else:
#
# Calcul de la valeur linéarisée de H en X appliqué à dX
# ------------------------------------------------------
- _dX = numpy.asmatrix(numpy.ravel( dX )).T
- _HtX = numpy.dot(_Jacobienne, _dX)
- if self.__mfEnabled: return [_HtX.A1,]
- else: return _HtX.A1
+ _HtX = self.TangentMatrix( X, dotWith = dX )
+ if self.__mfEnabled: return [_HtX,]
+ else: return _HtX
# ---------------------------------------------------------
- def AdjointOperator(self, paire ):
+ def AdjointOperator(self, paire, **extraArgs ):
"""
Calcul de l'adjoint à l'aide de la Jacobienne.
+
+ NB : les extraArgs sont là pour assurer la compatibilité d'appel, mais
+ ne doivent pas être données ici à la fonction utilisateur.
"""
if self.__mfEnabled:
- assert len(paire) == 1, "Incorrect lenght of arguments"
+ assert len(paire) == 1, "Incorrect length of arguments"
_paire = paire[0]
assert len(_paire) == 2, "Incorrect number of arguments"
else:
assert len(paire) == 2, "Incorrect number of arguments"
_paire = paire
X, Y = _paire
- _JacobienneT = self.TangentMatrix( X ).T
if Y is None or len(Y) == 0:
#
# Calcul de la forme matricielle si le second argument est None
# -------------------------------------------------------------
+ _JacobienneT = self.TangentMatrix( X ).T
if self.__mfEnabled: return [_JacobienneT,]
else: return _JacobienneT
else:
#
# Calcul de la valeur de l'adjoint en X appliqué à Y
# --------------------------------------------------
- _Y = numpy.asmatrix(numpy.ravel( Y )).T
- _HaY = numpy.dot(_JacobienneT, _Y)
- if self.__mfEnabled: return [_HaY.A1,]
- else: return _HaY.A1
+ _HaY = self.TangentMatrix( X, dotTWith = Y )
+ if self.__mfEnabled: return [_HaY,]
+ else: return _HaY
# ==============================================================================
-def mmqr(
- func = None,
- x0 = None,
- fprime = None,
- bounds = None,
- quantile = 0.5,
- maxfun = 15000,
- toler = 1.e-06,
- y = None,
- ):
- """
- Implémentation informatique de l'algorithme MMQR, basée sur la publication :
- David R. Hunter, Kenneth Lange, "Quantile Regression via an MM Algorithm",
- Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000.
- """
+def SetInitialDirection( __Direction = [], __Amplitude = 1., __Position = None ):
+ "Établit ou élabore une direction avec une amplitude"
#
- # Recuperation des donnees et informations initiales
- # --------------------------------------------------
- variables = numpy.ravel( x0 )
- mesures = numpy.ravel( y )
- increment = sys.float_info[0]
- p = variables.size
- n = mesures.size
- quantile = float(quantile)
+ if len(__Direction) == 0 and __Position is None:
+ raise ValueError("If initial direction is void, current position has to be given")
+ if abs(float(__Amplitude)) < mpr:
+ raise ValueError("Amplitude of perturbation can not be zero")
#
- # Calcul des parametres du MM
- # ---------------------------
- tn = float(toler) / n
- e0 = -tn / math.log(tn)
- epsilon = (e0-tn)/(1+math.log(e0))
+ if len(__Direction) > 0:
+ __dX0 = numpy.asarray(__Direction)
+ else:
+ __dX0 = []
+ __X0 = numpy.ravel(numpy.asarray(__Position))
+ __mX0 = numpy.mean( __X0, dtype=mfp )
+ if abs(__mX0) < 2*mpr: __mX0 = 1. # Évite le problème de position nulle
+ for v in __X0:
+ if abs(v) > 1.e-8:
+ __dX0.append( numpy.random.normal(0.,abs(v)) )
+ else:
+ __dX0.append( numpy.random.normal(0.,__mX0) )
#
- # Calculs d'initialisation
- # ------------------------
- residus = mesures - numpy.ravel( func( variables ) )
- poids = 1./(epsilon+numpy.abs(residus))
- veps = 1. - 2. * quantile - residus * poids
- lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus)
- iteration = 0
+ __dX0 = numpy.asarray(__dX0,float) # Évite le problème d'array de taille 1
+ __dX0 = numpy.ravel( __dX0 ) # Redresse les vecteurs
+ __dX0 = float(__Amplitude) * __dX0
#
- # Recherche iterative
- # -------------------
- while (increment > toler) and (iteration < maxfun) :
- iteration += 1
- #
- Derivees = numpy.array(fprime(variables))
- Derivees = Derivees.reshape(n,p) # Necessaire pour remettre en place la matrice si elle passe par des tuyaux YACS
- DeriveesT = Derivees.transpose()
- M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) )
- SM = numpy.transpose(numpy.dot( DeriveesT , veps ))
- step = - numpy.linalg.lstsq( M, SM, rcond=-1 )[0]
- #
- variables = variables + step
- if bounds is not None:
- # Attention : boucle infinie à éviter si un intervalle est trop petit
- while( (variables < numpy.ravel(numpy.asmatrix(bounds)[:,0])).any() or (variables > numpy.ravel(numpy.asmatrix(bounds)[:,1])).any() ):
- step = step/2.
- variables = variables - step
- residus = mesures - numpy.ravel( func(variables) )
- surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
- #
- while ( (surrogate > lastsurrogate) and ( max(list(numpy.abs(step))) > 1.e-16 ) ) :
- step = step/2.
- variables = variables - step
- residus = mesures - numpy.ravel( func(variables) )
- surrogate = numpy.sum(residus**2 * poids) + (4.*quantile-2.) * numpy.sum(residus)
- #
- increment = lastsurrogate-surrogate
- poids = 1./(epsilon+numpy.abs(residus))
- veps = 1. - 2. * quantile - residus * poids
- lastsurrogate = -numpy.sum(residus * veps) - (1.-2.*quantile)*numpy.sum(residus)
+ return __dX0
+
+# ==============================================================================
+def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
+ "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
#
- # Mesure d'écart
- # --------------
- Ecart = quantile * numpy.sum(residus) - numpy.sum( residus[residus<0] )
+ __bgCenter = numpy.ravel(__bgCenter)[:,None]
+ if __nbMembers < 1:
+ raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
+ #
+ if __bgCovariance is None:
+ _Perturbations = numpy.tile( __bgCenter, __nbMembers)
+ else:
+ _Z = numpy.random.multivariate_normal(numpy.zeros(__bgCenter.size), __bgCovariance, size=__nbMembers).T
+ _Perturbations = numpy.tile( __bgCenter, __nbMembers) + _Z
+ #
+ return _Perturbations
+
+# ==============================================================================
+def EnsembleOfBackgroundPerturbations( __bgCenter, __bgCovariance, __nbMembers, __withSVD = True):
+ "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"
+ def __CenteredRandomAnomalies(Zr, N):
+ """
+ Génère une matrice de N anomalies aléatoires centrées sur Zr selon les
+ notes manuscrites de MB et conforme au code de PS avec eps = -1
+ """
+ eps = -1
+ Q = numpy.identity(N-1)-numpy.ones((N-1,N-1))/numpy.sqrt(N)/(numpy.sqrt(N)-eps)
+ Q = numpy.concatenate((Q, [eps*numpy.ones(N-1)/numpy.sqrt(N)]), axis=0)
+ R, _ = numpy.linalg.qr(numpy.random.normal(size = (N-1,N-1)))
+ Q = numpy.dot(Q,R)
+ Zr = numpy.dot(Q,Zr)
+ return Zr.T
+ #
+ __bgCenter = numpy.ravel(__bgCenter).reshape((-1,1))
+ if __nbMembers < 1:
+ raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(__nbMembers),))
+ if __bgCovariance is None:
+ _Perturbations = numpy.tile( __bgCenter, __nbMembers)
+ else:
+ if __withSVD:
+ _U, _s, _V = numpy.linalg.svd(__bgCovariance, full_matrices=False)
+ _nbctl = __bgCenter.size
+ if __nbMembers > _nbctl:
+ _Z = numpy.concatenate((numpy.dot(
+ numpy.diag(numpy.sqrt(_s[:_nbctl])), _V[:_nbctl]),
+ numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1-_nbctl)), axis = 0)
+ else:
+ _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:__nbMembers-1])), _V[:__nbMembers-1])
+ _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
+ _Perturbations = __bgCenter + _Zca
+ else:
+ if max(abs(__bgCovariance.flatten())) > 0:
+ _nbctl = __bgCenter.size
+ _Z = numpy.random.multivariate_normal(numpy.zeros(_nbctl),__bgCovariance,__nbMembers-1)
+ _Zca = __CenteredRandomAnomalies(_Z, __nbMembers)
+ _Perturbations = __bgCenter + _Zca
+ else:
+ _Perturbations = numpy.tile( __bgCenter, __nbMembers)
+ #
+ return _Perturbations
+
+# ==============================================================================
+def EnsembleMean( __Ensemble ):
+ "Renvoie la moyenne empirique d'un ensemble"
+ return numpy.asarray(__Ensemble).mean(axis=1, dtype=mfp).astype('float').reshape((-1,1))
+
+# ==============================================================================
+def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.):
+ "Renvoie les anomalies centrées à partir d'un ensemble"
+ if __OptMean is None:
+ __Em = EnsembleMean( __Ensemble )
+ else:
+ __Em = numpy.ravel( __OptMean ).reshape((-1,1))
+ #
+ return __Normalisation * (numpy.asarray( __Ensemble ) - __Em)
+
+# ==============================================================================
+def EnsembleErrorCovariance( __Ensemble, __Quick = False ):
+ "Renvoie l'estimation empirique de la covariance d'ensemble"
+ if __Quick:
+ # Covariance rapide mais rarement définie positive
+ __Covariance = numpy.cov( __Ensemble )
+ else:
+ # Résultat souvent identique à numpy.cov, mais plus robuste
+ __n, __m = numpy.asarray( __Ensemble ).shape
+ __Anomalies = EnsembleOfAnomalies( __Ensemble )
+ # Estimation empirique
+ __Covariance = ( __Anomalies @ __Anomalies.T ) / (__m-1)
+ # Assure la symétrie
+ __Covariance = ( __Covariance + __Covariance.T ) * 0.5
+ # Assure la positivité
+ __epsilon = mpr*numpy.trace( __Covariance )
+ __Covariance = __Covariance + __epsilon * numpy.identity(__n)
+ #
+ return __Covariance
+
+# ==============================================================================
+def EnsemblePerturbationWithGivenCovariance(
+ __Ensemble,
+ __Covariance,
+ __Seed = None,
+ ):
+ "Ajout d'une perturbation à chaque membre d'un ensemble selon une covariance prescrite"
+ if hasattr(__Covariance,"assparsematrix"):
+ if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance.assparsematrix())/abs(__Ensemble).mean() < mpr).all():
+ # Traitement d'une covariance nulle ou presque
+ return __Ensemble
+ if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance.assparsematrix()) < mpr).all():
+ # Traitement d'une covariance nulle ou presque
+ return __Ensemble
+ else:
+ if (abs(__Ensemble).mean() > mpr) and (abs(__Covariance)/abs(__Ensemble).mean() < mpr).all():
+ # Traitement d'une covariance nulle ou presque
+ return __Ensemble
+ if (abs(__Ensemble).mean() <= mpr) and (abs(__Covariance) < mpr).all():
+ # Traitement d'une covariance nulle ou presque
+ return __Ensemble
+ #
+ __n, __m = __Ensemble.shape
+ if __Seed is not None: numpy.random.seed(__Seed)
+ #
+ if hasattr(__Covariance,"isscalar") and __Covariance.isscalar():
+ # Traitement d'une covariance multiple de l'identité
+ __zero = 0.
+ __std = numpy.sqrt(__Covariance.assparsematrix())
+ __Ensemble += numpy.random.normal(__zero, __std, size=(__m,__n)).T
+ #
+ elif hasattr(__Covariance,"isvector") and __Covariance.isvector():
+ # Traitement d'une covariance diagonale avec variances non identiques
+ __zero = numpy.zeros(__n)
+ __std = numpy.sqrt(__Covariance.assparsematrix())
+ __Ensemble += numpy.asarray([numpy.random.normal(__zero, __std) for i in range(__m)]).T
+ #
+ elif hasattr(__Covariance,"ismatrix") and __Covariance.ismatrix():
+ # Traitement d'une covariance pleine
+ __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance.asfullmatrix(__n), size=__m).T
+ #
+ elif isinstance(__Covariance, numpy.ndarray):
+ # Traitement d'une covariance numpy pleine, sachant qu'on arrive ici en dernier
+ __Ensemble += numpy.random.multivariate_normal(numpy.zeros(__n), __Covariance, size=__m).T
+ #
+ else:
+ raise ValueError("Error in ensemble perturbation with inadequate covariance specification")
#
- return variables, Ecart, [n,p,iteration,increment,0]
+ return __Ensemble
# ==============================================================================
def CovarianceInflation(
- InputCovOrEns,
- InflationType = None,
- InflationFactor = None,
- BackgroundCov = None,
+ __InputCovOrEns,
+ __InflationType = None,
+ __InflationFactor = None,
+ __BackgroundCov = None,
):
"""
Inflation applicable soit sur Pb ou Pa, soit sur les ensembles EXb ou EXa
Synthèse : Hunt 2007, section 2.3.5
"""
- if InflationFactor is None:
- return InputCovOrEns
+ if __InflationFactor is None:
+ return __InputCovOrEns
else:
- InflationFactor = float(InflationFactor)
+ __InflationFactor = float(__InflationFactor)
#
- if InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
- if InflationFactor < 1.:
- raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
- if InflationFactor < 1.+mpr:
- return InputCovOrEns
- OutputCovOrEns = InflationFactor**2 * InputCovOrEns
+ __InputCovOrEns = numpy.asarray(__InputCovOrEns)
+ if __InputCovOrEns.size == 0: return __InputCovOrEns
#
- elif InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
- if InflationFactor < 1.:
+ if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]:
+ if __InflationFactor < 1.:
raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
- if InflationFactor < 1.+mpr:
- return InputCovOrEns
- InputCovOrEnsMean = InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
- OutputCovOrEns = InputCovOrEnsMean[:,numpy.newaxis] \
- + InflationFactor * (InputCovOrEns - InputCovOrEnsMean[:,numpy.newaxis])
+ if __InflationFactor < 1.+mpr: # No inflation = 1
+ return __InputCovOrEns
+ __OutputCovOrEns = __InflationFactor**2 * __InputCovOrEns
#
- elif InflationType in ["AdditiveOnBackgroundCovariance", "AdditiveOnAnalysisCovariance"]:
- if InflationFactor < 0.:
+ elif __InflationType in ["MultiplicativeOnAnalysisAnomalies", "MultiplicativeOnBackgroundAnomalies"]:
+ if __InflationFactor < 1.:
+ raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.")
+ if __InflationFactor < 1.+mpr: # No inflation = 1
+ return __InputCovOrEns
+ __InputCovOrEnsMean = __InputCovOrEns.mean(axis=1, dtype=mfp).astype('float')
+ __OutputCovOrEns = __InputCovOrEnsMean[:,numpy.newaxis] \
+ + __InflationFactor * (__InputCovOrEns - __InputCovOrEnsMean[:,numpy.newaxis])
+ #
+ elif __InflationType in ["AdditiveOnAnalysisCovariance", "AdditiveOnBackgroundCovariance"]:
+ if __InflationFactor < 0.:
raise ValueError("Inflation factor for additive inflation has to be greater or equal than 0.")
- if InflationFactor < mpr:
- return InputCovOrEns
- __n, __m = numpy.asarray(InputCovOrEns).shape
+ if __InflationFactor < mpr: # No inflation = 0
+ return __InputCovOrEns
+ __n, __m = __InputCovOrEns.shape
if __n != __m:
raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
- OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * numpy.eye(__n)
+ __tr = __InputCovOrEns.trace()/__n
+ if __InflationFactor > __tr:
+ raise ValueError("Inflation factor for additive inflation has to be small over %.0e."%__tr)
+ __OutputCovOrEns = (1. - __InflationFactor)*__InputCovOrEns + __InflationFactor * numpy.identity(__n)
#
- elif InflationType == "HybridOnBackgroundCovariance":
- if InflationFactor < 0.:
+ elif __InflationType == "HybridOnBackgroundCovariance":
+ if __InflationFactor < 0.:
raise ValueError("Inflation factor for hybrid inflation has to be greater or equal than 0.")
- if InflationFactor < mpr:
- return InputCovOrEns
- __n, __m = numpy.asarray(InputCovOrEns).shape
+ if __InflationFactor < mpr: # No inflation = 0
+ return __InputCovOrEns
+ __n, __m = __InputCovOrEns.shape
if __n != __m:
raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.")
- if BackgroundCov is None:
+ if __BackgroundCov is None:
raise ValueError("Background covariance matrix B has to be given for hybrid inflation.")
- if InputCovOrEns.shape != BackgroundCov.shape:
+ if __InputCovOrEns.shape != __BackgroundCov.shape:
raise ValueError("Ensemble covariance matrix has to be of same size than background covariance matrix B.")
- OutputCovOrEns = (1. - InflationFactor) * InputCovOrEns + InflationFactor * BackgroundCov
+ __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov
#
- elif InflationType == "Relaxation":
- raise NotImplementedError("InflationType Relaxation")
+ elif __InflationType == "Relaxation":
+ raise NotImplementedError("Relaxation inflation type not implemented")
#
else:
- raise ValueError("Error in inflation type, '%s' is not a valid keyword."%InflationType)
+ raise ValueError("Error in inflation type, '%s' is not a valid keyword."%__InflationType)
#
- return OutputCovOrEns
+ return __OutputCovOrEns
# ==============================================================================
-def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q):
- """
- Stochastic EnKF (Envensen 1994, Burgers 1998)
+def HessienneEstimation(__selfA, __nb, __HaM, __HtM, __BI, __RI):
+ "Estimation de la Hessienne"
+ #
+ __HessienneI = []
+ for i in range(int(__nb)):
+ __ee = numpy.zeros((__nb,1))
+ __ee[i] = 1.
+ __HtEE = numpy.dot(__HtM,__ee).reshape((-1,1))
+ __HessienneI.append( numpy.ravel( __BI * __ee + __HaM * (__RI * __HtEE) ) )
+ #
+ __A = numpy.linalg.inv(numpy.array( __HessienneI ))
+ __A = (__A + __A.T) * 0.5 # Symétrie
+ __A = __A + mpr*numpy.trace( __A ) * numpy.identity(__nb) # Positivité
+ #
+ if min(__A.shape) != max(__A.shape):
+ raise ValueError("The %s a posteriori covariance matrix A is of shape %s, despites it has to be a squared matrix. There is an error in the observation operator, please check it."%(__selfA._name,str(__A.shape)))
+ if (numpy.diag(__A) < 0).any():
+ raise ValueError("The %s a posteriori covariance matrix A has at least one negative value on its diagonal. There is an error in the observation operator, please check it."%(__selfA._name,))
+ if logging.getLogger().level < logging.WARNING: # La vérification n'a lieu qu'en debug
+ try:
+ numpy.linalg.cholesky( __A )
+ except:
+ raise ValueError("The %s a posteriori covariance matrix A is not symmetric positive-definite. Please check your a priori covariances and your observation operator."%(__selfA._name,))
+ #
+ return __A
- selfA est identique au "self" d'algorithme appelant et contient les
- valeurs.
- """
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA._parameters["StoreInternalVariables"] = True
- #
- # Opérateurs
- # ----------
- H = HO["Direct"].appliedControledFormTo
- #
- if selfA._parameters["EstimationOf"] == "State":
- M = EM["Direct"].appliedControledFormTo
- #
- if CM is not None and "Tangent" in CM and U is not None:
- Cm = CM["Tangent"].asMatrix(Xb)
- else:
- Cm = None
- #
- # Nombre de pas identique au nombre de pas d'observations
- # -------------------------------------------------------
- if hasattr(Y,"stepnumber"):
- duration = Y.stepnumber()
- __p = numpy.cumprod(Y.shape())[-1]
+# ==============================================================================
+def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None):
+ "Estimation des quantiles a posteriori à partir de A>0 (selfA est modifié)"
+ nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"]
+ #
+ # Traitement des bornes
+ if "StateBoundsForQuantiles" in selfA._parameters:
+ LBounds = selfA._parameters["StateBoundsForQuantiles"] # Prioritaire
+ elif "Bounds" in selfA._parameters:
+ LBounds = selfA._parameters["Bounds"] # Défaut raisonnable
else:
- duration = 2
- __p = numpy.array(Y).size
- #
- # Précalcul des inversions de B et R
- # ----------------------------------
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- BI = B.getI()
- RI = R.getI()
- #
- # Initialisation
- # --------------
- __n = Xb.size
- __m = selfA._parameters["NumberOfMembers"]
- Xn = numpy.asmatrix(numpy.dot( Xb.reshape(__n,1), numpy.ones((1,__m)) ))
- if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
- else: Pn = B
- if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
- else: Rn = R
- if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
- else: Qn = Q
- #
- if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
- selfA.StoredVariables["Analysis"].store( numpy.ravel(Xb) )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- covarianceXa = Pn
- #
- previousJMinimum = numpy.finfo(float).max
- #
- # Predimensionnement
- Xn_predicted = numpy.asmatrix(numpy.zeros((__n,__m)))
- HX_predicted = numpy.asmatrix(numpy.zeros((__p,__m)))
- #
- for step in range(duration-1):
- if hasattr(Y,"store"):
- Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
+ LBounds = None
+ if LBounds is not None:
+ LBounds = ForceNumericBounds( LBounds )
+ __Xa = numpy.ravel(Xa)
+ #
+ # Échantillonnage des états
+ YfQ = None
+ EXr = None
+ for i in range(nbsamples):
+ if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None and HXa is not None:
+ dXr = (numpy.random.multivariate_normal(__Xa,A) - __Xa).reshape((-1,1))
+ if LBounds is not None: # "EstimateProjection" par défaut
+ dXr = numpy.max(numpy.hstack((dXr,LBounds[:,0].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
+ dXr = numpy.min(numpy.hstack((dXr,LBounds[:,1].reshape((-1,1))) - __Xa.reshape((-1,1))),axis=1)
+ dYr = HtM @ dXr
+ Yr = HXa.reshape((-1,1)) + dYr
+ if selfA._toStore("SampledStateForQuantiles"): Xr = __Xa + numpy.ravel(dXr)
+ elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None:
+ Xr = numpy.random.multivariate_normal(__Xa,A)
+ if LBounds is not None: # "EstimateProjection" par défaut
+ Xr = numpy.max(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,0].reshape((-1,1)))),axis=1)
+ Xr = numpy.min(numpy.hstack((Xr.reshape((-1,1)),LBounds[:,1].reshape((-1,1)))),axis=1)
+ Yr = numpy.asarray(Hm( Xr ))
else:
- Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
+ raise ValueError("Quantile simulations has only to be Linear or NonLinear.")
#
- if U is not None:
- if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
- elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
- else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ if YfQ is None:
+ YfQ = Yr.reshape((-1,1))
+ if selfA._toStore("SampledStateForQuantiles"): EXr = Xr.reshape((-1,1))
else:
- Un = None
- #
- if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
- Xn = CovarianceInflation( Xn,
- selfA._parameters["InflationType"],
- selfA._parameters["InflationFactor"],
- )
- #
- if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
- for i in range(__m):
- qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn, (1,1,1))).T
- Xn_predicted[:,i] = numpy.asmatrix(numpy.ravel( M((Xn[:,i], Un)) )).T + qi
- HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( H((Xn_predicted[:,i], Un)) )).T
- if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
- Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
- elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
- # --- > Par principe, M = Id, Q = 0
- Xn_predicted = Xn
- for i in range(__m):
- HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( H((Xn_predicted[:,i], Un)) )).T
- #
- # Mean of forecast and observation of forecast
- Xfm = numpy.asmatrix(numpy.ravel(Xn_predicted.mean(axis=1, dtype=mfp).astype('float'))).T
- Hfm = numpy.asmatrix(numpy.ravel(HX_predicted.mean(axis=1, dtype=mfp).astype('float'))).T
- #
- PfHT, HPfHT = 0., 0.
- for i in range(__m):
- Exfi = Xn_predicted[:,i] - Xfm
- Eyfi = HX_predicted[:,i] - Hfm
- PfHT += Exfi * Eyfi.T
- HPfHT += Eyfi * Eyfi.T
- PfHT = (1./(__m-1)) * PfHT
- HPfHT = (1./(__m-1)) * HPfHT
- K = PfHT * ( R + HPfHT ).I
- del PfHT, HPfHT
- #
- for i in range(__m):
- ri = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__p), Rn, (1,1,1))).T
- Xn[:,i] = Xn_predicted[:,i] + K * (Ynpu + ri - HX_predicted[:,i])
- #
- if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
- Xn = CovarianceInflation( Xn,
- selfA._parameters["InflationType"],
- selfA._parameters["InflationFactor"],
- )
- #
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float')
- #
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("APosterioriCovariance") \
- or selfA._toStore("InnovationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
- _Innovation = Ynpu - _HXa
- #
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- # ---> avec analysis
- selfA.StoredVariables["Analysis"].store( Xa )
- if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
- selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
- if selfA._toStore("InnovationAtCurrentAnalysis"):
- selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
- # ---> avec current state
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CurrentState"):
- selfA.StoredVariables["CurrentState"].store( Xn )
- if selfA._toStore("ForecastState"):
- selfA.StoredVariables["ForecastState"].store( Xn_predicted )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
- if selfA._toStore("InnovationAtCurrentState"):
- selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
- if selfA._toStore("SimulatedObservationAtCurrentState") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
- # ---> autres
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
- J = Jb + Jo
- selfA.StoredVariables["CostFunctionJb"].store( Jb )
- selfA.StoredVariables["CostFunctionJo"].store( Jo )
- selfA.StoredVariables["CostFunctionJ" ].store( J )
- #
- if selfA._toStore("IndexOfOptimum") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("CostFunctionJAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
- if selfA._toStore("IndexOfOptimum"):
- selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
- if selfA._toStore("CurrentOptimum"):
- selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
- if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
- if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
- if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
- if selfA._toStore("CostFunctionJAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
- if selfA._toStore("APosterioriCovariance"):
- Eai = (1/numpy.sqrt(__m-1)) * (Xn - Xa.reshape((__n,-1))) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- if selfA._parameters["EstimationOf"] == "Parameters" \
- and J < previousJMinimum:
- previousJMinimum = J
- XaMin = Xa
- if selfA._toStore("APosterioriCovariance"):
- covarianceXaMin = Pn
- #
- # Stockage final supplémentaire de l'optimum en estimation de paramètres
- # ----------------------------------------------------------------------
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- selfA.StoredVariables["Analysis"].store( XaMin )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
+ YfQ = numpy.hstack((YfQ,Yr.reshape((-1,1))))
+ if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.hstack((EXr,Xr.reshape((-1,1))))
+ #
+ # Extraction des quantiles
+ YfQ.sort(axis=-1)
+ YQ = None
+ for quantile in selfA._parameters["Quantiles"]:
+ if not (0. <= float(quantile) <= 1.): continue
+ indice = int(nbsamples * float(quantile) - 1./nbsamples)
+ if YQ is None: YQ = YfQ[:,indice].reshape((-1,1))
+ else: YQ = numpy.hstack((YQ,YfQ[:,indice].reshape((-1,1))))
+ if YQ is not None: # Liste non vide de quantiles
+ selfA.StoredVariables["SimulationQuantiles"].store( YQ )
+ if selfA._toStore("SampledStateForQuantiles"):
+ selfA.StoredVariables["SampledStateForQuantiles"].store( EXr )
#
return 0
# ==============================================================================
-def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, KorV="KalmanFilterFormula"):
- """
- Ensemble-Transform EnKF (ETKF or Deterministic EnKF: Bishop 2001, Hunt 2007)
+def ForceNumericBounds( __Bounds ):
+ "Force les bornes à être des valeurs numériques, sauf si globalement None"
+ # Conserve une valeur par défaut à None s'il n'y a pas de bornes
+ if __Bounds is None: return None
+ # Converti toutes les bornes individuelles None à +/- l'infini
+ __Bounds = numpy.asarray( __Bounds, dtype=float )
+ if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
+ raise ValueError("Incorrectly shaped bounds data")
+ __Bounds[numpy.isnan(__Bounds[:,0]),0] = -sys.float_info.max
+ __Bounds[numpy.isnan(__Bounds[:,1]),1] = sys.float_info.max
+ return __Bounds
- selfA est identique au "self" d'algorithme appelant et contient les
- valeurs.
- """
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA._parameters["StoreInternalVariables"] = True
- #
- # Opérateurs
- # ----------
- H = HO["Direct"].appliedControledFormTo
- #
- if selfA._parameters["EstimationOf"] == "State":
- M = EM["Direct"].appliedControledFormTo
- #
- if CM is not None and "Tangent" in CM and U is not None:
- Cm = CM["Tangent"].asMatrix(Xb)
+# ==============================================================================
+def RecentredBounds( __Bounds, __Center, __Scale = None):
+ "Recentre les bornes autour de 0, sauf si globalement None"
+ # Conserve une valeur par défaut à None s'il n'y a pas de bornes
+ if __Bounds is None: return None
+ if __Scale is None:
+ # Recentre les valeurs numériques de bornes
+ return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1))
else:
- Cm = None
- #
- # Nombre de pas identique au nombre de pas d'observations
- # -------------------------------------------------------
- if hasattr(Y,"stepnumber"):
- duration = Y.stepnumber()
- __p = numpy.cumprod(Y.shape())[-1]
+ # Recentre les valeurs numériques de bornes et change l'échelle par une matrice
+ return __Scale @ (ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1)))
+
+# ==============================================================================
+def ApplyBounds( __Vector, __Bounds, __newClip = True):
+ "Applique des bornes numériques à un point"
+ # Conserve une valeur par défaut s'il n'y a pas de bornes
+ if __Bounds is None: return __Vector
+ #
+ if not isinstance(__Vector, numpy.ndarray): # Is an array
+ raise ValueError("Incorrect array definition of vector data")
+ if not isinstance(__Bounds, numpy.ndarray): # Is an array
+ raise ValueError("Incorrect array definition of bounds data")
+ if 2*__Vector.size != __Bounds.size: # Is a 2 column array of vector lenght
+ raise ValueError("Incorrect bounds number (%i) to be applied for this vector (of size %i)"%(__Bounds.size,__Vector.size))
+ if len(__Bounds.shape) != 2 or min(__Bounds.shape) <= 0 or __Bounds.shape[1] != 2:
+ raise ValueError("Incorrectly shaped bounds data")
+ #
+ if __newClip:
+ __Vector = __Vector.clip(
+ __Bounds[:,0].reshape(__Vector.shape),
+ __Bounds[:,1].reshape(__Vector.shape),
+ )
else:
- duration = 2
- __p = numpy.array(Y).size
- #
- # Précalcul des inversions de B et R
- # ----------------------------------
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- BI = B.getI()
- RI = R.getI()
- elif KorV != "KalmanFilterFormula":
- RI = R.getI()
- if KorV == "KalmanFilterFormula":
- RIdemi = R.choleskyI()
+ __Vector = numpy.max(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,0])),axis=1)
+ __Vector = numpy.min(numpy.hstack((__Vector.reshape((-1,1)),numpy.asmatrix(__Bounds)[:,1])),axis=1)
+ __Vector = numpy.asarray(__Vector)
#
- # Initialisation
- # --------------
- __n = Xb.size
- __m = selfA._parameters["NumberOfMembers"]
- Xn = numpy.asmatrix(numpy.dot( Xb.reshape(__n,1), numpy.ones((1,__m)) ))
- if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
- else: Pn = B
- if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
- else: Rn = R
- if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
- else: Qn = Q
- #
- if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
- selfA.StoredVariables["Analysis"].store( numpy.ravel(Xb) )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- covarianceXa = Pn
- #
- previousJMinimum = numpy.finfo(float).max
- #
- # Predimensionnement
- Xn_predicted = numpy.asmatrix(numpy.zeros((__n,__m)))
- HX_predicted = numpy.asmatrix(numpy.zeros((__p,__m)))
- #
- for step in range(duration-1):
- if hasattr(Y,"store"):
- Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
- else:
- Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
- #
- if U is not None:
- if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
- elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
- else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
- else:
- Un = None
- #
- if selfA._parameters["InflationType"] == "MultiplicativeOnBackgroundAnomalies":
- Xn = CovarianceInflation( Xn,
- selfA._parameters["InflationType"],
- selfA._parameters["InflationFactor"],
- )
- #
- if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
- EMX = M( [(Xn[:,i], Un) for i in range(__m)], argsAsSerie = True )
- for i in range(__m):
- qi = numpy.asmatrix(numpy.random.multivariate_normal(numpy.zeros(__n), Qn, (1,1,1))).T
- Xn_predicted[:,i] = numpy.asmatrix(numpy.ravel( EMX[i] )).T + qi
- EHX = H( [(Xn_predicted[:,i], Un) for i in range(__m)], argsAsSerie = True )
- for i in range(__m):
- HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( EHX[i] )).T
- if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
- Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
- elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
- # --- > Par principe, M = Id, Q = 0
- Xn_predicted = Xn
- EHX = H( [(Xn_predicted[:,i], Un) for i in range(__m)], argsAsSerie = True )
- for i in range(__m):
- HX_predicted[:,i] = numpy.asmatrix(numpy.ravel( EHX[i] )).T
- #
- # Mean of forecast and observation of forecast
- Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float')
- Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float')
- #
- EaX = (Xn_predicted - Xfm.reshape((__n,-1)))
- EaHX = (HX_predicted - Hfm.reshape((__p,-1)))
- #
- #--------------------------
- if KorV == "KalmanFilterFormula":
- EaX = EaX / numpy.sqrt(__m-1)
- mS = RIdemi * EaHX / numpy.sqrt(__m-1)
- delta = RIdemi * ( Ynpu.reshape((__p,-1)) - Hfm.reshape((__p,-1)) )
- mT = numpy.linalg.inv( numpy.eye(__m) + mS.T @ mS )
- vw = mT @ mS.transpose() @ delta
- #
- Tdemi = numpy.real(scipy.linalg.sqrtm(mT))
- mU = numpy.eye(__m)
- #
- Xn = Xfm.reshape((__n,-1)) + EaX @ ( vw.reshape((__m,-1)) + numpy.sqrt(__m-1) * Tdemi @ mU )
- #--------------------------
- elif KorV == "Variational":
- HXfm = H((Xfm, Un)) # Eventuellement Hfm
- def CostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _Jo = 0.5 * _A.T * RI * _A
- _Jb = 0.5 * (__m-1) * w.T @ w
- _J = _Jo + _Jb
- return float(_J)
- def GradientOfCostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _GardJo = - EaHX.T * RI * _A
- _GradJb = (__m-1) * w.reshape((__m,1))
- _GradJ = _GardJo + _GradJb
- return numpy.ravel(_GradJ)
- vw = scipy.optimize.fmin_cg(
- f = CostFunction,
- x0 = numpy.zeros(__m),
- fprime = GradientOfCostFunction,
- args = (),
- disp = False,
- )
- #
- Hto = EaHX.T * RI * EaHX
- Htb = (__m-1) * numpy.eye(__m)
- Hta = Hto + Htb
- #
- Pta = numpy.linalg.inv( Hta )
- EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
- #
- Xn = Xfm.reshape((__n,-1)) + EaX @ (vw.reshape((__m,-1)) + EWa)
- #--------------------------
- elif KorV == "FiniteSize11": # Jauge Boc2011
- HXfm = H((Xfm, Un)) # Eventuellement Hfm
- def CostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _Jo = 0.5 * _A.T * RI * _A
- _Jb = 0.5 * __m * math.log(1 + 1/__m + w.T @ w)
- _J = _Jo + _Jb
- return float(_J)
- def GradientOfCostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _GardJo = - EaHX.T * RI * _A
- _GradJb = __m * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
- _GradJ = _GardJo + _GradJb
- return numpy.ravel(_GradJ)
- vw = scipy.optimize.fmin_cg(
- f = CostFunction,
- x0 = numpy.zeros(__m),
- fprime = GradientOfCostFunction,
- args = (),
- disp = False,
- )
- #
- Hto = EaHX.T * RI * EaHX
- Htb = __m * \
- ( (1 + 1/__m + vw.T @ vw) * numpy.eye(__m) - 2 * vw @ vw.T ) \
- / (1 + 1/__m + vw.T @ vw)**2
- Hta = Hto + Htb
- #
- Pta = numpy.linalg.inv( Hta )
- EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
- #
- Xn = Xfm.reshape((__n,-1)) + EaX @ (vw.reshape((__m,-1)) + EWa)
- #--------------------------
- elif KorV == "FiniteSize15": # Jauge Boc2015
- HXfm = H((Xfm, Un)) # Eventuellement Hfm
- def CostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _Jo = 0.5 * _A.T * RI * _A
- _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w)
- _J = _Jo + _Jb
- return float(_J)
- def GradientOfCostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _GardJo = - EaHX.T * RI * _A
- _GradJb = (__m+1) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w)
- _GradJ = _GardJo + _GradJb
- return numpy.ravel(_GradJ)
- vw = scipy.optimize.fmin_cg(
- f = CostFunction,
- x0 = numpy.zeros(__m),
- fprime = GradientOfCostFunction,
- args = (),
- disp = False,
- )
- #
- Hto = EaHX.T * RI * EaHX
- Htb = (__m+1) * \
- ( (1 + 1/__m + vw.T @ vw) * numpy.eye(__m) - 2 * vw @ vw.T ) \
- / (1 + 1/__m + vw.T @ vw)**2
- Hta = Hto + Htb
- #
- Pta = numpy.linalg.inv( Hta )
- EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
- #
- Xn = Xfm.reshape((__n,-1)) + EaX @ (vw.reshape((__m,-1)) + EWa)
- #--------------------------
- elif KorV == "FiniteSize16": # Jauge Boc2016
- HXfm = H((Xfm, Un)) # Eventuellement Hfm
- def CostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _Jo = 0.5 * _A.T * RI * _A
- _Jb = 0.5 * (__m+1) * math.log(1 + 1/__m + w.T @ w / (__m-1))
- _J = _Jo + _Jb
- return float(_J)
- def GradientOfCostFunction(w):
- _A = Ynpu.reshape((__p,-1)) - HXfm.reshape((__p,-1)) - (EaHX @ w).reshape((__p,-1))
- _GardJo = - EaHX.T * RI * _A
- _GradJb = ((__m+1) / (__m-1)) * w.reshape((__m,1)) / (1 + 1/__m + w.T @ w / (__m-1))
- _GradJ = _GardJo + _GradJb
- return numpy.ravel(_GradJ)
- vw = scipy.optimize.fmin_cg(
- f = CostFunction,
- x0 = numpy.zeros(__m),
- fprime = GradientOfCostFunction,
- args = (),
- disp = False,
- )
- #
- Hto = EaHX.T * RI * EaHX
- Htb = ((__m+1) / (__m-1)) * \
- ( (1 + 1/__m + vw.T @ vw / (__m-1)) * numpy.eye(__m) - 2 * vw @ vw.T / (__m-1) ) \
- / (1 + 1/__m + vw.T @ vw / (__m-1))**2
- Hta = Hto + Htb
- #
- Pta = numpy.linalg.inv( Hta )
- EWa = numpy.real(scipy.linalg.sqrtm((__m-1)*Pta)) # Partie imaginaire ~= 10^-18
- #
- Xn = Xfm.reshape((__n,-1)) + EaX @ (vw.reshape((__m,-1)) + EWa)
- #--------------------------
- else:
- raise ValueError("KorV has to be chosen in the authorized methods list.")
- #
- if selfA._parameters["InflationType"] == "MultiplicativeOnAnalysisAnomalies":
- Xn = CovarianceInflation( Xn,
- selfA._parameters["InflationType"],
- selfA._parameters["InflationFactor"],
- )
- #
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float')
- #--------------------------
- #
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("APosterioriCovariance") \
- or selfA._toStore("InnovationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
- _Innovation = Ynpu - _HXa
- #
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- # ---> avec analysis
- selfA.StoredVariables["Analysis"].store( Xa )
- if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
- selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
- if selfA._toStore("InnovationAtCurrentAnalysis"):
- selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
- # ---> avec current state
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CurrentState"):
- selfA.StoredVariables["CurrentState"].store( Xn )
- if selfA._toStore("ForecastState"):
- selfA.StoredVariables["ForecastState"].store( Xn_predicted )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
- if selfA._toStore("InnovationAtCurrentState"):
- selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
- if selfA._toStore("SimulatedObservationAtCurrentState") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
- # ---> autres
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
- J = Jb + Jo
- selfA.StoredVariables["CostFunctionJb"].store( Jb )
- selfA.StoredVariables["CostFunctionJo"].store( Jo )
- selfA.StoredVariables["CostFunctionJ" ].store( J )
- #
- if selfA._toStore("IndexOfOptimum") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("CostFunctionJAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
- if selfA._toStore("IndexOfOptimum"):
- selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
- if selfA._toStore("CurrentOptimum"):
- selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
- if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
- if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
- if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
- if selfA._toStore("CostFunctionJAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
- if selfA._toStore("APosterioriCovariance"):
- Eai = (1/numpy.sqrt(__m-1)) * (Xn - Xa.reshape((__n,-1))) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- if selfA._parameters["EstimationOf"] == "Parameters" \
- and J < previousJMinimum:
- previousJMinimum = J
- XaMin = Xa
- if selfA._toStore("APosterioriCovariance"):
- covarianceXaMin = Pn
- #
- # Stockage final supplémentaire de l'optimum en estimation de paramètres
- # ----------------------------------------------------------------------
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- selfA.StoredVariables["Analysis"].store( XaMin )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
- #
- return 0
+ return __Vector
# ==============================================================================
-def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, BnotT=False, _epsilon=1.e-1, _e=1.e-7, _jmax=15000):
- """
- Maximum Likelihood Ensemble Filter (EnKF/MLEF Zupanski 2005, Bocquet 2013)
+def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __Betaf):
+ "Recentre l'ensemble Xn autour de l'analyse 3DVAR"
+ #
+ Xf = EnsembleMean( __EnXf )
+ Pf = Covariance( asCovariance=EnsembleErrorCovariance(__EnXf) )
+ Pf = (1 - __Betaf) * __B.asfullmatrix(Xf.size) + __Betaf * Pf
+ #
+ selfB = PartialAlgorithm("3DVAR")
+ selfB._parameters["Minimizer"] = "LBFGSB"
+ selfB._parameters["MaximumNumberOfSteps"] = 15000
+ selfB._parameters["CostDecrementTolerance"] = 1.e-7
+ selfB._parameters["ProjectedGradientTolerance"] = -1
+ selfB._parameters["GradientNormTolerance"] = 1.e-05
+ selfB._parameters["StoreInternalVariables"] = False
+ selfB._parameters["optiprint"] = -1
+ selfB._parameters["optdisp"] = 0
+ selfB._parameters["Bounds"] = None
+ selfB._parameters["InitializationPoint"] = Xf
+ from daAlgorithms.Atoms import std3dvar
+ std3dvar.std3dvar(selfB, Xf, __Ynpu, None, __HO, None, __R, Pf)
+ Xa = selfB.get("Analysis")[-1].reshape((-1,1))
+ del selfB
+ #
+ return Xa + EnsembleOfAnomalies( __EnXn )
- selfA est identique au "self" d'algorithme appelant et contient les
- valeurs.
+# ==============================================================================
+def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle,
+ __CovForecast = False, __LinEvolution = False,
+ ):
+ """
+ Prévision multi-pas avec une correction par pas (multi-méthodes)
"""
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA._parameters["StoreInternalVariables"] = True
- #
- # Opérateurs
- # ----------
- H = HO["Direct"].appliedControledFormTo
#
+ # Initialisation
+ # --------------
if selfA._parameters["EstimationOf"] == "State":
- M = EM["Direct"].appliedControledFormTo
- #
- if CM is not None and "Tangent" in CM and U is not None:
- Cm = CM["Tangent"].asMatrix(Xb)
+ if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
+ Xn = numpy.asarray(Xb)
+ if __CovForecast: Pn = B
+ selfA.StoredVariables["Analysis"].store( Xn )
+ if selfA._toStore("APosterioriCovariance"):
+ if hasattr(B,"asfullmatrix"):
+ selfA.StoredVariables["APosterioriCovariance"].store( B.asfullmatrix(Xn.size) )
+ else:
+ selfA.StoredVariables["APosterioriCovariance"].store( B )
+ selfA._setInternalState("seed", numpy.random.get_state())
+ elif selfA._parameters["nextStep"]:
+ Xn = selfA._getInternalState("Xn")
+ if __CovForecast: Pn = selfA._getInternalState("Pn")
else:
- Cm = None
+ Xn = numpy.asarray(Xb)
+ if __CovForecast: Pn = B
#
- # Nombre de pas identique au nombre de pas d'observations
- # -------------------------------------------------------
if hasattr(Y,"stepnumber"):
duration = Y.stepnumber()
- __p = numpy.cumprod(Y.shape())[-1]
else:
duration = 2
- __p = numpy.array(Y).size
- #
- # Précalcul des inversions de B et R
- # ----------------------------------
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- BI = B.getI()
- RI = R.getI()
- #
- # Initialisation
- # --------------
- __n = Xb.size
- __m = selfA._parameters["NumberOfMembers"]
- Xn = numpy.asmatrix(numpy.dot( Xb.reshape(__n,1), numpy.ones((1,__m)) ))
- if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n)
- else: Pn = B
- if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p)
- else: Rn = R
- if hasattr(Q,"asfullmatrix"): Qn = Q.asfullmatrix(__n)
- else: Qn = Q
- #
- if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]:
- selfA.StoredVariables["Analysis"].store( numpy.ravel(Xb) )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- covarianceXa = Pn
- #
- previousJMinimum = numpy.finfo(float).max
- #
- # Predimensionnement
- Xn_predicted = numpy.asmatrix(numpy.zeros((__n,__m)))
#
+ # Multi-steps
+ # -----------
for step in range(duration-1):
+ selfA.StoredVariables["CurrentStepNumber"].store( len(selfA.StoredVariables["Analysis"]) )
+ #
if hasattr(Y,"store"):
- Ynpu = numpy.asmatrix(numpy.ravel( Y[step+1] )).T
+ Ynpu = numpy.asarray( Y[step+1] ).reshape((-1,1))
else:
- Ynpu = numpy.asmatrix(numpy.ravel( Y )).T
+ Ynpu = numpy.asarray( Y ).reshape((-1,1))
#
if U is not None:
if hasattr(U,"store") and len(U)>1:
- Un = numpy.asmatrix(numpy.ravel( U[step] )).T
+ Un = numpy.asarray( U[step] ).reshape((-1,1))
elif hasattr(U,"store") and len(U)==1:
- Un = numpy.asmatrix(numpy.ravel( U[0] )).T
+ Un = numpy.asarray( U[0] ).reshape((-1,1))
else:
- Un = numpy.asmatrix(numpy.ravel( U )).T
+ Un = numpy.asarray( U ).reshape((-1,1))
else:
Un = None
#
- if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast
- EMX = M( [(Xn[:,i], Un) for i in range(__m)], argsAsSerie = True )
- for i in range(__m):
- qi = numpy.random.multivariate_normal(numpy.zeros(__n), Qn)
- Xn_predicted[:,i] = (numpy.ravel( EMX[i] ) + qi).reshape((__n,-1))
- if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon !
- Cm = Cm.reshape(__n,Un.size) # ADAO & check shape
- Xn_predicted = Xn_predicted + Cm * Un
- elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast
+ # Predict (Time Update)
+ # ---------------------
+ if selfA._parameters["EstimationOf"] == "State":
+ if __CovForecast or __LinEvolution:
+ Mt = EM["Tangent"].asMatrix(Xn)
+ Mt = Mt.reshape(Xn.size,Xn.size) # ADAO & check shape
+ if __CovForecast:
+ Ma = EM["Adjoint"].asMatrix(Xn)
+ Ma = Ma.reshape(Xn.size,Xn.size) # ADAO & check shape
+ Pn_predicted = Q + Mt @ (Pn @ Ma)
+ if __LinEvolution:
+ Xn_predicted = Mt @ Xn
+ else:
+ M = EM["Direct"].appliedControledFormTo
+ Xn_predicted = M( (Xn, Un) )
+ if CM is not None and "Tangent" in CM and Un is not None: # Attention : si Cm est aussi dans M, doublon !
+ Cm = CM["Tangent"].asMatrix(Xn_predicted)
+ Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape
+ Xn_predicted = Xn_predicted + (Cm @ Un).reshape((-1,1))
+ elif selfA._parameters["EstimationOf"] == "Parameters": # No forecast
# --- > Par principe, M = Id, Q = 0
Xn_predicted = Xn
- #
- # Mean of forecast and observation of forecast
- Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float')
- #
- EaX = (Xn_predicted - Xfm.reshape((__n,-1))) / numpy.sqrt(__m-1)
- #
- #--------------------------
- Ua = numpy.eye(__m)
- Ta = numpy.eye(__m)
- #
- __j = 0 # 4:
- vw = numpy.zeros(__m) # 4:
- Deltaw = 1
- while numpy.linalg.norm(Deltaw) >= _e or __j >= _jmax: # 5: et 19:
- vx = numpy.ravel(Xfm) + EaX @ vw # 6:
- #
- if BnotT:
- EE = vx.reshape((__n,-1)) + _epsilon * EaX # 7:
- else:
- EE = vx.reshape((__n,-1)) + numpy.sqrt(__m-1) * EaX @ Ta # 8:
- #
- EZ = H( [(EE[:,i], Un) for i in range(__m)],
- argsAsSerie = True,
- returnSerieAsArrayMatrix = True )
- #
- ybar = EZ.mean(axis=1, dtype=mfp).astype('float').reshape((__p,-1)) # 10: Observation mean
- #
- if BnotT:
- EY = (EZ - ybar) / _epsilon # 11:
+ if __CovForecast: Pn_predicted = Pn
+ Xn_predicted = numpy.asarray(Xn_predicted).reshape((-1,1))
+ if selfA._toStore("ForecastState"):
+ selfA.StoredVariables["ForecastState"].store( Xn_predicted )
+ if __CovForecast:
+ if hasattr(Pn_predicted,"asfullmatrix"):
+ Pn_predicted = Pn_predicted.asfullmatrix(Xn.size)
else:
- EY = ( (EZ - ybar) @ numpy.linalg.inv(Ta) ) / numpy.sqrt(__m-1) # 12:
- #
- GradJ = numpy.ravel(vw.reshape((__m,1)) - EY.transpose() @ (RI * (Ynpu - ybar))) # 13:
- mH = numpy.eye(__m) + EY.transpose() @ (RI * EY) # 14:
- Deltaw = numpy.linalg.solve(mH,GradJ) # 15:
- vw = vw - Deltaw # 16:
- if not BnotT:
- Ta = numpy.linalg.inv(numpy.real(scipy.linalg.sqrtm( mH ))) # 17:
- __j = __j + 1 # 18:
- #
- if BnotT:
- Ta = numpy.linalg.inv(numpy.real(scipy.linalg.sqrtm( mH ))) # 20:
- #
- Xn = vx.reshape((__n,-1)) + numpy.sqrt(__m-1) * EaX @ Ta @ Ua # 21:
+ Pn_predicted = numpy.asarray(Pn_predicted).reshape((Xn.size,Xn.size))
+ if selfA._toStore("ForecastCovariance"):
+ selfA.StoredVariables["ForecastCovariance"].store( Pn_predicted )
+ #
+ # Correct (Measurement Update)
+ # ----------------------------
+ if __CovForecast:
+ oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, Pn_predicted, True)
+ else:
+ oneCycle(selfA, Xn_predicted, Ynpu, Un, HO, CM, R, B, True)
#
- Xa = Xn.mean(axis=1, dtype=mfp).astype('float')
#--------------------------
- #
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("APosterioriCovariance") \
- or selfA._toStore("InnovationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T
- _Innovation = Ynpu - _HXa
- #
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- # ---> avec analysis
- selfA.StoredVariables["Analysis"].store( Xa )
- if selfA._toStore("SimulatedObservationAtCurrentAnalysis"):
- selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"].store( _HXa )
- if selfA._toStore("InnovationAtCurrentAnalysis"):
- selfA.StoredVariables["InnovationAtCurrentAnalysis"].store( _Innovation )
- # ---> avec current state
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CurrentState"):
- selfA.StoredVariables["CurrentState"].store( Xn )
- if selfA._toStore("ForecastState"):
- selfA.StoredVariables["ForecastState"].store( Xn_predicted )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( Xn_predicted - Xa )
- #~ if selfA._toStore("InnovationAtCurrentState"):
- #~ selfA.StoredVariables["InnovationAtCurrentState"].store( - HX_predicted + Ynpu )
- #~ if selfA._toStore("SimulatedObservationAtCurrentState") \
- #~ or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- #~ selfA.StoredVariables["SimulatedObservationAtCurrentState"].store( HX_predicted )
- # ---> autres
- if selfA._parameters["StoreInternalVariables"] \
- or selfA._toStore("CostFunctionJ") \
- or selfA._toStore("CostFunctionJb") \
- or selfA._toStore("CostFunctionJo") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("APosterioriCovariance"):
- Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) )
- Jo = float( 0.5 * _Innovation.T * RI * _Innovation )
- J = Jb + Jo
- selfA.StoredVariables["CostFunctionJb"].store( Jb )
- selfA.StoredVariables["CostFunctionJo"].store( Jo )
- selfA.StoredVariables["CostFunctionJ" ].store( J )
- #
- if selfA._toStore("IndexOfOptimum") \
- or selfA._toStore("CurrentOptimum") \
- or selfA._toStore("CostFunctionJAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJbAtCurrentOptimum") \
- or selfA._toStore("CostFunctionJoAtCurrentOptimum") \
- or selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- IndexMin = numpy.argmin( selfA.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
- if selfA._toStore("IndexOfOptimum"):
- selfA.StoredVariables["IndexOfOptimum"].store( IndexMin )
- if selfA._toStore("CurrentOptimum"):
- selfA.StoredVariables["CurrentOptimum"].store( selfA.StoredVariables["Analysis"][IndexMin] )
- if selfA._toStore("SimulatedObservationAtCurrentOptimum"):
- selfA.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( selfA.StoredVariables["SimulatedObservationAtCurrentAnalysis"][IndexMin] )
- if selfA._toStore("CostFunctionJbAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJb"][IndexMin] )
- if selfA._toStore("CostFunctionJoAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( selfA.StoredVariables["CostFunctionJo"][IndexMin] )
- if selfA._toStore("CostFunctionJAtCurrentOptimum"):
- selfA.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( selfA.StoredVariables["CostFunctionJ" ][IndexMin] )
- if selfA._toStore("APosterioriCovariance"):
- Eai = (1/numpy.sqrt(__m-1)) * (Xn - Xa.reshape((__n,-1))) # Anomalies
- Pn = Eai @ Eai.T
- Pn = 0.5 * (Pn + Pn.T)
- selfA.StoredVariables["APosterioriCovariance"].store( Pn )
- if selfA._parameters["EstimationOf"] == "Parameters" \
- and J < previousJMinimum:
- previousJMinimum = J
- XaMin = Xa
- if selfA._toStore("APosterioriCovariance"):
- covarianceXaMin = Pn
- #
- # Stockage final supplémentaire de l'optimum en estimation de paramètres
- # ----------------------------------------------------------------------
- if selfA._parameters["EstimationOf"] == "Parameters":
- selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) )
- selfA.StoredVariables["Analysis"].store( XaMin )
- if selfA._toStore("APosterioriCovariance"):
- selfA.StoredVariables["APosterioriCovariance"].store( covarianceXaMin )
- if selfA._toStore("BMA"):
- selfA.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(XaMin) )
+ Xn = selfA._getInternalState("Xn")
+ if __CovForecast: Pn = selfA._getInternalState("Pn")
#
return 0
