X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FdaComposant%2FdaCore%2FNumericObjects.py;h=aa6c683592e0c1a146b1ed356675970ae5d8defd;hb=ba52f8705ad24049d6fa5957da1cbe00d7bcb4d6;hp=85ecff2048846829367362bc70444f4bf66b1b06;hpb=018c6093cf0f1f41fb2a99f3c32cfa4c7cdb9d12;p=modules%2Fadao.git diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 85ecff2..aa6c683 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -1,6 +1,6 @@ # -*- 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 @@ -27,7 +27,7 @@ __author__ = "Jean-Philippe ARGAUD" import os, time, copy, types, sys, logging import math, numpy, scipy, scipy.optimize, scipy.version -from daCore.BasicObjects import Operator +from daCore.BasicObjects import Operator, Covariance, PartialAlgorithm from daCore.PlatformInfo import PlatformInfo mpr = PlatformInfo().MachinePrecision() mfp = PlatformInfo().MaximumPrecision() @@ -66,6 +66,7 @@ class FDApproximation(object): increment = 0.01, dX = None, extraArguments = None, + reducingMemoryUse = False, avoidingRedundancy = True, toleranceInRedundancy = 1.e-18, lenghtOfRedundancy = -1, @@ -75,6 +76,7 @@ class FDApproximation(object): ): self.__name = str(name) self.__extraArgs = extraArguments + # if mpEnabled: try: import multiprocessing @@ -88,12 +90,12 @@ class FDApproximation(object): 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) @@ -105,6 +107,9 @@ class FDApproximation(object): 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): @@ -150,9 +155,6 @@ class FDApproximation(object): self.__dX = None else: self.__dX = numpy.ravel( dX ) - 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) # --------------------------------------------------------- def __doublon__(self, e, l, n, v=None): @@ -164,6 +166,29 @@ class FDApproximation(object): break return __ac, __iac + # --------------------------------------------------------- + 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 ): """ @@ -181,7 +206,7 @@ class FDApproximation(object): 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 @@ -240,6 +265,11 @@ class FDApproximation(object): 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: @@ -359,24 +389,29 @@ class FDApproximation(object): _HX_plus_dXi = self.DirectOperator( _X_plus_dXi ) # _Jacobienne.append( numpy.ravel(( _HX_plus_dXi - _HX ) / _dXi) ) - # # - _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 (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 @@ -389,26 +424,25 @@ class FDApproximation(object): 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.ravel( dX ) - _HtX = numpy.dot(_Jacobienne, _dX) + _HtX = self.TangentMatrix( X, dotWith = dX ) if self.__mfEnabled: return [_HtX,] else: return _HtX @@ -421,26 +455,25 @@ class FDApproximation(object): 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.ravel( Y ) - _HaY = numpy.dot(_JacobienneT, _Y) + _HaY = self.TangentMatrix( X, dotTWith = Y ) if self.__mfEnabled: return [_HaY,] else: return _HaY @@ -453,12 +486,12 @@ def EnsembleOfCenteredPerturbations( _bgcenter, _bgcovariance, _nbmembers ): raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),)) # if _bgcovariance is None: - BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Perturbations = numpy.tile( _bgcenter, _nbmembers) else: _Z = numpy.random.multivariate_normal(numpy.zeros(_bgcenter.size), _bgcovariance, size=_nbmembers).T - BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Z + _Perturbations = numpy.tile( _bgcenter, _nbmembers) + _Z # - return BackgroundEnsemble + return _Perturbations # ============================================================================== def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _withSVD = True): @@ -480,29 +513,29 @@ def EnsembleOfBackgroundPerturbations( _bgcenter, _bgcovariance, _nbmembers, _wi if _nbmembers < 1: raise ValueError("Number of members has to be strictly more than 1 (given number: %s)."%(str(_nbmembers),)) if _bgcovariance is None: - BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Perturbations = numpy.tile( _bgcenter, _nbmembers) else: if _withSVD: - U, s, V = numpy.linalg.svd(_bgcovariance, full_matrices=False) + _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.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]) + _Z = numpy.dot(numpy.diag(numpy.sqrt(_s[:_nbmembers-1])), _V[:_nbmembers-1]) _Zca = __CenteredRandomAnomalies(_Z, _nbmembers) - BackgroundEnsemble = _bgcenter + _Zca + _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) - BackgroundEnsemble = _bgcenter + _Zca + _Perturbations = _bgcenter + _Zca else: - BackgroundEnsemble = numpy.tile( _bgcenter, _nbmembers) + _Perturbations = numpy.tile( _bgcenter, _nbmembers) # - return BackgroundEnsemble + return _Perturbations # ============================================================================== def EnsembleMean( __Ensemble ): @@ -708,7 +741,7 @@ def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None): 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 = Hm( Xr ) + Yr = numpy.asarray(Hm( Xr )) else: raise ValueError("Quantile simulations has only to be Linear or NonLinear.") # @@ -753,7 +786,7 @@ def RecentredBounds( __Bounds, __Center): # Conserve une valeur par défaut à None s'il n'y a pas de bornes if __Bounds is None: return None # Recentre les valeurs numériques de bornes - return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).transpose((-1,1)) + return ForceNumericBounds( __Bounds ) - numpy.ravel( __Center ).reshape((-1,1)) # ============================================================================== def ApplyBounds( __Vector, __Bounds, __newClip = True): @@ -766,7 +799,7 @@ def ApplyBounds( __Vector, __Bounds, __newClip = True): 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 to be applied for this vector") + 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") # @@ -782,6 +815,31 @@ def ApplyBounds( __Vector, __Bounds, __newClip = True): # return __Vector +# ============================================================================== +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 + __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 + std3dvar(selfB, Xf, __Ynpu, None, __HO, None, None, __R, Pf, None) + Xa = selfB.get("Analysis")[-1].reshape((-1,1)) + del selfB + # + return Xa + EnsembleOfAnomalies( __EnXn ) + # ============================================================================== def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): """ @@ -844,6 +902,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): RI = R.getI() # __n = Xb.size + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = Xb @@ -893,7 +952,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1)) if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon ! Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape - XEtnnpi = XEtnnpi + Cm * Un + XEtnnpi = XEtnnpi + Cm @ Un if selfA._parameters["Bounds"] is not None and selfA._parameters["ConstrainedBy"] == "EstimateProjection": XEtnnpi = ApplyBounds( XEtnnpi, selfA._parameters["Bounds"] ) elif selfA._parameters["EstimationOf"] == "Parameters": @@ -943,7 +1002,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): _Innovation = Ynpu - Yncm.reshape((-1,1)) if selfA._parameters["EstimationOf"] == "Parameters": if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! - _Innovation = _Innovation - Cm * Un + _Innovation = _Innovation - Cm @ Un # Kn = Pxyn * Pyyn.I Xn = Xncm.reshape((-1,1)) + Kn * _Innovation @@ -987,7 +1046,7 @@ def c2ukf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): or selfA._toStore("CostFunctionJo") \ or selfA._toStore("CurrentOptimum") \ or selfA._toStore("APosterioriCovariance"): - Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) ) + 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 ) @@ -1073,6 +1132,7 @@ def cekf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): RI = R.getI() # __n = Xb.size + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = Xb @@ -1127,7 +1187,7 @@ def cekf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).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 + Xn_predicted = Xn_predicted + Cm @ Un Pn_predicted = Q + Mt * (Pn * Ma) elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 @@ -1144,7 +1204,7 @@ def cekf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1)) _Innovation = Ynpu - HX_predicted if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! - _Innovation = _Innovation - Cm * Un + _Innovation = _Innovation - Cm @ Un # Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha)) Xn = Xn_predicted + Kn * _Innovation @@ -1320,7 +1380,7 @@ def enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterForm returnSerieAsArrayMatrix = True ) 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 - EZ = EZ + Cm * Un + EZ = EZ + Cm @ Un elif selfA._parameters["EstimationOf"] == "Parameters": # --- > Par principe, M = Id, Q = 0 EZ = H( [(EL[:,i], Un) for i in range(__m)], @@ -1376,7 +1436,10 @@ def enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterForm return 0 # ============================================================================== -def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): +def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, + VariantM="KalmanFilterFormula", + Hybrid=None, + ): """ Ensemble-Transform EnKF """ @@ -1418,6 +1481,8 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): # __n = Xb.size __m = selfA._parameters["NumberOfMembers"] + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) + previousJMinimum = numpy.finfo(float).max # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m ) @@ -1431,8 +1496,6 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): elif selfA._parameters["nextStep"]: Xn = selfA._getInternalState("Xn") # - previousJMinimum = numpy.finfo(float).max - # for step in range(duration-1): numpy.random.set_state(selfA._getInternalState("seed")) if hasattr(Y,"store"): @@ -1466,7 +1529,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): returnSerieAsArrayMatrix = True ) 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 + Xn_predicted = Xn_predicted + Cm @ Un elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 Xn_predicted = EMX = Xn @@ -1475,8 +1538,8 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): returnSerieAsArrayMatrix = True ) # # Mean of forecast and observation of forecast - Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) - Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1)) + Xfm = EnsembleMean( Xn_predicted ) + Hfm = EnsembleMean( HX_predicted ) # # Anomalies EaX = EnsembleOfAnomalies( Xn_predicted, Xfm ) @@ -1564,7 +1627,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): HXfm = H((Xfm[:,None], Un)) # Eventuellement Hfm def CostFunction(w): _A = Ynpu - HXfm.reshape((__p,1)) - (EaHX @ w).reshape((__p,1)) - _Jo = 0.5 * _A.T * RI * _A + _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) @@ -1635,7 +1698,11 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): selfA._parameters["InflationFactor"], ) # - Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) + if Hybrid == "E3DVAR": + betaf = selfA._parameters["HybridCovarianceEquilibrium"] + Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf) + # + Xa = EnsembleMean( Xn ) #-------------------------- selfA._setInternalState("Xn", Xn) selfA._setInternalState("seed", numpy.random.get_state()) @@ -1649,7 +1716,7 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): or selfA._toStore("InnovationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentOptimum"): - _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T + _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1)) _Innovation = Ynpu - _HXa # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) ) @@ -1681,8 +1748,8 @@ def etkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula"): 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 ) + 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 ) @@ -1769,6 +1836,7 @@ def exkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): RI = R.getI() # __n = Xb.size + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = Xb @@ -1820,7 +1888,7 @@ def exkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): Xn_predicted = numpy.ravel( M( (Xn, Un) ) ).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 + Xn_predicted = Xn_predicted + Cm @ Un Pn_predicted = Q + Mt * (Pn * Ma) elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 @@ -1834,7 +1902,7 @@ def exkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): HX_predicted = numpy.ravel( H( (Xn_predicted, Un) ) ).reshape((__p,1)) _Innovation = Ynpu - HX_predicted if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! - _Innovation = _Innovation - Cm * Un + _Innovation = _Innovation - Cm @ Un # Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha)) Xn = Xn_predicted + Kn * _Innovation @@ -1962,6 +2030,8 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", # __n = Xb.size __m = selfA._parameters["NumberOfMembers"] + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) + previousJMinimum = numpy.finfo(float).max # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: if hasattr(B,"asfullmatrix"): Pn = B.asfullmatrix(__n) @@ -1977,8 +2047,6 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", elif selfA._parameters["nextStep"]: Xn = selfA._getInternalState("Xn") # - previousJMinimum = numpy.finfo(float).max - # for step in range(duration-1): numpy.random.set_state(selfA._getInternalState("seed")) if hasattr(Y,"store"): @@ -1988,11 +2056,11 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", # if U is not None: if hasattr(U,"store") and len(U)>1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # @@ -2067,7 +2135,7 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", selfA._parameters["InflationFactor"], ) # - Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) + Xa = EnsembleMean( Xn ) #-------------------------- selfA._setInternalState("Xn", Xn) selfA._setInternalState("seed", numpy.random.get_state()) @@ -2081,7 +2149,7 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", or selfA._toStore("InnovationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentOptimum"): - _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T + _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1)) _Innovation = Ynpu - _HXa # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) ) @@ -2113,8 +2181,8 @@ def ienkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="IEnKF12", 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 ) + 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 ) @@ -2171,18 +2239,12 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): # # Initialisations # --------------- - # - # Opérateur non-linéaire pour la boucle externe Hm = HO["Direct"].appliedTo # - # Précalcul des inversions de B et R BI = B.getI() RI = R.getI() # - # Point de démarrage de l'optimisation - Xini = selfA._parameters["InitializationPoint"] - # - HXb = numpy.asmatrix(numpy.ravel( Hm( Xb ) )).T + HXb = numpy.asarray(Hm( Xb )).reshape((-1,1)) Innovation = Y - HXb # # Outer Loop @@ -2190,7 +2252,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): iOuter = 0 J = 1./mpr DeltaJ = 1./mpr - Xr = Xini.reshape((-1,1)) + Xr = numpy.asarray(selfA._parameters["InitializationPoint"]).reshape((-1,1)) while abs(DeltaJ) >= selfA._parameters["CostDecrementTolerance"] and iOuter <= selfA._parameters["MaximumNumberOfSteps"]: # # Inner Loop @@ -2201,13 +2263,12 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): # Définition de la fonction-coût # ------------------------------ def CostFunction(dx): - _dX = numpy.asmatrix(numpy.ravel( dx )).T + _dX = numpy.asarray(dx).reshape((-1,1)) if selfA._parameters["StoreInternalVariables"] or \ selfA._toStore("CurrentState") or \ selfA._toStore("CurrentOptimum"): selfA.StoredVariables["CurrentState"].store( Xb + _dX ) - _HdX = Ht * _dX - _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T + _HdX = (Ht @ _dX).reshape((-1,1)) _dInnovation = Innovation - _HdX if selfA._toStore("SimulatedObservationAtCurrentState") or \ selfA._toStore("SimulatedObservationAtCurrentOptimum"): @@ -2215,8 +2276,8 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): if selfA._toStore("InnovationAtCurrentState"): selfA.StoredVariables["InnovationAtCurrentState"].store( _dInnovation ) # - Jb = float( 0.5 * _dX.T * BI * _dX ) - Jo = float( 0.5 * _dInnovation.T * RI * _dInnovation ) + Jb = float( 0.5 * _dX.T * (BI * _dX) ) + Jo = float( 0.5 * _dInnovation.T * (RI * _dInnovation) ) J = Jb + Jo # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) ) @@ -2245,11 +2306,10 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return J # def GradientOfCostFunction(dx): - _dX = numpy.asmatrix(numpy.ravel( dx )).T - _HdX = Ht * _dX - _HdX = numpy.asmatrix(numpy.ravel( _HdX )).T + _dX = numpy.ravel( dx ) + _HdX = (Ht @ _dX).reshape((-1,1)) _dInnovation = Innovation - _HdX - GradJb = BI * _dX + GradJb = BI @ _dX GradJo = - Ht.T @ (RI * _dInnovation) GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo ) return GradJ @@ -2266,7 +2326,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( func = CostFunction, - x0 = numpy.zeros(Xini.size), + x0 = numpy.zeros(Xb.size), fprime = GradientOfCostFunction, args = (), bounds = RecentredBounds(selfA._parameters["Bounds"], Xb), @@ -2280,7 +2340,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): elif selfA._parameters["Minimizer"] == "TNC": Minimum, nfeval, rc = scipy.optimize.fmin_tnc( func = CostFunction, - x0 = numpy.zeros(Xini.size), + x0 = numpy.zeros(Xb.size), fprime = GradientOfCostFunction, args = (), bounds = RecentredBounds(selfA._parameters["Bounds"], Xb), @@ -2292,7 +2352,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): elif selfA._parameters["Minimizer"] == "CG": Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg( f = CostFunction, - x0 = numpy.zeros(Xini.size), + x0 = numpy.zeros(Xb.size), fprime = GradientOfCostFunction, args = (), maxiter = selfA._parameters["MaximumNumberOfSteps"], @@ -2303,7 +2363,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): elif selfA._parameters["Minimizer"] == "NCG": Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg( f = CostFunction, - x0 = numpy.zeros(Xini.size), + x0 = numpy.zeros(Xb.size), fprime = GradientOfCostFunction, args = (), maxiter = selfA._parameters["MaximumNumberOfSteps"], @@ -2314,7 +2374,7 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): elif selfA._parameters["Minimizer"] == "BFGS": Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs( f = CostFunction, - x0 = numpy.zeros(Xini.size), + x0 = numpy.zeros(Xb.size), fprime = GradientOfCostFunction, args = (), maxiter = selfA._parameters["MaximumNumberOfSteps"], @@ -2406,8 +2466,10 @@ def incr3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return 0 # ============================================================================== -def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", - BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000): +def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, + VariantM="MLEF13", BnotT=False, _epsilon=1.e-3, _e=1.e-7, _jmax=15000, + Hybrid=None, + ): """ Maximum Likelihood Ensemble Filter """ @@ -2445,6 +2507,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", # __n = Xb.size __m = selfA._parameters["NumberOfMembers"] + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) + previousJMinimum = numpy.finfo(float).max # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = EnsembleOfBackgroundPerturbations( Xb, None, __m ) @@ -2458,8 +2522,6 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", elif selfA._parameters["nextStep"]: Xn = selfA._getInternalState("Xn") # - previousJMinimum = numpy.finfo(float).max - # for step in range(duration-1): numpy.random.set_state(selfA._getInternalState("seed")) if hasattr(Y,"store"): @@ -2469,11 +2531,11 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", # if U is not None: if hasattr(U,"store") and len(U)>1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # @@ -2490,7 +2552,7 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", Xn_predicted = EnsemblePerturbationWithGivenCovariance( EMX, Q ) 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 + Xn_predicted = Xn_predicted + Cm @ Un elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 Xn_predicted = EMX = Xn @@ -2548,7 +2610,11 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", selfA._parameters["InflationFactor"], ) # - Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) + if Hybrid == "E3DVAR": + betaf = selfA._parameters["HybridCovarianceEquilibrium"] + Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf) + # + Xa = EnsembleMean( Xn ) #-------------------------- selfA._setInternalState("Xn", Xn) selfA._setInternalState("seed", numpy.random.get_state()) @@ -2562,7 +2628,7 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", or selfA._toStore("InnovationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentOptimum"): - _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T + _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1)) _Innovation = Ynpu - _HXa # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) ) @@ -2594,8 +2660,8 @@ def mlef(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="MLEF13", 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 ) + 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 ) @@ -2690,7 +2756,7 @@ def mmqr( 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 + Derivees = Derivees.reshape(n,p) # ADAO & check shape DeriveesT = Derivees.transpose() M = numpy.dot( DeriveesT , (numpy.array(numpy.matrix(p*[poids,]).T)*Derivees) ) SM = numpy.transpose(numpy.dot( DeriveesT , veps )) @@ -2729,8 +2795,13 @@ def multi3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle): """ # # Initialisation + # -------------- if selfA._parameters["EstimationOf"] == "State": - M = EM["Direct"].appliedTo + 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 # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = numpy.ravel(Xb).reshape((-1,1)) @@ -2759,16 +2830,29 @@ def multi3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle): else: Ynpu = numpy.ravel( Y ).reshape((-1,1)) # + if U is not None: + if hasattr(U,"store") and len(U)>1: + Un = numpy.ravel( U[step] ).reshape((-1,1)) + elif hasattr(U,"store") and len(U)==1: + Un = numpy.ravel( U[0] ).reshape((-1,1)) + else: + Un = numpy.ravel( U ).reshape((-1,1)) + else: + Un = None + # if selfA._parameters["EstimationOf"] == "State": # Forecast - Xn_predicted = M( Xn ) + Xn_predicted = M( (Xn, Un) ) if selfA._toStore("ForecastState"): selfA.StoredVariables["ForecastState"].store( Xn_predicted ) + 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": # No forecast # --- > Par principe, M = Id, Q = 0 Xn_predicted = Xn Xn_predicted = numpy.ravel(Xn_predicted).reshape((-1,1)) # - oneCycle(selfA, Xn_predicted, Ynpu, U, HO, None, None, R, B, None) + oneCycle(selfA, Xn_predicted, Ynpu, None, HO, None, None, R, B, None) # Xn = selfA.StoredVariables["Analysis"][-1] #-------------------------- @@ -2784,16 +2868,13 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): # # Initialisations # --------------- - # - # Opérateurs Hm = HO["Direct"].appliedTo # - # Utilisation éventuelle d'un vecteur H(Xb) précalculé if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]: - HXb = Hm( Xb, HO["AppliedInX"]["HXb"] ) + HXb = numpy.asarray(Hm( Xb, HO["AppliedInX"]["HXb"] )) else: - HXb = Hm( Xb ) - HXb = numpy.asmatrix(numpy.ravel( HXb )).T + HXb = numpy.asarray(Hm( Xb )) + HXb = numpy.ravel( HXb ).reshape((-1,1)) if Y.size != HXb.size: raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size)) if max(Y.shape) != max(HXb.shape): @@ -2809,13 +2890,12 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): HBHTpR = R + Ht * BHT Innovation = Y - HXb # - # Point de démarrage de l'optimisation - Xini = numpy.zeros(Xb.shape) + Xini = numpy.zeros(Y.size) # # Définition de la fonction-coût # ------------------------------ def CostFunction(w): - _W = w.reshape((-1,1)) + _W = numpy.asarray(w).reshape((-1,1)) if selfA._parameters["StoreInternalVariables"] or \ selfA._toStore("CurrentState") or \ selfA._toStore("CurrentOptimum"): @@ -2856,7 +2936,7 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return J # def GradientOfCostFunction(w): - _W = w.reshape((-1,1)) + _W = numpy.asarray(w).reshape((-1,1)) GradJb = HBHTpR @ _W GradJo = - Innovation GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo ) @@ -2876,7 +2956,6 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): x0 = Xini, fprime = GradientOfCostFunction, args = (), - bounds = RecentredBounds(selfA._parameters["Bounds"], Xb), maxfun = selfA._parameters["MaximumNumberOfSteps"]-1, factr = selfA._parameters["CostDecrementTolerance"]*1.e14, pgtol = selfA._parameters["ProjectedGradientTolerance"], @@ -2890,7 +2969,6 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): x0 = Xini, fprime = GradientOfCostFunction, args = (), - bounds = RecentredBounds(selfA._parameters["Bounds"], Xb), maxfun = selfA._parameters["MaximumNumberOfSteps"], pgtol = selfA._parameters["ProjectedGradientTolerance"], ftol = selfA._parameters["CostDecrementTolerance"], @@ -3013,7 +3091,10 @@ def psas3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return 0 # ============================================================================== -def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16"): +def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, + VariantM="KalmanFilterFormula16", + Hybrid=None, + ): """ Stochastic EnKF """ @@ -3051,6 +3132,8 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" # __n = Xb.size __m = selfA._parameters["NumberOfMembers"] + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) + previousJMinimum = numpy.finfo(float).max # if hasattr(R,"asfullmatrix"): Rn = R.asfullmatrix(__p) else: Rn = R @@ -3066,8 +3149,6 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" elif selfA._parameters["nextStep"]: Xn = selfA._getInternalState("Xn") # - previousJMinimum = numpy.finfo(float).max - # for step in range(duration-1): numpy.random.set_state(selfA._getInternalState("seed")) if hasattr(Y,"store"): @@ -3077,11 +3158,11 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" # if U is not None: if hasattr(U,"store") and len(U)>1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # @@ -3101,7 +3182,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" returnSerieAsArrayMatrix = True ) 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 + Xn_predicted = Xn_predicted + Cm @ Un elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 Xn_predicted = EMX = Xn @@ -3110,8 +3191,8 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" returnSerieAsArrayMatrix = True ) # # Mean of forecast and observation of forecast - Xfm = Xn_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) - Hfm = HX_predicted.mean(axis=1, dtype=mfp).astype('float').reshape((__p,1)) + Xfm = EnsembleMean( Xn_predicted ) + Hfm = EnsembleMean( HX_predicted ) # #-------------------------- if VariantM == "KalmanFilterFormula05": @@ -3151,7 +3232,11 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" selfA._parameters["InflationFactor"], ) # - Xa = Xn.mean(axis=1, dtype=mfp).astype('float').reshape((__n,1)) + if Hybrid == "E3DVAR": + betaf = selfA._parameters["HybridCovarianceEquilibrium"] + Xn = Apply3DVarRecentringOnEnsemble(Xn, EMX, Ynpu, HO, R, B, betaf) + # + Xa = EnsembleMean( Xn ) #-------------------------- selfA._setInternalState("Xn", Xn) selfA._setInternalState("seed", numpy.random.get_state()) @@ -3165,7 +3250,7 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" or selfA._toStore("InnovationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentAnalysis") \ or selfA._toStore("SimulatedObservationAtCurrentOptimum"): - _HXa = numpy.asmatrix(numpy.ravel( H((Xa, Un)) )).T + _HXa = numpy.ravel( H((Xa, Un)) ).reshape((-1,1)) _Innovation = Ynpu - _HXa # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["Analysis"]) ) @@ -3197,8 +3282,8 @@ def senkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="KalmanFilterFormula16" 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 ) + 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 ) @@ -3255,16 +3340,13 @@ def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): # # Initialisations # --------------- - # - # Opérateurs Hm = HO["Direct"].appliedTo Ha = HO["Adjoint"].appliedInXTo # - # Utilisation éventuelle d'un vecteur H(Xb) précalculé if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]: - HXb = Hm( Xb, HO["AppliedInX"]["HXb"] ) + HXb = numpy.asarray(Hm( Xb, HO["AppliedInX"]["HXb"] )) else: - HXb = Hm( Xb ) + HXb = numpy.asarray(Hm( Xb )) HXb = HXb.reshape((-1,1)) if Y.size != HXb.size: raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size)) @@ -3276,22 +3358,20 @@ def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): HtMb = HtMb.reshape(Y.size,Xb.size) # ADAO & check shape selfA.StoredVariables["JacobianMatrixAtBackground"].store( HtMb ) # - # Précalcul des inversions de B et R BI = B.getI() RI = R.getI() # - # Point de démarrage de l'optimisation Xini = selfA._parameters["InitializationPoint"] # # Définition de la fonction-coût # ------------------------------ def CostFunction(x): - _X = numpy.ravel( x ).reshape((-1,1)) + _X = numpy.asarray(x).reshape((-1,1)) if selfA._parameters["StoreInternalVariables"] or \ selfA._toStore("CurrentState") or \ selfA._toStore("CurrentOptimum"): selfA.StoredVariables["CurrentState"].store( _X ) - _HX = Hm( _X ).reshape((-1,1)) + _HX = numpy.asarray(Hm( _X )).reshape((-1,1)) _Innovation = Y - _HX if selfA._toStore("SimulatedObservationAtCurrentState") or \ selfA._toStore("SimulatedObservationAtCurrentOptimum"): @@ -3329,8 +3409,8 @@ def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return J # def GradientOfCostFunction(x): - _X = x.reshape((-1,1)) - _HX = Hm( _X ).reshape((-1,1)) + _X = numpy.asarray(x).reshape((-1,1)) + _HX = numpy.asarray(Hm( _X )).reshape((-1,1)) GradJb = BI * (_X - Xb) GradJo = - Ha( (_X, RI * (Y - _HX)) ) GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo ) @@ -3503,18 +3583,18 @@ def std4dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): def Un(_step): if U is not None: if hasattr(U,"store") and 1<=_step1: - Un = numpy.asmatrix(numpy.ravel( U[step] )).T + Un = numpy.ravel( U[step] ).reshape((-1,1)) elif hasattr(U,"store") and len(U)==1: - Un = numpy.asmatrix(numpy.ravel( U[0] )).T + Un = numpy.ravel( U[0] ).reshape((-1,1)) else: - Un = numpy.asmatrix(numpy.ravel( U )).T + Un = numpy.ravel( U ).reshape((-1,1)) else: Un = None # if selfA._parameters["EstimationOf"] == "State": # Forecast + Q and observation of forecast - Xn_predicted = Mt * Xn + Xn_predicted = Mt @ Xn 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 + Xn_predicted = Xn_predicted + Cm @ Un Pn_predicted = Q + Mt * (Pn * Ma) elif selfA._parameters["EstimationOf"] == "Parameters": # Observation of forecast # --- > Par principe, M = Id, Q = 0 @@ -3804,13 +3885,13 @@ def stdkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): Pn_predicted = Pn # if selfA._parameters["EstimationOf"] == "State": - HX_predicted = Ht * Xn_predicted + HX_predicted = Ht @ Xn_predicted _Innovation = Ynpu - HX_predicted elif selfA._parameters["EstimationOf"] == "Parameters": - HX_predicted = Ht * Xn_predicted + HX_predicted = Ht @ Xn_predicted _Innovation = Ynpu - HX_predicted if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! - _Innovation = _Innovation - Cm * Un + _Innovation = _Innovation - Cm @ Un # Kn = Pn_predicted * Ha * numpy.linalg.inv(R + numpy.dot(Ht, Pn_predicted * Ha)) Xn = Xn_predicted + Kn * _Innovation @@ -3851,8 +3932,8 @@ def stdkf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): 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 ) + 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 ) @@ -3959,6 +4040,7 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): RI = R.getI() # __n = Xb.size + nbPreviousSteps = len(selfA.StoredVariables["Analysis"]) # if len(selfA.StoredVariables["Analysis"])==0 or not selfA._parameters["nextStep"]: Xn = Xb @@ -4004,7 +4086,7 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): XEtnnpi = numpy.asarray( Mm( (Xnp[:,point], Un) ) ).reshape((-1,1)) if Cm is not None and Un is not None: # Attention : si Cm est aussi dans M, doublon ! Cm = Cm.reshape(Xn.size,Un.size) # ADAO & check shape - XEtnnpi = XEtnnpi + Cm * Un + XEtnnpi = XEtnnpi + Cm @ Un elif selfA._parameters["EstimationOf"] == "Parameters": # --- > Par principe, M = Id, Q = 0 XEtnnpi = Xnp[:,point] @@ -4042,7 +4124,7 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): _Innovation = Ynpu - Yncm.reshape((-1,1)) if selfA._parameters["EstimationOf"] == "Parameters": if Cm is not None and Un is not None: # Attention : si Cm est aussi dans H, doublon ! - _Innovation = _Innovation - Cm * Un + _Innovation = _Innovation - Cm @ Un # Kn = Pxyn * Pyyn.I Xn = Xncm.reshape((-1,1)) + Kn * _Innovation @@ -4083,7 +4165,7 @@ def uskf(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): or selfA._toStore("CostFunctionJo") \ or selfA._toStore("CurrentOptimum") \ or selfA._toStore("APosterioriCovariance"): - Jb = float( 0.5 * (Xa - Xb).T * BI * (Xa - Xb) ) + 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 ) @@ -4138,29 +4220,24 @@ def van3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): # # Initialisations # --------------- - # - # Opérateurs Hm = HO["Direct"].appliedTo Ha = HO["Adjoint"].appliedInXTo # - # Précalcul des inversions de B et R BT = B.getT() RI = R.getI() # - # Point de démarrage de l'optimisation - Xini = numpy.zeros(Xb.shape) + Xini = numpy.zeros(Xb.size) # # Définition de la fonction-coût # ------------------------------ def CostFunction(v): - _V = numpy.asmatrix(numpy.ravel( v )).T - _X = Xb + B * _V + _V = numpy.asarray(v).reshape((-1,1)) + _X = Xb + (B @ _V).reshape((-1,1)) if selfA._parameters["StoreInternalVariables"] or \ selfA._toStore("CurrentState") or \ selfA._toStore("CurrentOptimum"): selfA.StoredVariables["CurrentState"].store( _X ) - _HX = Hm( _X ) - _HX = numpy.asmatrix(numpy.ravel( _HX )).T + _HX = numpy.asarray(Hm( _X )).reshape((-1,1)) _Innovation = Y - _HX if selfA._toStore("SimulatedObservationAtCurrentState") or \ selfA._toStore("SimulatedObservationAtCurrentOptimum"): @@ -4168,8 +4245,8 @@ def van3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): if selfA._toStore("InnovationAtCurrentState"): selfA.StoredVariables["InnovationAtCurrentState"].store( _Innovation ) # - Jb = float( 0.5 * _V.T * BT * _V ) - Jo = float( 0.5 * _Innovation.T * RI * _Innovation ) + Jb = float( 0.5 * _V.T * (BT * _V) ) + Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) ) J = Jb + Jo # selfA.StoredVariables["CurrentIterationNumber"].store( len(selfA.StoredVariables["CostFunctionJ"]) ) @@ -4198,9 +4275,9 @@ def van3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): return J # def GradientOfCostFunction(v): - _V = v.reshape((-1,1)) + _V = numpy.asarray(v).reshape((-1,1)) _X = Xb + (B @ _V).reshape((-1,1)) - _HX = Hm( _X ).reshape((-1,1)) + _HX = numpy.asarray(Hm( _X )).reshape((-1,1)) GradJb = BT * _V GradJo = - Ha( (_X, RI * (Y - _HX)) ) GradJ = numpy.ravel( GradJb ) + numpy.ravel( GradJo )