Minor internal improvements

diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py
index df1dc90af9fb6be93bb27b4cb056fc62716f6776..20d6626bc7ee4d5042159f0f7f660cf8aaa5d22a 100644 (file)
--- 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 3f68c4ee800b2930aa748340914bbf6ef1072b93..e6a8fcd7c8c80b7ca643ae65469272f517871ced 100644 (file)
--- 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"):