From dc2e5a3fe0f5d12041a360d1dd19dda0a8157b04 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Thu, 10 Feb 2022 21:29:07 +0100 Subject: [PATCH] Minor documentation and code review corrections (21) --- doc/en/snippets/CurrentIterationNumber.rst | 4 +- doc/en/snippets/Observation.rst | 11 +-- doc/fr/snippets/Observation.rst | 7 +- src/daComposant/daAlgorithms/3DVAR.py | 1 - src/daComposant/daCore/BasicObjects.py | 8 ++- src/daComposant/daCore/NumericObjects.py | 80 ++++++++++++---------- 6 files changed, 62 insertions(+), 49 deletions(-) diff --git a/doc/en/snippets/CurrentIterationNumber.rst b/doc/en/snippets/CurrentIterationNumber.rst index d1eea00..8409421 100644 --- a/doc/en/snippets/CurrentIterationNumber.rst +++ b/doc/en/snippets/CurrentIterationNumber.rst @@ -1,8 +1,8 @@ .. index:: single: CurrentIterationNumber CurrentIterationNumber - *List of integers*. Each element is the iteration index at the current step during the - iterative algorithm procedure. + *List of integers*. Each element is the iteration index at the current step + during the iterative algorithm procedure. Example: ``i = ADD.get("CurrentIterationNumber")[-1]`` diff --git a/doc/en/snippets/Observation.rst b/doc/en/snippets/Observation.rst index 53c095c..49518fd 100644 --- a/doc/en/snippets/Observation.rst +++ b/doc/en/snippets/Observation.rst @@ -1,8 +1,9 @@ .. index:: single: Observation Observation - *Vector*. The variable indicates the observation vector used for data - assimilation or optimization, usually noted as :math:`\mathbf{y}^o`. Its - value is defined as a "*Vector*" or "*VectorSerie*" type object. Its - availability in output is conditioned by the boolean "*Stored*" associated - with input. + *List of vectors*. The variable indicates the observation vector used for + data assimilation or optimization, and usually noted :math:`\mathbf{y}^o`. + Its value is defined as an object of type "*Vector*" if it is a single + observation (temporal or not) or "*VectorSeries*" if it is a succession of + observations. Its availability in output is conditioned by the boolean + "*Stored*" associated in input. diff --git a/doc/fr/snippets/Observation.rst b/doc/fr/snippets/Observation.rst index 81bcb3b..da59d99 100644 --- a/doc/fr/snippets/Observation.rst +++ b/doc/fr/snippets/Observation.rst @@ -1,8 +1,9 @@ .. index:: single: Observation Observation - *Vecteur*. La variable désigne le vecteur d'observation utilisé en + *Liste de vecteurs*. La variable désigne le vecteur d'observation utilisé en assimilation de données ou en optimisation, et usuellement noté :math:`\mathbf{y}^o`. Sa valeur est définie comme un objet de type "*Vector*" - ou "*VectorSerie*". Sa disponibilité en sortie est conditionnée par le - booléen "*Stored*" associé en entrée. + si c'est une unique observation (temporelle ou pas) ou "*VectorSerie*" si + c'est une succession d'observations. Sa disponibilité en sortie est + conditionnée par le booléen "*Stored*" associé en entrée. diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py index 28ae3bc..1500028 100644 --- a/src/daComposant/daAlgorithms/3DVAR.py +++ b/src/daComposant/daAlgorithms/3DVAR.py @@ -43,7 +43,6 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): listadv = [ "3DVAR-Std", "Incr3DVAR", - "OneCorrection3DVAR-Std", ], ) self.defineRequiredParameter( diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index 47b8c37..09f4bdb 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -625,6 +625,7 @@ class Algorithm(object): - CurrentIterationNumber : numéro courant d'itération dans les algorithmes itératifs, à partir de 0 - CurrentOptimum : état optimal courant lors d'itérations - CurrentState : état courant lors d'itérations + - CurrentStepNumber : numéro courant de pas de mesure dans les algorithmes temporels - GradientOfCostFunctionJ : gradient de la fonction-coût globale - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-coût - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-coût @@ -684,6 +685,7 @@ class Algorithm(object): self.StoredVariables["CurrentIterationNumber"] = Persistence.OneIndex(name = "CurrentIterationNumber") self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum") self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState") + self.StoredVariables["CurrentStepNumber"] = Persistence.OneIndex(name = "CurrentStepNumber") self.StoredVariables["ForecastCovariance"] = Persistence.OneMatrix(name = "ForecastCovariance") self.StoredVariables["ForecastState"] = Persistence.OneVector(name = "ForecastState") self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ") @@ -1094,6 +1096,7 @@ class PartialAlgorithm(object): self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb") self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo") self.StoredVariables["CurrentIterationNumber"] = Persistence.OneIndex(name = "CurrentIterationNumber") + self.StoredVariables["CurrentStepNumber"] = Persistence.OneIndex(name = "CurrentStepNumber") # self.__canonical_stored_name = {} for k in self.StoredVariables: @@ -2099,7 +2102,10 @@ class Covariance(object): return self.__C + numpy.asmatrix(other) elif self.isvector() or self.isscalar(): _A = numpy.asarray(other) - _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C + if len(_A.shape) == 1: + _A.reshape((-1,1))[::2] += self.__C + else: + _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C return numpy.asmatrix(_A) def __radd__(self, other): diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index c1feb22..d1fbf96 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -477,24 +477,24 @@ class FDApproximation(object): else: return _HaY # ============================================================================== -def EnsembleOfCenteredPerturbations( __bgcenter, __bgcovariance, __nbmembers ): - "Génération d'un ensemble de taille __nbmembers-1 d'états aléatoires centrés" +def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ): + "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés" # - __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),)) + __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) + 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 + _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 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 @@ -508,31 +508,31 @@ def EnsembleOfBackgroundPerturbations( __bgcenter, __bgcovariance, __nbmembers, 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) + __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: + _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) + 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 + _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 + 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) + _Perturbations = numpy.tile( __bgCenter, __nbMembers) # return _Perturbations @@ -552,9 +552,9 @@ def EnsembleOfAnomalies( __Ensemble, __OptMean = None, __Normalisation = 1.): return __Normalisation * (numpy.asarray( __Ensemble ) - __Em) # ============================================================================== -def EnsembleErrorCovariance( __Ensemble, __quick = False ): +def EnsembleErrorCovariance( __Ensemble, __Quick = False ): "Renvoie l'estimation empirique de la covariance d'ensemble" - if __quick: + if __Quick: # Covariance rapide mais rarement définie positive __Covariance = numpy.cov( __Ensemble ) else: @@ -638,6 +638,9 @@ def CovarianceInflation( else: __InflationFactor = float(__InflationFactor) # + __InputCovOrEns = numpy.asarray(__InputCovOrEns) + if __InputCovOrEns.size == 0: return __InputCovOrEns + # if __InflationType in ["MultiplicativeOnAnalysisCovariance", "MultiplicativeOnBackgroundCovariance"]: if __InflationFactor < 1.: raise ValueError("Inflation factor for multiplicative inflation has to be greater or equal than 1.") @@ -659,17 +662,20 @@ def CovarianceInflation( 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 + __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.identity(__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.: 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 + __n, __m = __InputCovOrEns.shape if __n != __m: raise ValueError("Additive inflation can only be applied to squared (covariance) matrix.") if __BackgroundCov is None: @@ -679,10 +685,10 @@ def CovarianceInflation( __OutputCovOrEns = (1. - __InflationFactor) * __InputCovOrEns + __InflationFactor * __BackgroundCov # elif __InflationType == "Relaxation": - raise NotImplementedError("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 @@ -853,7 +859,7 @@ def Apply3DVarRecentringOnEnsemble(__EnXn, __EnXf, __Ynpu, __HO, __R, __B, __Bet # ============================================================================== def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle): """ - Prévision multi-pas avec une correction par pas en X (multi-méthodes) + Prévision multi-pas avec une correction par pas (multi-méthodes) """ # # Initialisation -- 2.39.2