src/daComposant/daAlgorithms/ParticleSwarmOptimization.py

   1 #-*-coding:iso-8859-1-*-
   2 #
   3 #  Copyright (C) 2008-2012 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
  22 import logging
  23 from daCore import BasicObjects, PlatformInfo
  24 m = PlatformInfo.SystemUsage()
  25
  26 import numpy
  27 import copy
  28
  29 # ==============================================================================
  30 class ElementaryAlgorithm(BasicObjects.Algorithm):
  31     def __init__(self):
  32         BasicObjects.Algorithm.__init__(self, "PARTICLESWARMOPTIMIZATION")
  33         self.defineRequiredParameter(
  34             name     = "MaximumNumberOfSteps",
  35             default  = 50,
  36             typecast = int,
  37             message  = "Nombre maximal de pas d'optimisation",
  38             minval   = 1,
  39             )
  40         self.defineRequiredParameter(
  41             name     = "SetSeed",
  42             typecast = numpy.random.seed,
  43             message  = "Graine fixée pour le générateur aléatoire",
  44             )
  45         self.defineRequiredParameter(
  46             name     = "NumberOfInsects",
  47             default  = 100,
  48             typecast = int,
  49             message  = "Nombre d'insectes dans l'essaim",
  50             minval   = -1,
  51             )
  52         self.defineRequiredParameter(
  53             name     = "SwarmVelocity",
  54             default  = 1.,
  55             typecast = float,
  56             message  = "Vitesse de groupe imposée par l'essaim",
  57             minval   = 0.,
  58             )
  59         self.defineRequiredParameter(
  60             name     = "GroupRecallRate",
  61             default  = 0.5,
  62             typecast = float,
  63             message  = "Taux de rappel au meilleur insecte du groupe (entre 0 et 1)",
  64             minval   = 0.,
  65             maxval   = 1.,
  66             )
  67         self.defineRequiredParameter(
  68             name     = "QualityCriterion",
  69             default  = "AugmentedPonderatedLeastSquares",
  70             typecast = str,
  71             message  = "Critère de qualité utilisé",
  72             listval  = ["AugmentedPonderatedLeastSquares","APLS","DA",
  73                         "PonderatedLeastSquares","PLS",
  74                         "LeastSquares","LS","L2",
  75                         "AbsoluteValue","L1",
  76                         "MaximumError","ME"],
  77             )
  78         self.defineRequiredParameter(
  79             name     = "StoreInternalVariables",
  80             default  = False,
  81             typecast = bool,
  82             message  = "Stockage des variables internes ou intermédiaires du calcul",
  83             )
  84
  85     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
  86         """
  87         Calcul de l'estimateur
  88         """
  89         logging.debug("%s Lancement"%self._name)
  90         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("Mo")))
  91         #
  92         # Paramètres de pilotage
  93         # ----------------------
  94         self.setParameters(Parameters)
  95         #
  96         if self._parameters.has_key("BoxBounds") and (type(self._parameters["BoxBounds"]) is type([]) or type(self._parameters["BoxBounds"]) is type(())) and (len(self._parameters["BoxBounds"]) > 0):
  97             BoxBounds = self._parameters["BoxBounds"]
  98             logging.debug("%s Prise en compte des bornes d'incréments de paramètres effectuee"%(self._name,))
  99         else:
 100             raise ValueError("Particle Swarm Optimization requires bounds on all variables to be given.")
 101         BoxBounds   = numpy.array(BoxBounds)
 102         if numpy.isnan(BoxBounds).any():
 103             raise ValueError("Particle Swarm Optimization requires bounds on all variables increments to be truly given, \"None\" is not allowed. The actual increments bounds are:\n%s"%BoxBounds)
 104         #
 105         Phig = float( self._parameters["GroupRecallRate"] )
 106         Phip = 1. - Phig
 107         logging.debug("%s Taux de rappel au meilleur insecte du groupe (entre 0 et 1) = %s et à la meilleure position précédente (son complémentaire à 1) = %s"%(self._name, str(Phig), str(Phip)))
 108         #
 109         # Opérateur d'observation
 110         # -----------------------
 111         Hm = H["Direct"].appliedTo
 112         #
 113         # Précalcul des inversions de B et R
 114         # ----------------------------------
 115         if B is not None:
 116             BI = B.I
 117         elif self._parameters["B_scalar"] is not None:
 118             BI = 1.0 / self._parameters["B_scalar"]
 119         else:
 120             BI = None
 121         #
 122         if R is not None:
 123             RI = R.I
 124         elif self._parameters["R_scalar"] is not None:
 125             RI = 1.0 / self._parameters["R_scalar"]
 126         else:
 127             RI = None
 128         #
 129         # Définition de la fonction-coût
 130         # ------------------------------
 131         def CostFunction(x, QualityMeasure="AugmentedPonderatedLeastSquares"):
 132             _X  = numpy.asmatrix(x).flatten().T
 133             logging.debug("%s CostFunction X  = %s"%(self._name, _X.A1))
 134             _HX = Hm( _X )
 135             _HX = numpy.asmatrix(_HX).flatten().T
 136             #
 137             if QualityMeasure in ["AugmentedPonderatedLeastSquares","APLS","DA"]:
 138                 if BI is None or RI is None:
 139                     raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
 140                 Jb  = 0.5 * (_X - Xb).T * BI * (_X - Xb)
 141                 Jo  = 0.5 * (Y - _HX).T * RI * (Y - _HX)
 142                 J   = float( Jb ) + float( Jo )
 143             elif QualityMeasure in ["PonderatedLeastSquares","PLS"]:
 144                 if RI is None:
 145                     raise ValueError("Observation error covariance matrix has to be properly defined!")
 146                 Jb  = 0.
 147                 Jo  = 0.5 * (Y - _HX).T * RI * (Y - _HX)
 148                 J   = float( Jb ) + float( Jo )
 149             elif QualityMeasure in ["LeastSquares","LS","L2"]:
 150                 Jb  = 0.
 151                 Jo  = 0.5 * (Y - _HX).T * (Y - _HX)
 152                 J   = float( Jb ) + float( Jo )
 153             elif QualityMeasure in ["AbsoluteValue","L1"]:
 154                 Jb  = 0.
 155                 Jo  = numpy.sum( numpy.abs(Y - _HX) )
 156                 J   = float( Jb ) + float( Jo )
 157             elif QualityMeasure in ["MaximumError","ME"]:
 158                 Jb  = 0.
 159                 Jo  = numpy.max( numpy.abs(Y - _HX) )
 160                 J   = float( Jb ) + float( Jo )
 161             #
 162             logging.debug("%s CostFunction Jb = %s"%(self._name, Jb))
 163             logging.debug("%s CostFunction Jo = %s"%(self._name, Jo))
 164             logging.debug("%s CostFunction J  = %s"%(self._name, J))
 165             return J
 166         #
 167         # Point de démarrage de l'optimisation : Xini = Xb
 168         # ------------------------------------
 169         if type(Xb) is type(numpy.matrix([])):
 170             Xini = Xb.A1.tolist()
 171         elif Xb is not None:
 172             Xini = list(Xb)
 173         else:
 174             Xini = numpy.zeros(len(BoxBounds[:,0]))
 175         logging.debug("%s Point de démarrage Xini = %s"%(self._name, Xini))
 176         #
 177         # Initialisation des bornes
 178         # -------------------------
 179         SpaceUp  = BoxBounds[:,1] + Xini
 180         Spacelow = BoxBounds[:,0] + Xini
 181         nbparam  = len(SpaceUp)
 182         #
 183         # Initialisation de l'essaim
 184         # --------------------------
 185         LimitVelocity = numpy.abs(SpaceUp-Spacelow)
 186         #
 187         PosInsect = []
 188         VelocityInsect = []
 189         for i in range(nbparam) :
 190             PosInsect.append(numpy.random.uniform(low=Spacelow[i], high=SpaceUp[i], size=self._parameters["NumberOfInsects"]))
 191             VelocityInsect.append(numpy.random.uniform(low=-LimitVelocity[i], high=LimitVelocity[i], size=self._parameters["NumberOfInsects"]))
 192         VelocityInsect = numpy.matrix(VelocityInsect)
 193         PosInsect = numpy.matrix(PosInsect)
 194         #
 195         BestPosInsect = numpy.array(PosInsect)
 196         qBestPosInsect = []
 197         Best = copy.copy(Spacelow)
 198         qBest = CostFunction(Best,self._parameters["QualityCriterion"])
 199         #
 200         for i in range(self._parameters["NumberOfInsects"]):
 201             insect  = numpy.array(PosInsect[:,i].A1)
 202             quality = CostFunction(insect,self._parameters["QualityCriterion"])
 203             qBestPosInsect.append(quality)
 204             if quality < qBest:
 205                 Best  = insect
 206                 qBest = quality
 207         #
 208         # Minimisation de la fonctionnelle
 209         # --------------------------------
 210         for n in range(self._parameters["MaximumNumberOfSteps"]):
 211             for i in range(self._parameters["NumberOfInsects"]) :
 212                 insect  = PosInsect[:,i]
 213                 rp = numpy.random.uniform(size=nbparam)
 214                 rg = numpy.random.uniform(size=nbparam)
 215                 for j in range(nbparam) :
 216                     VelocityInsect[j,i] = self._parameters["SwarmVelocity"]*VelocityInsect[j,i] +  Phip*rp[j]*(BestPosInsect[j,i]-PosInsect[j,i]) +  Phig*rg[j]*(Best[j]-PosInsect[j,i])
 217                     PosInsect[j,i] = PosInsect[j,i]+VelocityInsect[j,i]
 218                 quality = CostFunction(insect,self._parameters["QualityCriterion"])
 219                 if quality < qBestPosInsect[i]:
 220                     BestPosInsect[:,i] = numpy.asmatrix(insect).flatten().A1
 221                     if quality < qBest :
 222                         Best  = numpy.asmatrix(insect).flatten().A1
 223                         qBest = quality
 224             logging.debug("%s Iteration %i : qBest = %.5f, Best = %s"%(self._name, n+1,qBest,Best))
 225             #
 226             if self._parameters["StoreInternalVariables"]:
 227                 self.StoredVariables["CurrentState"].store( Best )
 228             self.StoredVariables["CostFunctionJb"].store( 0. )
 229             self.StoredVariables["CostFunctionJo"].store( 0. )
 230             self.StoredVariables["CostFunctionJ" ].store( qBest )
 231         #
 232         logging.debug("%s %s Step of min cost  = %s"%(self._name, self._parameters["QualityCriterion"], self._parameters["MaximumNumberOfSteps"]))
 233         logging.debug("%s %s Minimum cost      = %s"%(self._name, self._parameters["QualityCriterion"], qBest))
 234         logging.debug("%s %s Minimum state     = %s"%(self._name, self._parameters["QualityCriterion"], Best))
 235         logging.debug("%s %s Nb of F           = %s"%(self._name, self._parameters["QualityCriterion"], (self._parameters["MaximumNumberOfSteps"]+1)*self._parameters["NumberOfInsects"]+1))
 236         logging.debug("%s %s RetCode           = %s"%(self._name, self._parameters["QualityCriterion"], 0))
 237         #
 238         # Obtention de l'analyse
 239         # ----------------------
 240         Xa = numpy.asmatrix(Best).flatten().T
 241         logging.debug("%s Analyse Xa = %s"%(self._name, Xa))
 242         #
 243         self.StoredVariables["Analysis"].store( Xa.A1 )
 244         #
 245         logging.debug("%s Taille mémoire utilisée de %.1f Mo"%(self._name, m.getUsedMemory("MB")))
 246         logging.debug("%s Terminé"%self._name)
 247         #
 248         return 0
 249
 250 # ==============================================================================
 251 if __name__ == "__main__":
 252     print '\n AUTODIAGNOSTIC \n'