Salome HOME
Minor internal modifications and bounds corrections
[modules/adao.git] / src / daComposant / daAlgorithms / NonLinearLeastSquares.py
index a2e4433ea88b2ca8697ee5d7873e55cf2e15557f..aad66fc737b92a4dac264863d8522024f8a2e278 100644 (file)
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 #
-# Copyright (C) 2008-2017 EDF R&D
+# Copyright (C) 2008-2021 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
@@ -22,7 +22,7 @@
 
 import logging
 from daCore import BasicObjects
-import numpy, scipy.optimize
+import numpy, scipy.optimize, scipy.version
 
 # ==============================================================================
 class ElementaryAlgorithm(BasicObjects.Algorithm):
@@ -46,7 +46,8 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             name     = "CostDecrementTolerance",
             default  = 1.e-7,
             typecast = float,
-            message  = "Diminution relative minimale du cout lors de l'arrêt",
+            message  = "Diminution relative minimale du coût lors de l'arrêt",
+            minval   = 0.,
             )
         self.defineRequiredParameter(
             name     = "ProjectedGradientTolerance",
@@ -60,6 +61,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = 1.e-05,
             typecast = float,
             message  = "Maximum des composantes du gradient lors de l'arrêt",
+            minval   = 0.,
             )
         self.defineRequiredParameter(
             name     = "StoreInternalVariables",
@@ -72,19 +74,49 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             default  = [],
             typecast = tuple,
             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
-            listval  = ["BMA", "OMA", "OMB", "CostFunctionJ", "CostFunctionJb", "CostFunctionJo", "CurrentState", "CurrentOptimum", "IndexOfOptimum", "Innovation", "InnovationAtCurrentState", "CostFunctionJAtCurrentOptimum", "SimulatedObservationAtBackground", "SimulatedObservationAtCurrentState", "SimulatedObservationAtOptimum", "SimulatedObservationAtCurrentOptimum"]
+            listval  = [
+                "Analysis",
+                "BMA",
+                "CostFunctionJ",
+                "CostFunctionJAtCurrentOptimum",
+                "CostFunctionJb",
+                "CostFunctionJbAtCurrentOptimum",
+                "CostFunctionJo",
+                "CostFunctionJoAtCurrentOptimum",
+                "CurrentIterationNumber",
+                "CurrentOptimum",
+                "CurrentState",
+                "IndexOfOptimum",
+                "Innovation",
+                "InnovationAtCurrentState",
+                "OMA",
+                "OMB",
+                "SimulatedObservationAtBackground",
+                "SimulatedObservationAtCurrentOptimum",
+                "SimulatedObservationAtCurrentState",
+                "SimulatedObservationAtOptimum",
+                ]
             )
         self.defineRequiredParameter( # Pas de type
             name     = "Bounds",
             message  = "Liste des valeurs de bornes",
             )
+        self.defineRequiredParameter(
+            name     = "InitializationPoint",
+            typecast = numpy.ravel,
+            message  = "État initial imposé (par défaut, c'est l'ébauche si None)",
+            )
+        self.requireInputArguments(
+            mandatory= ("Xb", "Y", "HO", "R"),
+            )
+        self.setAttributes(tags=(
+            "Optimization",
+            "NonLinear",
+            "Variational",
+            ))
 
     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
-        self._pre_run(Parameters)
-        #
-        # Correction pour pallier a un bug de TNC sur le retour du Minimum
-        if "Minimizer" in self._parameters and self._parameters["Minimizer"] == "TNC":
-            self.setParameterValue("StoreInternalVariables",True)
+        self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
         #
         # Opérateurs
         # ----------
@@ -94,7 +126,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         # Utilisation éventuelle d'un vecteur H(Xb) précalculé
         # ----------------------------------------------------
         if HO["AppliedInX"] is not None and "HXb" in HO["AppliedInX"]:
-            HXb = Hm( Xb, HO["AppliedInX"]["HXb"])
+            HXb = Hm( Xb, HO["AppliedInX"]["HXb"] )
         else:
             HXb = Hm( Xb )
         HXb = numpy.asmatrix(numpy.ravel( HXb )).T
@@ -114,39 +146,44 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         def CostFunction(x):
             _X  = numpy.asmatrix(numpy.ravel( x )).T
             if self._parameters["StoreInternalVariables"] or \
-                "CurrentState" in self._parameters["StoreSupplementaryCalculations"] or \
-                "CurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+                self._toStore("CurrentState") or \
+                self._toStore("CurrentOptimum"):
                 self.StoredVariables["CurrentState"].store( _X )
             _HX = Hm( _X )
             _HX = numpy.asmatrix(numpy.ravel( _HX )).T
             _Innovation = Y - _HX
-            if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"] or \
-               "SimulatedObservationAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("SimulatedObservationAtCurrentState") or \
+                self._toStore("SimulatedObservationAtCurrentOptimum"):
                 self.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
-            if "InnovationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("InnovationAtCurrentState"):
                 self.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
             #
             Jb  = 0.
             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
             J   = Jb + Jo
             #
+            self.StoredVariables["CurrentIterationNumber"].store( len(self.StoredVariables["CostFunctionJ"]) )
             self.StoredVariables["CostFunctionJb"].store( Jb )
             self.StoredVariables["CostFunctionJo"].store( Jo )
             self.StoredVariables["CostFunctionJ" ].store( J )
-            if "IndexOfOptimum" in self._parameters["StoreSupplementaryCalculations"] or \
-               "CurrentOptimum" in self._parameters["StoreSupplementaryCalculations"] or \
-               "CostFunctionJAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"] or \
-               "SimulatedObservationAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("IndexOfOptimum") or \
+                self._toStore("CurrentOptimum") or \
+                self._toStore("CostFunctionJAtCurrentOptimum") or \
+                self._toStore("CostFunctionJbAtCurrentOptimum") or \
+                self._toStore("CostFunctionJoAtCurrentOptimum") or \
+                self._toStore("SimulatedObservationAtCurrentOptimum"):
                 IndexMin = numpy.argmin( self.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
-            if "IndexOfOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("IndexOfOptimum"):
                 self.StoredVariables["IndexOfOptimum"].store( IndexMin )
-            if "CurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("CurrentOptimum"):
                 self.StoredVariables["CurrentOptimum"].store( self.StoredVariables["CurrentState"][IndexMin] )
-            if "SimulatedObservationAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("SimulatedObservationAtCurrentOptimum"):
                 self.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( self.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
-            if "CostFunctionJAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            if self._toStore("CostFunctionJbAtCurrentOptimum"):
                 self.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( self.StoredVariables["CostFunctionJb"][IndexMin] )
+            if self._toStore("CostFunctionJoAtCurrentOptimum"):
                 self.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( self.StoredVariables["CostFunctionJo"][IndexMin] )
+            if self._toStore("CostFunctionJAtCurrentOptimum"):
                 self.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( self.StoredVariables["CostFunctionJ" ][IndexMin] )
             return J
         #
@@ -168,7 +205,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             Jo  = float( 0.5 * _Innovation.T * RI * _Innovation )
             J   = Jb + Jo
             if self._parameters["StoreInternalVariables"] or \
-                "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
+                self._toStore("CurrentState"):
                 self.StoredVariables["CurrentState"].store( _X )
             self.StoredVariables["CostFunctionJb"].store( Jb )
             self.StoredVariables["CostFunctionJo"].store( Jo )
@@ -187,14 +224,19 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Point de démarrage de l'optimisation : Xini = Xb
         # ------------------------------------
-        Xini = numpy.ravel(Xb)
+        Xini = self._parameters["InitializationPoint"]
         #
         # Minimisation de la fonctionnelle
         # --------------------------------
         nbPreviousSteps = self.StoredVariables["CostFunctionJ"].stepnumber()
         #
         if self._parameters["Minimizer"] == "LBFGSB":
-            Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
+            # Minimum, J_optimal, Informations = scipy.optimize.fmin_l_bfgs_b(
+            if "0.19" <= scipy.version.version <= "1.1.0":
+                import lbfgsbhlt as optimiseur
+            else:
+                import scipy.optimize as optimiseur
+            Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b(
                 func        = CostFunction,
                 x0          = Xini,
                 fprime      = GradientOfCostFunction,
@@ -272,7 +314,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         # Correction pour pallier a un bug de TNC sur le retour du Minimum
         # ----------------------------------------------------------------
-        if self._parameters["StoreInternalVariables"] or "CurrentState" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._parameters["StoreInternalVariables"] or self._toStore("CurrentState"):
             Minimum = self.StoredVariables["CurrentState"][IndexMin]
         #
         # Obtention de l'analyse
@@ -281,32 +323,32 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         #
         self.StoredVariables["Analysis"].store( Xa.A1 )
         #
-        if "OMA"                           in self._parameters["StoreSupplementaryCalculations"] or \
-           "SimulatedObservationAtOptimum" in self._parameters["StoreSupplementaryCalculations"]:
-            if "SimulatedObservationAtCurrentState" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("OMA") or \
+            self._toStore("SimulatedObservationAtOptimum"):
+            if self._toStore("SimulatedObservationAtCurrentState"):
                 HXa = self.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
-            elif "SimulatedObservationAtCurrentOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+            elif self._toStore("SimulatedObservationAtCurrentOptimum"):
                 HXa = self.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
             else:
-                HXa = Hm(Xa)
+                HXa = Hm( Xa )
         #
         #
         # Calculs et/ou stockages supplémentaires
         # ---------------------------------------
-        if "Innovation" in self._parameters["StoreSupplementaryCalculations"] or \
-            "OMB" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("Innovation") or \
+            self._toStore("OMB"):
             d  = Y - HXb
-        if "Innovation" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("Innovation"):
             self.StoredVariables["Innovation"].store( numpy.ravel(d) )
-        if "BMA" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("BMA"):
             self.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
-        if "OMA" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("OMA"):
             self.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
-        if "OMB" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("OMB"):
             self.StoredVariables["OMB"].store( numpy.ravel(d) )
-        if "SimulatedObservationAtBackground" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("SimulatedObservationAtBackground"):
             self.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(HXb) )
-        if "SimulatedObservationAtOptimum" in self._parameters["StoreSupplementaryCalculations"]:
+        if self._toStore("SimulatedObservationAtOptimum"):
             self.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
         #
         self._post_run(HO)
@@ -314,4 +356,4 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
 
 # ==============================================================================
 if __name__ == "__main__":
-    print('\n AUTODIAGNOSTIC \n')
+    print('\n AUTODIAGNOSTIC\n')