From 3f84140cee5261d3ede72b355f41897ec54201a1 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Wed, 14 Apr 2021 07:16:43 +0200 Subject: [PATCH] Minor internal improvements --- src/daComposant/daAlgorithms/Blue.py | 32 +---------- src/daComposant/daCore/NumericObjects.py | 70 ++++++++++++++---------- 2 files changed, 44 insertions(+), 58 deletions(-) diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py index df1dc90..20d6626 100644 --- a/src/daComposant/daAlgorithms/Blue.py +++ b/src/daComposant/daAlgorithms/Blue.py @@ -21,7 +21,7 @@ # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D import logging -from daCore import BasicObjects +from daCore import BasicObjects, NumericObjects import numpy # ============================================================================== @@ -215,34 +215,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): if self._toStore("MahalanobisConsistency"): self.StoredVariables["MahalanobisConsistency"].store( float( 2.*J/d.size ) ) if self._toStore("SimulationQuantiles"): - nech = self._parameters["NumberOfSamplesForQuantiles"] - EXr = None - YfQ = None - for i in range(nech): - if self._parameters["SimulationForQuantiles"] == "Linear": - dXr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A) - Xa.A1).T - dYr = numpy.matrix(numpy.ravel( Hm * dXr )).T - Yr = HXa + dYr - if self._toStore("SampledStateForQuantiles"): Xr = Xa+dXr - elif self._parameters["SimulationForQuantiles"] == "NonLinear": - Xr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A)).T - Yr = numpy.matrix(numpy.ravel( Hm * Xr )).T - if YfQ is None: - YfQ = Yr - if self._toStore("SampledStateForQuantiles"): EXr = numpy.ravel(Xr) - else: - YfQ = numpy.hstack((YfQ,Yr)) - if self._toStore("SampledStateForQuantiles"): EXr = numpy.vstack((EXr,numpy.ravel(Xr))) - YfQ.sort(axis=-1) - YQ = None - for quantile in self._parameters["Quantiles"]: - if not (0. <= float(quantile) <= 1.): continue - indice = int(nech * float(quantile) - 1./nech) - if YQ is None: YQ = YfQ[:,indice] - else: YQ = numpy.hstack((YQ,YfQ[:,indice])) - self.StoredVariables["SimulationQuantiles"].store( YQ ) - if self._toStore("SampledStateForQuantiles"): - self.StoredVariables["SampledStateForQuantiles"].store( EXr.T ) + H = HO["Direct"].appliedTo + NumericObjects.QuantilesEstimations(self, A, Xa, HXa, H, Hm) if self._toStore("SimulatedObservationAtBackground"): self.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) ) if self._toStore("SimulatedObservationAtCurrentState"): diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index 3f68c4e..e6a8fcd 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -644,6 +644,46 @@ def CovarianceInflation( # return OutputCovOrEns +# ============================================================================== +def QuantilesEstimations(selfA, A, Xa, HXa = None, Hm = None, HtM = None): + "Estimation des quantiles a posteriori (selfA est modifié)" + nbsamples = selfA._parameters["NumberOfSamplesForQuantiles"] + # + # Échantillonnage des états + YfQ = None + EXr = None + if selfA._parameters["SimulationForQuantiles"] == "Linear": + HXa = numpy.matrix(numpy.ravel( HXa )).T + for i in range(nbsamples): + if selfA._parameters["SimulationForQuantiles"] == "Linear" and HtM is not None: + dXr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A) - Xa.A1).T + dYr = numpy.matrix(numpy.ravel( HtM * dXr )).T + Yr = HXa + dYr + if selfA._toStore("SampledStateForQuantiles"): Xr = Xa+dXr + elif selfA._parameters["SimulationForQuantiles"] == "NonLinear" and Hm is not None: + Xr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A)).T + Yr = numpy.matrix(numpy.ravel( Hm( Xr ) )).T + if YfQ is None: + YfQ = Yr + if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.ravel(Xr) + else: + YfQ = numpy.hstack((YfQ,Yr)) + if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.vstack((EXr,numpy.ravel(Xr))) + # + # 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] + else: YQ = numpy.hstack((YQ,YfQ[:,indice])) + selfA.StoredVariables["SimulationQuantiles"].store( YQ ) + if selfA._toStore("SampledStateForQuantiles"): + selfA.StoredVariables["SampledStateForQuantiles"].store( EXr.T ) + # + return 0 + # ============================================================================== def enks(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, VariantM="EnKS16-KalmanFilterFormula"): """ @@ -2784,35 +2824,7 @@ def std3dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): if selfA._toStore("MahalanobisConsistency"): selfA.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) ) if selfA._toStore("SimulationQuantiles"): - nech = selfA._parameters["NumberOfSamplesForQuantiles"] - HXa = numpy.matrix(numpy.ravel( HXa )).T - EXr = None - YfQ = None - for i in range(nech): - if selfA._parameters["SimulationForQuantiles"] == "Linear": - dXr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A) - Xa.A1).T - dYr = numpy.matrix(numpy.ravel( HtM * dXr )).T - Yr = HXa + dYr - if selfA._toStore("SampledStateForQuantiles"): Xr = Xa+dXr - elif selfA._parameters["SimulationForQuantiles"] == "NonLinear": - Xr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A)).T - Yr = numpy.matrix(numpy.ravel( Hm( Xr ) )).T - if YfQ is None: - YfQ = Yr - if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.ravel(Xr) - else: - YfQ = numpy.hstack((YfQ,Yr)) - if selfA._toStore("SampledStateForQuantiles"): EXr = numpy.vstack((EXr,numpy.ravel(Xr))) - YfQ.sort(axis=-1) - YQ = None - for quantile in selfA._parameters["Quantiles"]: - if not (0. <= float(quantile) <= 1.): continue - indice = int(nech * float(quantile) - 1./nech) - if YQ is None: YQ = YfQ[:,indice] - else: YQ = numpy.hstack((YQ,YfQ[:,indice])) - selfA.StoredVariables["SimulationQuantiles"].store( YQ ) - if selfA._toStore("SampledStateForQuantiles"): - selfA.StoredVariables["SampledStateForQuantiles"].store( EXr.T ) + QuantilesEstimations(selfA, A, Xa, HXa, Hm, HtM) if selfA._toStore("SimulatedObservationAtBackground"): selfA.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) ) if selfA._toStore("SimulatedObservationAtOptimum"): -- 2.39.2