]> SALOME platform Git repositories - modules/adao.git/commitdiff
Salome HOME
Ajout de l'algorithme NonLinearLeastSquares et correction des autres
authorJean-Philippe ARGAUD <ahbhhjp@cli23jp>
Fri, 3 Feb 2012 13:19:05 +0000 (14:19 +0100)
committerJean-Philippe ARGAUD <jean-philippe.argaud@edf.fr>
Fri, 3 Feb 2012 13:27:04 +0000 (14:27 +0100)
src/daComposant/daAlgorithms/3DVAR.py
src/daComposant/daAlgorithms/Blue.py
src/daComposant/daAlgorithms/NonLinearLeastSquares.py [new file with mode: 0644]

index bc17491790e7891d82b7bd29b755c63fc56ac7f8..48127504105879317a9069d90e2ceb6ce22c0ad7 100644 (file)
@@ -66,13 +66,24 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de l'innovation
         # ----------------------
+        if Y.size != HXb.size:
+            raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
+        if max(Y.shape) != max(HXb.shape):
+            raise ValueError("The shapes %s of observations Y and %s of observed calculation H(X) are different, they have to be identical."%(Y.shape,HXb.shape))
         d  = Y - HXb
         logging.debug("%s Innovation d = %s"%(self._name, d))
         #
         # Précalcul des inversion appellée dans les fonction-coût et gradient
         # -------------------------------------------------------------------
-        BI = B.I
-        RI = R.I
+        if B is not None:
+            BI = B.I
+        elif Parameters["B_scalar"] is not None:
+            BI = 1.0 / Parameters["B_scalar"]
+        #
+        if R is not None:
+            RI = R.I
+        elif Parameters["R_scalar"] is not None:
+            RI = 1.0 / Parameters["R_scalar"]
         #
         # Définition de la fonction-coût
         # ------------------------------
@@ -135,10 +146,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         logging.debug("%s Nombre maximal de pas d'optimisation = %s"%(self._name, str(maxiter)))
         if Parameters.has_key("CostDecrementTolerance") and (Parameters["CostDecrementTolerance"] > 0):
             ftol  = float(Parameters["CostDecrementTolerance"])
-            factr = 1./ftol
+            factr = ftol * 1.e14
         else:
             ftol  = 1.e-7
-            factr = 1./ftol
+            factr = ftol * 1.e14
         logging.debug("%s Diminution relative minimale du cout lors de l'arret = %s"%(self._name, str(1./factr)))
         if Parameters.has_key("ProjectedGradientTolerance") and (Parameters["ProjectedGradientTolerance"] > -1):
             pgtol = float(Parameters["ProjectedGradientTolerance"])
@@ -152,11 +163,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         logging.debug("%s Maximum des composantes du gradient lors de l'arret = %s"%(self._name, str(gtol)))
         InnerMinimizerList = ["CG", "NCG", "BFGS"]
         if Parameters.has_key("InnerMinimizer") and (Parameters["InnerMinimizer"] in InnerMinimizerList):
-            InnerMinimizer = str( Parameters["Minimizer"] )
+            InnerMinimizer = str( Parameters["InnerMinimizer"] )
         else:
             InnerMinimizer = "BFGS"
         logging.debug("%s Minimiseur interne utilisé = %s"%(self._name, InnerMinimizer))
-        logging.debug("%s Norme du gradient lors de l'arret = %s"%(self._name, str(gtol)))
         #
         # Minimisation de la fonctionnelle
         # --------------------------------
@@ -167,7 +177,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
                 fprime      = GradientOfCostFunction,
                 args        = (),
                 bounds      = Bounds,
-                maxfun      = maxiter,
+                maxfun      = maxiter-1,
                 factr       = factr,
                 pgtol       = pgtol,
                 iprint      = iprint,
index 0dcd6d44a4c75f08d6435e4ef4502f6748413937..ae15c23d57a6a868155d790a0e6e04a0a742d8e5 100644 (file)
@@ -63,6 +63,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Calcul de l'innovation et de l'analyse
         # --------------------------------------
+        if Y.size != HXb.size:
+            raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
+        if max(Y.shape) != max(HXb.shape):
+            raise ValueError("The shapes %s of observations Y and %s of observed calculation H(X) are different, they have to be identical."%(Y.shape,HXb.shape))
         d  = Y - HXb
         logging.debug("%s Innovation d = %s"%(self._name, d))
         Xa = Xb + K*d
diff --git a/src/daComposant/daAlgorithms/NonLinearLeastSquares.py b/src/daComposant/daAlgorithms/NonLinearLeastSquares.py
new file mode 100644 (file)
index 0000000..27528d4
--- /dev/null
@@ -0,0 +1,263 @@
+#-*-coding:iso-8859-1-*-
+#
+#  Copyright (C) 2008-2011  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
+#  License as published by the Free Software Foundation; either
+#  version 2.1 of the License.
+#
+#  This library is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+#  Lesser General Public License for more details.
+#
+#  You should have received a copy of the GNU Lesser General Public
+#  License along with this library; if not, write to the Free Software
+#  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+#
+#  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+
+import logging
+from daCore import BasicObjects, PlatformInfo
+m = PlatformInfo.SystemUsage()
+
+import numpy
+import scipy.optimize
+
+if logging.getLogger().level < 30:
+    iprint  = 1
+    message = scipy.optimize.tnc.MSG_ALL
+    disp    = 1
+else:
+    iprint  = -1
+    message = scipy.optimize.tnc.MSG_NONE
+    disp    = 0
+
+# ==============================================================================
+class ElementaryAlgorithm(BasicObjects.Algorithm):
+    def __init__(self):
+        BasicObjects.Algorithm.__init__(self)
+        self._name = "NONLINEARLEASTSQUARES"
+        logging.debug("%s Initialisation"%self._name)
+
+    def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
+        """
+        Calcul de l'estimateur moindres carrés pondérés non linéaires
+        (assimilation variationnelle sans ébauche)
+        """
+        logging.debug("%s Lancement"%self._name)
+        logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
+        #
+        # Opérateur d'observation
+        # -----------------------
+        Hm = H["Direct"].appliedTo
+        Ht = H["Adjoint"].appliedInXTo
+        #
+        # Utilisation éventuelle d'un vecteur H(Xb) précalculé
+        # ----------------------------------------------------
+        if H["AppliedToX"] is not None and H["AppliedToX"].has_key("HXb"):
+            logging.debug("%s Utilisation de HXb"%self._name)
+            HXb = H["AppliedToX"]["HXb"]
+        else:
+            logging.debug("%s Calcul de Hm(Xb)"%self._name)
+            HXb = Hm( Xb )
+        HXb = numpy.asmatrix(HXb).flatten().T
+        #
+        # Calcul de l'innovation
+        # ----------------------
+        if Y.size != HXb.size:
+            raise ValueError("The size %i of observations Y and %i of observed calculation H(X) are different, they have to be identical."%(Y.size,HXb.size))
+        if max(Y.shape) != max(HXb.shape):
+            raise ValueError("The shapes %s of observations Y and %s of observed calculation H(X) are different, they have to be identical."%(Y.shape,HXb.shape))
+        d  = Y - HXb
+        logging.debug("%s Innovation d = %s"%(self._name, d))
+        #
+        # Précalcul des inversion appellée dans les fonction-coût et gradient
+        # -------------------------------------------------------------------
+        # if B is not None:
+        #     BI = B.I
+        # elif Parameters["B_scalar"] is not None:
+        #     BI = 1.0 / Parameters["B_scalar"]
+        #
+        if R is not None:
+            RI = R.I
+        elif Parameters["R_scalar"] is not None:
+            RI = 1.0 / Parameters["R_scalar"]
+        #
+        # Définition de la fonction-coût
+        # ------------------------------
+        def CostFunction(x):
+            _X  = numpy.asmatrix(x).flatten().T
+            logging.info("%s CostFunction X  = %s"%(self._name, numpy.asmatrix( _X ).flatten()))
+            _HX = Hm( _X )
+            _HX = numpy.asmatrix(_HX).flatten().T
+            Jb  = 0.
+            Jo  = 0.5 * (Y - _HX).T * RI * (Y - _HX)
+            J   = float( Jb ) + float( Jo )
+            logging.info("%s CostFunction Jb = %s"%(self._name, Jb))
+            logging.info("%s CostFunction Jo = %s"%(self._name, Jo))
+            logging.info("%s CostFunction J  = %s"%(self._name, J))
+            self.StoredVariables["CurrentState"].store( _X.A1 )
+            self.StoredVariables["CostFunctionJb"].store( Jb )
+            self.StoredVariables["CostFunctionJo"].store( Jo )
+            self.StoredVariables["CostFunctionJ" ].store( J )
+            return float( J )
+        #
+        def GradientOfCostFunction(x):
+            _X      = numpy.asmatrix(x).flatten().T
+            logging.info("%s GradientOfCostFunction X      = %s"%(self._name, numpy.asmatrix( _X ).flatten()))
+            _HX     = Hm( _X )
+            _HX     = numpy.asmatrix(_HX).flatten().T
+            GradJb  = 0.
+            GradJo  = - Ht( (_X, RI * (Y - _HX)) )
+            GradJ   = numpy.asmatrix( GradJb ).flatten().T + numpy.asmatrix( GradJo ).flatten().T
+            logging.debug("%s GradientOfCostFunction GradJb = %s"%(self._name, numpy.asmatrix( GradJb ).flatten()))
+            logging.debug("%s GradientOfCostFunction GradJo = %s"%(self._name, numpy.asmatrix( GradJo ).flatten()))
+            logging.debug("%s GradientOfCostFunction GradJ  = %s"%(self._name, numpy.asmatrix( GradJ  ).flatten()))
+            return GradJ.A1
+        #
+        # Point de démarrage de l'optimisation : Xini = Xb
+        # ------------------------------------
+        if type(Xb) is type(numpy.matrix([])):
+            Xini = Xb.A1.tolist()
+        else:
+            Xini = list(Xb)
+        logging.debug("%s Point de démarrage Xini = %s"%(self._name, Xini))
+        #
+        # Paramètres de pilotage
+        # ----------------------
+        # Potentiels : "Bounds", "Minimizer", "MaximumNumberOfSteps", "ProjectedGradientTolerance", "GradientNormTolerance", "InnerMinimizer"
+        if Parameters.has_key("Bounds") and (type(Parameters["Bounds"]) is type([]) or type(Parameters["Bounds"]) is type(())) and (len(Parameters["Bounds"]) > 0):
+            Bounds = Parameters["Bounds"]
+        else:
+            Bounds = None
+        MinimizerList = ["LBFGSB","TNC", "CG", "NCG", "BFGS"]
+        if Parameters.has_key("Minimizer") and (Parameters["Minimizer"] in MinimizerList):
+            Minimizer = str( Parameters["Minimizer"] )
+        else:
+            Minimizer = "LBFGSB"
+            logging.warning("%s Unknown or undefined minimizer, replaced by the default one \"%s\""%(self._name,Minimizer))
+        logging.debug("%s Minimiseur utilisé = %s"%(self._name, Minimizer))
+        if Parameters.has_key("MaximumNumberOfSteps") and (Parameters["MaximumNumberOfSteps"] > -1):
+            maxiter = int( Parameters["MaximumNumberOfSteps"] )
+        else:
+            maxiter = 15000
+        logging.debug("%s Nombre maximal de pas d'optimisation = %s"%(self._name, str(maxiter)))
+        if Parameters.has_key("CostDecrementTolerance") and (Parameters["CostDecrementTolerance"] > 0):
+            ftol  = float(Parameters["CostDecrementTolerance"])
+            factr = ftol * 1.e14
+        else:
+            ftol  = 1.e-7
+            factr = ftol * 1.e14
+        logging.debug("%s Diminution relative minimale du cout lors de l'arret = %s"%(self._name, str(1./factr)))
+        if Parameters.has_key("ProjectedGradientTolerance") and (Parameters["ProjectedGradientTolerance"] > -1):
+            pgtol = float(Parameters["ProjectedGradientTolerance"])
+        else:
+            pgtol = -1
+        logging.debug("%s Maximum des composantes du gradient projete lors de l'arret = %s"%(self._name, str(pgtol)))
+        if Parameters.has_key("GradientNormTolerance") and (Parameters["GradientNormTolerance"] > -1):
+            gtol = float(Parameters["GradientNormTolerance"])
+        else:
+            gtol = 1.e-05
+        logging.debug("%s Maximum des composantes du gradient lors de l'arret = %s"%(self._name, str(gtol)))
+        InnerMinimizerList = ["CG", "NCG", "BFGS"]
+        if Parameters.has_key("InnerMinimizer") and (Parameters["InnerMinimizer"] in InnerMinimizerList):
+            InnerMinimizer = str( Parameters["InnerMinimizer"] )
+        else:
+            InnerMinimizer = "BFGS"
+        logging.debug("%s Minimiseur interne utilisé = %s"%(self._name, InnerMinimizer))
+        #
+        # Minimisation de la fonctionnelle
+        # --------------------------------
+        if Minimizer == "LBFGSB":
+            Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
+                func        = CostFunction,
+                x0          = Xini,
+                fprime      = GradientOfCostFunction,
+                args        = (),
+                bounds      = Bounds,
+                maxfun      = maxiter-1,
+                factr       = factr,
+                pgtol       = pgtol,
+                iprint      = iprint,
+                )
+            nfeval = Informations['funcalls']
+            rc     = Informations['warnflag']
+        elif Minimizer == "TNC":
+            Minimum, nfeval, rc = scipy.optimize.fmin_tnc(
+                func        = CostFunction,
+                x0          = Xini,
+                fprime      = GradientOfCostFunction,
+                args        = (),
+                bounds      = Bounds,
+                maxfun      = maxiter,
+                pgtol       = pgtol,
+                ftol        = ftol,
+                messages    = message,
+                )
+        elif Minimizer == "CG":
+            Minimum, fopt, nfeval, grad_calls, rc = scipy.optimize.fmin_cg(
+                f           = CostFunction,
+                x0          = Xini,
+                fprime      = GradientOfCostFunction,
+                args        = (),
+                maxiter     = maxiter,
+                gtol        = gtol,
+                disp        = disp,
+                full_output = True,
+                )
+        elif Minimizer == "NCG":
+            Minimum, fopt, nfeval, grad_calls, hcalls, rc = scipy.optimize.fmin_ncg(
+                f           = CostFunction,
+                x0          = Xini,
+                fprime      = GradientOfCostFunction,
+                args        = (),
+                maxiter     = maxiter,
+                avextol     = ftol,
+                disp        = disp,
+                full_output = True,
+                )
+        elif Minimizer == "BFGS":
+            Minimum, fopt, gopt, Hopt, nfeval, grad_calls, rc = scipy.optimize.fmin_bfgs(
+                f           = CostFunction,
+                x0          = Xini,
+                fprime      = GradientOfCostFunction,
+                args        = (),
+                maxiter     = maxiter,
+                gtol        = gtol,
+                disp        = disp,
+                full_output = True,
+                )
+        else:
+            raise ValueError("Error in Minimizer name: %s"%Minimizer)
+        #
+        # Correction pour pallier a un bug de TNC sur le retour du Minimum
+        # ----------------------------------------------------------------
+        StepMin = numpy.argmin( self.StoredVariables["CostFunctionJ"].valueserie() )
+        MinJ    = self.StoredVariables["CostFunctionJ"].valueserie(step = StepMin)
+        Minimum = self.StoredVariables["CurrentState"].valueserie(step = StepMin)
+        #
+        logging.debug("%s %s Step of min cost  = %s"%(self._name, Minimizer, StepMin))
+        logging.debug("%s %s Minimum cost      = %s"%(self._name, Minimizer, MinJ))
+        logging.debug("%s %s Minimum state     = %s"%(self._name, Minimizer, Minimum))
+        logging.debug("%s %s Nb of F           = %s"%(self._name, Minimizer, nfeval))
+        logging.debug("%s %s RetCode           = %s"%(self._name, Minimizer, rc))
+        #
+        # Calcul  de l'analyse
+        # --------------------
+        Xa = numpy.asmatrix(Minimum).T
+        logging.debug("%s Analyse Xa = %s"%(self._name, Xa))
+        #
+        self.StoredVariables["Analysis"].store( Xa.A1 )
+        self.StoredVariables["Innovation"].store( d.A1 )
+        #
+        logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
+        logging.debug("%s Terminé"%self._name)
+        #
+        return 0
+
+# ==============================================================================
+if __name__ == "__main__":
+    print '\n AUTODIAGNOSTIC \n'