Internal quantiles estimations

diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py
index d47a32af197e15fe67fd6e75fcfd3912ac55066b..57fa0786eacfcdd01452d478e47d2fc9777321a3 100644 (file)
--- a/src/daComposant/daAlgorithms/3DVAR.py
+++ b/src/daComposant/daAlgorithms/3DVAR.py
@@ -84,7 +84,32 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = [],
             typecast = tuple,
             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
-            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaObs2", "MahalanobisConsistency"]
+            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaObs2", "MahalanobisConsistency", "SimulationQuantiles"]
+            )
+        self.defineRequiredParameter(
+            name     = "Quantiles",
+            default  = [],
+            typecast = tuple,
+            message  = "Liste des valeurs de quantiles",
+            )
+        self.defineRequiredParameter(
+            name     = "SetSeed",
+            typecast = numpy.random.seed,
+            message  = "Graine fixée pour le générateur aléatoire",
+            )
+        self.defineRequiredParameter(
+            name     = "NumberOfSamplesForQuantiles",
+            default  = 100,
+            typecast = int,
+            message  = "Nombre d'échantillons simulés pour le calcul des quantiles",
+            minval   = 1,
+            )
+        self.defineRequiredParameter(
+            name     = "SimulationForQuantiles",
+            default  = "Linear",
+            typecast = str,
+            message  = "Type de simulation pour l'estimation des quantiles",
+            listval  = ["Linear", "NonLinear"]
             )
 
     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
@@ -243,9 +268,15 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         self.StoredVariables["Analysis"].store( Xa.A1 )
         #
+        if  "OMA"                   in self._parameters["StoreSupplementaryCalculations"] or \
+            "SigmaObs2"             in self._parameters["StoreSupplementaryCalculations"] or \
+            "SimulationQuantiles"   in self._parameters["StoreSupplementaryCalculations"]:
+            HXa = Hm(Xa)
+        #
         # Calcul de la covariance d'analyse
         # ---------------------------------
-        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
+        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
             HtM = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
             HtM = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
             HaM = HO["Adjoint"].asMatrix(ValueForMethodForm = Xa)
@@ -278,14 +309,39 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         if "BMA" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
         if "OMA" in self._parameters["StoreSupplementaryCalculations"]:
-            self.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(Hm(Xa)) )
+            self.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
         if "OMB" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["OMB"].store( numpy.ravel(d) )
         if "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"]:
             TraceR = R.trace(Y.size)
-            self.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(Hm(Xa))).T)) ) / TraceR )
+            self.StoredVariables["SigmaObs2"].store( float( (d.T * (numpy.asmatrix(numpy.ravel(Y)).T-numpy.asmatrix(numpy.ravel(HXa)).T)) ) / TraceR )
         if "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["MahalanobisConsistency"].store( float( 2.*MinJ/d.size ) )
+        if "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
+            Qtls = self._parameters["Quantiles"]
+            nech = self._parameters["NumberOfSamplesForQuantiles"]
+            HXa  = numpy.matrix(numpy.ravel( HXa )).T
+            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( HtM * dXr )).T
+                    Yr = HXa + dYr
+                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
+                else:
+                    YfQ = numpy.hstack((YfQ,Yr))
+            YfQ.sort(axis=-1)
+            YQ = None
+            for quantile in Qtls:
+                if not (0. <= quantile <= 1.): continue
+                indice = int(nech * quantile - 1./nech)
+                if YQ is None: YQ = YfQ[:,indice]
+                else:          YQ = numpy.hstack((YQ,YfQ[:,indice]))
+            self.StoredVariables["SimulationQuantiles"].store( YQ )
         #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
         logging.debug("%s Terminé"%self._name)

diff --git a/src/daComposant/daAlgorithms/Blue.py b/src/daComposant/daAlgorithms/Blue.py
index 7cbe364679471bdd8ef89e0cf2af668a9079b3d4..d83be73512a1579526c82abda4fb91fc8ab3b3b2 100644 (file)
--- a/src/daComposant/daAlgorithms/Blue.py
+++ b/src/daComposant/daAlgorithms/Blue.py
@@ -41,7 +41,32 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = [],
             typecast = tuple,
             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
-            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaBck2", "SigmaObs2", "MahalanobisConsistency"]
+            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaBck2", "SigmaObs2", "MahalanobisConsistency", "SimulationQuantiles"]
+            )
+        self.defineRequiredParameter(
+            name     = "Quantiles",
+            default  = [],
+            typecast = tuple,
+            message  = "Liste des valeurs de quantiles",
+            )
+        self.defineRequiredParameter(
+            name     = "SetSeed",
+            typecast = numpy.random.seed,
+            message  = "Graine fixée pour le générateur aléatoire",
+            )
+        self.defineRequiredParameter(
+            name     = "NumberOfSamplesForQuantiles",
+            default  = 100,
+            typecast = int,
+            message  = "Nombre d'échantillons simulés pour le calcul des quantiles",
+            minval   = 1,
+            )
+        self.defineRequiredParameter(
+            name     = "SimulationForQuantiles",
+            default  = "Linear",
+            typecast = str,
+            message  = "Type de simulation pour l'estimation des quantiles",
+            listval  = ["Linear", "NonLinear"]
             )
 
     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
@@ -102,9 +127,15 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de la fonction coût
         # --------------------------
-        if self._parameters["StoreInternalVariables"] or "OMA" in self._parameters["StoreSupplementaryCalculations"] or "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"] or "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
-            oma = Y - Hm * Xa
-        if self._parameters["StoreInternalVariables"] or "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._parameters["StoreInternalVariables"] or \
+           "OMA" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"] or \
+           "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
+            HXa = Hm * Xa
+            oma = Y - HXa
+        if self._parameters["StoreInternalVariables"] or \
+           "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
             Jb  = 0.5 * (Xa - Xb).T * BI * (Xa - Xb)
             Jo  = 0.5 * oma.T * RI * oma
             J   = float( Jb ) + float( Jo )
@@ -114,7 +145,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de la covariance d'analyse
         # ---------------------------------
-        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
+        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
             A = B - K * Hm * B
             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."%(self._name,str(A.shape)))
@@ -144,6 +176,30 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             self.StoredVariables["SigmaBck2"].store( float( (d.T * Hm * (Xa - Xb))/(Hm * B * Hm.T).trace() ) )
         if "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["MahalanobisConsistency"].store( float( 2.*J/d.size ) )
+        if "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
+            Qtls = self._parameters["Quantiles"]
+            nech = self._parameters["NumberOfSamplesForQuantiles"]
+            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
+                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
+                else:
+                    YfQ = numpy.hstack((YfQ,Yr))
+            YfQ.sort(axis=-1)
+            YQ = None
+            for quantile in Qtls:
+                if not (0. <= quantile <= 1.): continue
+                indice = int(nech * quantile - 1./nech)
+                if YQ is None: YQ = YfQ[:,indice]
+                else:          YQ = numpy.hstack((YQ,YfQ[:,indice]))
+            self.StoredVariables["SimulationQuantiles"].store( YQ )
         #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
         logging.debug("%s Terminé"%self._name)

diff --git a/src/daComposant/daAlgorithms/ExtendedBlue.py b/src/daComposant/daAlgorithms/ExtendedBlue.py
index 196b526ee5fa24a8fc5bcaf992501f569ba2a3a6..ee14628412bdc368c893fa97a218be8eb6dd5212 100644 (file)
--- a/src/daComposant/daAlgorithms/ExtendedBlue.py
+++ b/src/daComposant/daAlgorithms/ExtendedBlue.py
@@ -41,7 +41,32 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = [],
             typecast = tuple,
             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
-            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaBck2", "SigmaObs2", "MahalanobisConsistency"]
+            listval  = ["APosterioriCovariance", "BMA", "OMA", "OMB", "Innovation", "SigmaBck2", "SigmaObs2", "MahalanobisConsistency", "SimulationQuantiles"]
+            )
+        self.defineRequiredParameter(
+            name     = "Quantiles",
+            default  = [],
+            typecast = tuple,
+            message  = "Liste des valeurs de quantiles",
+            )
+        self.defineRequiredParameter(
+            name     = "SetSeed",
+            typecast = numpy.random.seed,
+            message  = "Graine fixée pour le générateur aléatoire",
+            )
+        self.defineRequiredParameter(
+            name     = "NumberOfSamplesForQuantiles",
+            default  = 100,
+            typecast = int,
+            message  = "Nombre d'échantillons simulés pour le calcul des quantiles",
+            minval   = 1,
+            )
+        self.defineRequiredParameter(
+            name     = "SimulationForQuantiles",
+            default  = "Linear",
+            typecast = str,
+            message  = "Type de simulation pour l'estimation des quantiles",
+            listval  = ["Linear", "NonLinear"]
             )
 
     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
@@ -58,13 +83,13 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         Hm = Hm.reshape(Y.size,Xb.size) # ADAO & check shape
         Ha = HO["Adjoint"].asMatrix(Xb)
         Ha = Ha.reshape(Xb.size,Y.size) # ADAO & check shape
+        H  = HO["Direct"].appliedTo
         #
         # Utilisation éventuelle d'un vecteur H(Xb) précalculé
         # ----------------------------------------------------
         if HO["AppliedToX"] is not None and HO["AppliedToX"].has_key("HXb"):
             HXb = HO["AppliedToX"]["HXb"]
         else:
-            H  = HO["Direct"].appliedTo
             HXb = H( Xb )
         HXb = numpy.asmatrix(numpy.ravel( HXb )).T
         #
@@ -103,9 +128,15 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de la fonction coût
         # --------------------------
-        if self._parameters["StoreInternalVariables"] or "OMA" in self._parameters["StoreSupplementaryCalculations"] or "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"] or "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
-            oma = Y - Hm * Xa
-        if self._parameters["StoreInternalVariables"] or "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._parameters["StoreInternalVariables"] or \
+           "OMA" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"] or \
+           "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
+            HXa  = numpy.matrix(numpy.ravel( H( Xa ) )).T
+            oma = Y - HXa
+        if self._parameters["StoreInternalVariables"] or \
+           "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
             Jb  = 0.5 * (Xa - Xb).T * BI * (Xa - Xb)
             Jo  = 0.5 * oma.T * RI * oma
             J   = float( Jb ) + float( Jo )
@@ -115,7 +146,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de la covariance d'analyse
         # ---------------------------------
-        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
+        if "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"] or \
+           "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
             A = B - K * Hm * B
             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."%(self._name,str(A.shape)))
@@ -145,6 +177,32 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             self.StoredVariables["SigmaBck2"].store( float( (d.T * Hm * (Xa - Xb))/(Hm * B * Hm.T).trace() ) )
         if "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"]:
             self.StoredVariables["MahalanobisConsistency"].store( float( 2.*J/d.size ) )
+        if "SimulationQuantiles" in self._parameters["StoreSupplementaryCalculations"]:
+            Qtls = self._parameters["Quantiles"]
+            nech = self._parameters["NumberOfSamplesForQuantiles"]
+            HtM  = HO["Tangent"].asMatrix(ValueForMethodForm = Xa)
+            HtM  = HtM.reshape(Y.size,Xa.size) # ADAO & check shape
+            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( HtM * dXr )).T
+                    Yr = HXa + dYr
+                elif self._parameters["SimulationForQuantiles"] == "NonLinear":
+                    Xr = numpy.matrix(numpy.random.multivariate_normal(Xa.A1,A)).T
+                    Yr = numpy.matrix(numpy.ravel( H( Xr ) )).T
+                if YfQ is None:
+                    YfQ = Yr
+                else:
+                    YfQ = numpy.hstack((YfQ,Yr))
+            YfQ.sort(axis=-1)
+            YQ = None
+            for quantile in Qtls:
+                if not (0. <= quantile <= 1.): continue
+                indice = int(nech * quantile - 1./nech)
+                if YQ is None: YQ = YfQ[:,indice]
+                else:          YQ = numpy.hstack((YQ,YfQ[:,indice]))
+            self.StoredVariables["SimulationQuantiles"].store( YQ )
         #
         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M")))
         logging.debug("%s Terminé"%self._name)

diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py
index e774e4f5a2cdfb4e251c55ed0c4db0cd0f45bd3a..e4b6e90ddfacbea1c4f928e7f0fbb44cf64fc990 100644 (file)
--- a/src/daComposant/daCore/BasicObjects.py
+++ b/src/daComposant/daCore/BasicObjects.py
@@ -190,6 +190,7 @@ class Algorithm:
         self.StoredVariables["OMB"]                      = Persistence.OneVector(name = "OMB")
         self.StoredVariables["BMA"]                      = Persistence.OneVector(name = "BMA")
         self.StoredVariables["APosterioriCovariance"]    = Persistence.OneMatrix(name = "APosterioriCovariance")
+        self.StoredVariables["SimulationQuantiles"]      = Persistence.OneMatrix(name = "SimulationQuantiles")
 
     def get(self, key=None):
         """