src/daComposant/daAlgorithms/DifferentialEvolution.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2020 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 import logging
  24 from daCore import BasicObjects
  25 import numpy, scipy.optimize
  26
  27 # ==============================================================================
  28 class ElementaryAlgorithm(BasicObjects.Algorithm):
  29     def __init__(self):
  30         BasicObjects.Algorithm.__init__(self, "DIFFERENTIALEVOLUTION")
  31         self.defineRequiredParameter(
  32             name     = "Minimizer",
  33             default  = "BEST1BIN",
  34             typecast = str,
  35             message  = "Stratégie de minimisation utilisée",
  36             listval  = [
  37                 "BEST1BIN",
  38                 "BEST1EXP",
  39                 "RAND1EXP",
  40                 "RANDTOBEST1EXP",
  41                 # "CURRENTTOBEST1EXP",
  42                 "BEST2EXP",
  43                 "RAND2EXP",
  44                 "RANDTOBEST1BIN",
  45                 # "CURRENTTOBEST1BIN",
  46                 "BEST2BIN",
  47                 "RAND2BIN",
  48                 "RAND1BIN",
  49                 ],
  50             )
  51         self.defineRequiredParameter(
  52             name     = "MaximumNumberOfSteps",
  53             default  = 15000,
  54             typecast = int,
  55             message  = "Nombre maximal de générations",
  56             minval   = 0,
  57             )
  58         self.defineRequiredParameter(
  59             name     = "MaximumNumberOfFunctionEvaluations",
  60             default  = 15000,
  61             typecast = int,
  62             message  = "Nombre maximal d'évaluations de la fonction",
  63             minval   = -1,
  64             )
  65         self.defineRequiredParameter(
  66             name     = "PopulationSize",
  67             default  = 100,
  68             typecast = int,
  69             message  = "Taille approximative de la population à chaque génération",
  70             minval   = 1,
  71             )
  72         self.defineRequiredParameter(
  73             name     = "MutationDifferentialWeight_F",
  74             default  = (0.5, 1),
  75             typecast = tuple,
  76             message  = "Poids différentiel de mutation, constant ou aléatoire dans l'intervalle, noté F",
  77             minval   = 0.,
  78             maxval   = 2.,
  79             )
  80         self.defineRequiredParameter(
  81             name     = "CrossOverProbability_CR",
  82             default  = 0.7,
  83             typecast = float,
  84             message  = "Probabilité de recombinaison ou de croisement, notée CR",
  85             minval   = 0.,
  86             maxval   = 1.,
  87             )
  88         self.defineRequiredParameter(
  89             name     = "QualityCriterion",
  90             default  = "AugmentedWeightedLeastSquares",
  91             typecast = str,
  92             message  = "Critère de qualité utilisé",
  93             listval  = ["AugmentedWeightedLeastSquares","AWLS","DA",
  94                         "WeightedLeastSquares","WLS",
  95                         "LeastSquares","LS","L2",
  96                         "AbsoluteValue","L1",
  97                         "MaximumError","ME"],
  98             )
  99         self.defineRequiredParameter(
 100             name     = "StoreInternalVariables",
 101             default  = False,
 102             typecast = bool,
 103             message  = "Stockage des variables internes ou intermédiaires du calcul",
 104             )
 105         self.defineRequiredParameter(
 106             name     = "StoreSupplementaryCalculations",
 107             default  = [],
 108             typecast = tuple,
 109             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
 110             listval  = [
 111                 "Analysis",
 112                 "BMA",
 113                 "CostFunctionJ",
 114                 "CostFunctionJb",
 115                 "CostFunctionJo",
 116                 "CostFunctionJAtCurrentOptimum",
 117                 "CostFunctionJbAtCurrentOptimum",
 118                 "CostFunctionJoAtCurrentOptimum",
 119                 "CurrentOptimum",
 120                 "CurrentState",
 121                 "IndexOfOptimum",
 122                 "Innovation",
 123                 "InnovationAtCurrentState",
 124                 "OMA",
 125                 "OMB",
 126                 "SimulatedObservationAtBackground",
 127                 "SimulatedObservationAtCurrentOptimum",
 128                 "SimulatedObservationAtCurrentState",
 129                 "SimulatedObservationAtOptimum",
 130                 ]
 131             )
 132         self.defineRequiredParameter(
 133             name     = "SetSeed",
 134             typecast = numpy.random.seed,
 135             message  = "Graine fixée pour le générateur aléatoire",
 136             )
 137         self.defineRequiredParameter( # Pas de type
 138             name     = "Bounds",
 139             message  = "Liste des valeurs de bornes",
 140             )
 141         self.requireInputArguments(
 142             mandatory= ("Xb", "Y", "HO", "R", "B" ),
 143             )
 144
 145     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
 146         self._pre_run(Parameters, Xb, Y, R, B, Q)
 147         #
 148         len_X = numpy.asarray(Xb).size
 149         popsize = round(self._parameters["PopulationSize"]/len_X)
 150         maxiter = min(self._parameters["MaximumNumberOfSteps"],round(self._parameters["MaximumNumberOfFunctionEvaluations"]/(popsize*len_X) - 1))
 151         logging.debug("%s Nombre maximal de générations = %i, taille de la population à chaque génération = %i"%(self._name, maxiter, popsize*len_X))
 152         #
 153         # Opérateurs
 154         # ----------
 155         Hm = HO["Direct"].appliedTo
 156         #
 157         # Précalcul des inversions de B et R
 158         # ----------------------------------
 159         BI = B.getI()
 160         RI = R.getI()
 161         #
 162         # Définition de la fonction-coût
 163         # ------------------------------
 164         def CostFunction(x, QualityMeasure="AugmentedWeightedLeastSquares"):
 165             _X  = numpy.asmatrix(numpy.ravel( x )).T
 166             self.StoredVariables["CurrentState"].store( _X )
 167             _HX = Hm( _X )
 168             _HX = numpy.asmatrix(numpy.ravel( _HX )).T
 169             _Innovation = Y - _HX
 170             if self._toStore("SimulatedObservationAtCurrentState") or \
 171                 self._toStore("SimulatedObservationAtCurrentOptimum"):
 172                 self.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
 173             if self._toStore("InnovationAtCurrentState"):
 174                 self.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
 175             #
 176             if QualityMeasure in ["AugmentedWeightedLeastSquares","AWLS","DA"]:
 177                 if BI is None or RI is None:
 178                     raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
 179                 Jb  = 0.5 * (_X - Xb).T * BI * (_X - Xb)
 180                 Jo  = 0.5 * (_Innovation).T * RI * (_Innovation)
 181             elif QualityMeasure in ["WeightedLeastSquares","WLS"]:
 182                 if RI is None:
 183                     raise ValueError("Observation error covariance matrix has to be properly defined!")
 184                 Jb  = 0.
 185                 Jo  = 0.5 * (_Innovation).T * RI * (_Innovation)
 186             elif QualityMeasure in ["LeastSquares","LS","L2"]:
 187                 Jb  = 0.
 188                 Jo  = 0.5 * (_Innovation).T * (_Innovation)
 189             elif QualityMeasure in ["AbsoluteValue","L1"]:
 190                 Jb  = 0.
 191                 Jo  = numpy.sum( numpy.abs(_Innovation) )
 192             elif QualityMeasure in ["MaximumError","ME"]:
 193                 Jb  = 0.
 194                 Jo  = numpy.max( numpy.abs(_Innovation) )
 195             #
 196             J   = float( Jb ) + float( Jo )
 197             #
 198             self.StoredVariables["CostFunctionJb"].store( Jb )
 199             self.StoredVariables["CostFunctionJo"].store( Jo )
 200             self.StoredVariables["CostFunctionJ" ].store( J )
 201             if self._toStore("IndexOfOptimum") or \
 202                 self._toStore("CurrentOptimum") or \
 203                 self._toStore("CostFunctionJAtCurrentOptimum") or \
 204                 self._toStore("CostFunctionJbAtCurrentOptimum") or \
 205                 self._toStore("CostFunctionJoAtCurrentOptimum") or \
 206                 self._toStore("SimulatedObservationAtCurrentOptimum"):
 207                 IndexMin = numpy.argmin( self.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
 208             if self._toStore("IndexOfOptimum"):
 209                 self.StoredVariables["IndexOfOptimum"].store( IndexMin )
 210             if self._toStore("CurrentOptimum"):
 211                 self.StoredVariables["CurrentOptimum"].store( self.StoredVariables["CurrentState"][IndexMin] )
 212             if self._toStore("SimulatedObservationAtCurrentOptimum"):
 213                 self.StoredVariables["SimulatedObservationAtCurrentOptimum"].store( self.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin] )
 214             if self._toStore("CostFunctionJAtCurrentOptimum"):
 215                 self.StoredVariables["CostFunctionJAtCurrentOptimum" ].store( self.StoredVariables["CostFunctionJ" ][IndexMin] )
 216             if self._toStore("CostFunctionJbAtCurrentOptimum"):
 217                 self.StoredVariables["CostFunctionJbAtCurrentOptimum"].store( self.StoredVariables["CostFunctionJb"][IndexMin] )
 218             if self._toStore("CostFunctionJoAtCurrentOptimum"):
 219                 self.StoredVariables["CostFunctionJoAtCurrentOptimum"].store( self.StoredVariables["CostFunctionJo"][IndexMin] )
 220             return J
 221         #
 222         # Point de démarrage de l'optimisation : Xini = Xb
 223         # ------------------------------------
 224         Xini = numpy.ravel(Xb)
 225         #
 226         # Minimisation de la fonctionnelle
 227         # --------------------------------
 228         nbPreviousSteps = self.StoredVariables["CostFunctionJ"].stepnumber()
 229         #
 230         optResults = scipy.optimize.differential_evolution(
 231             CostFunction,
 232             self._parameters["Bounds"],
 233             strategy      = str(self._parameters["Minimizer"]).lower(),
 234             maxiter       = maxiter,
 235             popsize       = popsize,
 236             mutation      = self._parameters["MutationDifferentialWeight_F"],
 237             recombination = self._parameters["CrossOverProbability_CR"],
 238             disp          = self._parameters["optdisp"],
 239             )
 240         #
 241         IndexMin = numpy.argmin( self.StoredVariables["CostFunctionJ"][nbPreviousSteps:] ) + nbPreviousSteps
 242         MinJ     = self.StoredVariables["CostFunctionJ"][IndexMin]
 243         Minimum  = self.StoredVariables["CurrentState"][IndexMin]
 244         #
 245         # Obtention de l'analyse
 246         # ----------------------
 247         Xa = numpy.ravel( Minimum )
 248         #
 249         self.StoredVariables["Analysis"].store( Xa )
 250         #
 251         # Calculs et/ou stockages supplémentaires
 252         # ---------------------------------------
 253         if self._toStore("OMA") or self._toStore("SimulatedObservationAtOptimum"):
 254             if self._toStore("SimulatedObservationAtCurrentState"):
 255                 HXa = self.StoredVariables["SimulatedObservationAtCurrentState"][IndexMin]
 256             elif self._toStore("SimulatedObservationAtCurrentOptimum"):
 257                 HXa = self.StoredVariables["SimulatedObservationAtCurrentOptimum"][-1]
 258             else:
 259                 HXa = Hm(Xa)
 260         if self._toStore("Innovation") or \
 261            self._toStore("OMB"):
 262             d  = Y - HXb
 263         if self._toStore("Innovation"):
 264             self.StoredVariables["Innovation"].store( numpy.ravel(d) )
 265         if self._toStore("OMB"):
 266             self.StoredVariables["OMB"].store( numpy.ravel(d) )
 267         if self._toStore("BMA"):
 268             self.StoredVariables["BMA"].store( numpy.ravel(Xb) - numpy.ravel(Xa) )
 269         if self._toStore("OMA"):
 270             self.StoredVariables["OMA"].store( numpy.ravel(Y) - numpy.ravel(HXa) )
 271         if self._toStore("SimulatedObservationAtBackground"):
 272             self.StoredVariables["SimulatedObservationAtBackground"].store( numpy.ravel(Hm(Xb)) )
 273         if self._toStore("SimulatedObservationAtOptimum"):
 274             self.StoredVariables["SimulatedObservationAtOptimum"].store( numpy.ravel(HXa) )
 275         #
 276         self._post_run()
 277         return 0
 278
 279 # ==============================================================================
 280 if __name__ == "__main__":
 281     print('\n AUTODIAGNOSTIC\n')