From 077902bbf06af7094230531679186ea388d804de Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Sat, 23 Nov 2013 19:00:02 +0100 Subject: [PATCH] Internal quantiles estimations --- src/daComposant/daAlgorithms/3DVAR.py | 64 ++++++++++++++++-- src/daComposant/daAlgorithms/Blue.py | 66 ++++++++++++++++-- src/daComposant/daAlgorithms/ExtendedBlue.py | 70 ++++++++++++++++++-- src/daComposant/daCore/BasicObjects.py | 1 + 4 files changed, 186 insertions(+), 15 deletions(-) diff --git a/src/daComposant/daAlgorithms/3DVAR.py b/src/daComposant/daAlgorithms/3DVAR.py index d47a32a..57fa078 100644 --- 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 7cbe364..d83be73 100644 --- 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 196b526..ee14628 100644 --- 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 e774e4f..e4b6e90 100644 --- 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): """ -- 2.39.2