X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FdaComposant%2FdaAlgorithms%2FExtendedBlue.py;h=ea93f62c7465c34a16a7d8cec7d3ee5bdf9dd396;hb=e0a87290afe3094c8daa93597c8bec55df0ba5b0;hp=c4085c5e70d82dc088462d7f65a834e3d59c8f8b;hpb=99bbc9e6f43597af2d5012cf4e4615d918015349;p=modules%2Fadao.git diff --git a/src/daComposant/daAlgorithms/ExtendedBlue.py b/src/daComposant/daAlgorithms/ExtendedBlue.py index c4085c5..ea93f62 100644 --- a/src/daComposant/daAlgorithms/ExtendedBlue.py +++ b/src/daComposant/daAlgorithms/ExtendedBlue.py @@ -1,6 +1,6 @@ #-*-coding:iso-8859-1-*- # -# Copyright (C) 2008-2013 EDF R&D +# Copyright (C) 2008-2015 EDF R&D # # This library is free software; you can redistribute it and/or # modify it under the terms of the GNU Lesser General Public @@ -21,9 +21,7 @@ # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D import logging -from daCore import BasicObjects, PlatformInfo -m = PlatformInfo.SystemUsage() - +from daCore import BasicObjects import numpy # ============================================================================== @@ -41,12 +39,38 @@ 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", "CostFunctionJ", "Innovation", "SigmaBck2", "SigmaObs2", "MahalanobisConsistency", "SimulationQuantiles", "SimulatedObservationAtBackground", "SimulatedObservationAtOptimum"] + ) + self.defineRequiredParameter( + name = "Quantiles", + default = [], + typecast = tuple, + message = "Liste des valeurs de quantiles", + minval = 0., + maxval = 1., + ) + 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): - logging.debug("%s Lancement"%self._name) - logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M"))) + self._pre_run() # # Paramètres de pilotage # ---------------------- @@ -54,11 +78,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm): # # Opérateur d'observation # ----------------------- - H = HO["Direct"].appliedTo Hm = HO["Tangent"].asMatrix(Xb) 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é # ---------------------------------------------------- @@ -103,9 +127,18 @@ 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 \ + "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] or \ + "OMA" in self._parameters["StoreSupplementaryCalculations"] or \ + "SigmaObs2" in self._parameters["StoreSupplementaryCalculations"] or \ + "MahalanobisConsistency" in self._parameters["StoreSupplementaryCalculations"] or \ + "SimulatedObservationAtOptimum" 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 \ + "CostFunctionJ" in self._parameters["StoreSupplementaryCalculations"] 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 +148,11 @@ 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"]: + if (Y.size <= Xb.size) and (Y.size > 100): K = B * Ha * (R + Hm * B * Ha).I + elif (Y.size > Xb.size) and (Y.size > 100): K = (BI + Ha * RI * Hm).I * Ha * RI + else: pass # K deja calcule 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,10 +182,38 @@ 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 ) + if "SimulatedObservationAtBackground" in self._parameters["StoreSupplementaryCalculations"]: + self.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) ) + if "SimulatedObservationAtOptimum" in self._parameters["StoreSupplementaryCalculations"]: + self.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) ) # - logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("M"))) - logging.debug("%s Terminé"%self._name) - # + self._post_run(HO) return 0 # ==============================================================================