1 #-*-coding:iso-8859-1-*-
3 # Copyright (C) 2008-2017 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
24 Définit les outils généraux élémentaires.
26 Ce module est destiné à être appelée par AssimilationStudy.
28 __author__ = "Jean-Philippe ARGAUD"
31 import logging, copy, os
33 from daCore import Persistence
34 from daCore import PlatformInfo
36 # ==============================================================================
37 class CacheManager(object):
39 Classe générale de gestion d'un cache de calculs
42 toleranceInRedundancy = 1.e-18,
43 lenghtOfRedundancy = -1,
46 Les caractéristiques de tolérance peuvent être modifées à la création.
48 self.__tolerBP = float(toleranceInRedundancy)
49 self.__lenghtOR = int(lenghtOfRedundancy)
50 self.__initlnOR = self.__lenghtOR
55 self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
56 # logging.debug("CM Tolerance de determination des doublons : %.2e", self.__tolerBP)
58 def wasCalculatedIn(self, xValue ): #, info="" ):
59 "Vérifie l'existence d'un calcul correspondant à la valeur"
62 for i in range(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
63 if xValue.size != self.__listOPCV[i][0].size:
64 # logging.debug("CM Différence de la taille %s de X et de celle %s du point %i déjà calculé", xValue.shape,i,self.__listOPCP[i].shape)
66 if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
68 __HxV = self.__listOPCV[i][1]
69 # logging.debug("CM Cas%s déja calculé, portant le numéro %i", info, i)
73 def storeValueInX(self, xValue, HxValue ):
74 "Stocke un calcul correspondant à la valeur"
75 if self.__lenghtOR < 0:
76 self.__lenghtOR = 2 * xValue.size + 2
77 self.__initlnOR = self.__lenghtOR
78 while len(self.__listOPCV) > self.__lenghtOR:
79 # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier", self.__lenghtOR)
80 self.__listOPCV.pop(0)
81 self.__listOPCV.append( (
82 copy.copy(numpy.ravel(xValue)),
84 numpy.linalg.norm(xValue),
89 self.__initlnOR = self.__lenghtOR
94 self.__lenghtOR = self.__initlnOR
96 # ==============================================================================
97 class Operator(object):
99 Classe générale d'interface de type opérateur
106 def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
108 On construit un objet de ce type en fournissant à l'aide de l'un des
109 deux mots-clé, soit une fonction python, soit une matrice.
111 - fromMethod : argument de type fonction Python
112 - fromMatrix : argument adapté au constructeur numpy.matrix
113 - avoidingRedundancy : évite ou pas les calculs redondants
115 self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
116 self.__AvoidRC = bool( avoidingRedundancy )
117 if fromMethod is not None:
118 self.__Method = fromMethod
120 self.__Type = "Method"
121 elif fromMatrix is not None:
123 self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
124 self.__Type = "Matrix"
130 def disableAvoidingRedundancy(self):
132 Operator.CM.disable()
134 def enableAvoidingRedundancy(self):
139 Operator.CM.disable()
145 def appliedTo(self, xValue, HValue = None):
147 Permet de restituer le résultat de l'application de l'opérateur à un
148 argument xValue. Cette méthode se contente d'appliquer, son argument
149 devant a priori être du bon type.
151 - xValue : argument adapté pour appliquer l'opérateur
153 if HValue is not None:
154 HxValue = numpy.asmatrix( numpy.ravel( HValue ) ).T
156 Operator.CM.storeValueInX(xValue,HxValue)
159 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
161 __alreadyCalculated = False
163 if __alreadyCalculated:
164 self.__addOneCacheCall()
167 if self.__Matrix is not None:
168 self.__addOneMatrixCall()
169 HxValue = self.__Matrix * xValue
171 self.__addOneMethodCall()
172 HxValue = self.__Method( xValue )
174 Operator.CM.storeValueInX(xValue,HxValue)
178 def appliedControledFormTo(self, paire ):
180 Permet de restituer le résultat de l'application de l'opérateur à une
181 paire (xValue, uValue). Cette méthode se contente d'appliquer, son
182 argument devant a priori être du bon type. Si la uValue est None,
183 on suppose que l'opérateur ne s'applique qu'à xValue.
185 - xValue : argument X adapté pour appliquer l'opérateur
186 - uValue : argument U adapté pour appliquer l'opérateur
188 assert len(paire) == 2, "Incorrect number of arguments"
189 xValue, uValue = paire
190 if self.__Matrix is not None:
191 self.__addOneMatrixCall()
192 return self.__Matrix * xValue
193 elif uValue is not None:
194 self.__addOneMethodCall()
195 return self.__Method( (xValue, uValue) )
197 self.__addOneMethodCall()
198 return self.__Method( xValue )
200 def appliedInXTo(self, paire ):
202 Permet de restituer le résultat de l'application de l'opérateur à un
203 argument xValue, sachant que l'opérateur est valable en xNominal.
204 Cette méthode se contente d'appliquer, son argument devant a priori
205 être du bon type. Si l'opérateur est linéaire car c'est une matrice,
206 alors il est valable en tout point nominal et il n'est pas nécessaire
208 Arguments : une liste contenant
209 - xNominal : argument permettant de donner le point où l'opérateur
210 est construit pour etre ensuite appliqué
211 - xValue : argument adapté pour appliquer l'opérateur
213 assert len(paire) == 2, "Incorrect number of arguments"
214 xNominal, xValue = paire
215 if self.__Matrix is not None:
216 self.__addOneMatrixCall()
217 return self.__Matrix * xValue
219 self.__addOneMethodCall()
220 return self.__Method( (xNominal, xValue) )
222 def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
224 Permet de renvoyer l'opérateur sous la forme d'une matrice
226 if self.__Matrix is not None:
227 self.__addOneMatrixCall()
229 elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
230 self.__addOneMethodCall()
231 return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
233 raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
237 Renvoie la taille sous forme numpy si l'opérateur est disponible sous
238 la forme d'une matrice
240 if self.__Matrix is not None:
241 return self.__Matrix.shape
243 raise ValueError("Matrix form of the operator is not available, nor the shape")
245 def nbcalls(self, which=None):
247 Renvoie les nombres d'évaluations de l'opérateur
250 self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
251 self.__NbCallsAsMatrix,
252 self.__NbCallsAsMethod,
253 self.__NbCallsOfCached,
254 Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
255 Operator.NbCallsAsMatrix,
256 Operator.NbCallsAsMethod,
257 Operator.NbCallsOfCached,
259 if which is None: return __nbcalls
260 else: return __nbcalls[which]
262 def __addOneMatrixCall(self):
263 "Comptabilise un appel"
264 self.__NbCallsAsMatrix += 1 # Decompte local
265 Operator.NbCallsAsMatrix += 1 # Decompte global
267 def __addOneMethodCall(self):
268 "Comptabilise un appel"
269 self.__NbCallsAsMethod += 1 # Decompte local
270 Operator.NbCallsAsMethod += 1 # Decompte global
272 def __addOneCacheCall(self):
273 "Comptabilise un appel"
274 self.__NbCallsOfCached += 1 # Decompte local
275 Operator.NbCallsOfCached += 1 # Decompte global
277 # ==============================================================================
278 class Algorithm(object):
280 Classe générale d'interface de type algorithme
282 Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
283 d'assimilation, en fournissant un container (dictionnaire) de variables
284 persistantes initialisées, et des méthodes d'accès à ces variables stockées.
286 Une classe élémentaire d'algorithme doit implémenter la méthode "run".
288 def __init__(self, name):
290 L'initialisation présente permet de fabriquer des variables de stockage
291 disponibles de manière générique dans les algorithmes élémentaires. Ces
292 variables de stockage sont ensuite conservées dans un dictionnaire
293 interne à l'objet, mais auquel on accède par la méthode "get".
295 Les variables prévues sont :
296 - CostFunctionJ : fonction-cout globale, somme des deux parties suivantes
297 - CostFunctionJb : partie ébauche ou background de la fonction-cout
298 - CostFunctionJo : partie observations de la fonction-cout
299 - GradientOfCostFunctionJ : gradient de la fonction-cout globale
300 - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-cout
301 - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-cout
302 - CurrentState : état courant lors d'itérations
303 - Analysis : l'analyse Xa
304 - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
305 - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
306 - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
307 - Innovation : l'innovation : d = Y - H(X)
308 - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
309 - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
310 - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
311 - MahalanobisConsistency : indicateur de consistance des covariances
312 - OMA : Observation moins Analysis : Y - Xa
313 - OMB : Observation moins Background : Y - Xb
314 - AMB : Analysis moins Background : Xa - Xb
315 - APosterioriCovariance : matrice A
316 - APosterioriVariances : variances de la matrice A
317 - APosterioriStandardDeviations : écart-types de la matrice A
318 - APosterioriCorrelations : correlations de la matrice A
319 - Residu : dans le cas des algorithmes de vérification
320 On peut rajouter des variables à stocker dans l'initialisation de
321 l'algorithme élémentaire qui va hériter de cette classe
323 logging.debug("%s Initialisation", str(name))
324 self._m = PlatformInfo.SystemUsage()
326 self._name = str( name )
327 self._parameters = {"StoreSupplementaryCalculations":[]}
328 self.__required_parameters = {}
329 self.StoredVariables = {}
331 self.StoredVariables["CostFunctionJ"] = Persistence.OneScalar(name = "CostFunctionJ")
332 self.StoredVariables["CostFunctionJb"] = Persistence.OneScalar(name = "CostFunctionJb")
333 self.StoredVariables["CostFunctionJo"] = Persistence.OneScalar(name = "CostFunctionJo")
334 self.StoredVariables["CostFunctionJAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
335 self.StoredVariables["CostFunctionJbAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
336 self.StoredVariables["CostFunctionJoAtCurrentOptimum"] = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
337 self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
338 self.StoredVariables["GradientOfCostFunctionJb"] = Persistence.OneVector(name = "GradientOfCostFunctionJb")
339 self.StoredVariables["GradientOfCostFunctionJo"] = Persistence.OneVector(name = "GradientOfCostFunctionJo")
340 self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
341 self.StoredVariables["Analysis"] = Persistence.OneVector(name = "Analysis")
342 self.StoredVariables["IndexOfOptimum"] = Persistence.OneIndex(name = "IndexOfOptimum")
343 self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
344 self.StoredVariables["SimulatedObservationAtBackground"] = Persistence.OneVector(name = "SimulatedObservationAtBackground")
345 self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
346 self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
347 self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
348 self.StoredVariables["Innovation"] = Persistence.OneVector(name = "Innovation")
349 self.StoredVariables["InnovationAtCurrentState"] = Persistence.OneVector(name = "InnovationAtCurrentState")
350 self.StoredVariables["SigmaObs2"] = Persistence.OneScalar(name = "SigmaObs2")
351 self.StoredVariables["SigmaBck2"] = Persistence.OneScalar(name = "SigmaBck2")
352 self.StoredVariables["MahalanobisConsistency"] = Persistence.OneScalar(name = "MahalanobisConsistency")
353 self.StoredVariables["OMA"] = Persistence.OneVector(name = "OMA")
354 self.StoredVariables["OMB"] = Persistence.OneVector(name = "OMB")
355 self.StoredVariables["BMA"] = Persistence.OneVector(name = "BMA")
356 self.StoredVariables["APosterioriCovariance"] = Persistence.OneMatrix(name = "APosterioriCovariance")
357 self.StoredVariables["APosterioriVariances"] = Persistence.OneVector(name = "APosterioriVariances")
358 self.StoredVariables["APosterioriStandardDeviations"] = Persistence.OneVector(name = "APosterioriStandardDeviations")
359 self.StoredVariables["APosterioriCorrelations"] = Persistence.OneMatrix(name = "APosterioriCorrelations")
360 self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
361 self.StoredVariables["Residu"] = Persistence.OneScalar(name = "Residu")
363 def _pre_run(self, Parameters ):
365 logging.debug("%s Lancement", self._name)
366 logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
368 # Mise a jour de self._parameters avec Parameters
369 self.__setParameters(Parameters)
371 # Corrections et complements
372 if "Bounds" in self._parameters and (type(self._parameters["Bounds"]) is type([]) or type(self._parameters["Bounds"]) is type(())) and (len(self._parameters["Bounds"]) > 0):
373 logging.debug("%s Prise en compte des bornes effectuee"%(self._name,))
375 self._parameters["Bounds"] = None
377 if logging.getLogger().level < logging.WARNING:
378 self._parameters["optiprint"], self._parameters["optdisp"] = 1, 1
379 if PlatformInfo.has_scipy:
380 import scipy.optimize
381 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_ALL
383 self._parameters["optmessages"] = 15
385 self._parameters["optiprint"], self._parameters["optdisp"] = -1, 0
386 if PlatformInfo.has_scipy:
387 import scipy.optimize
388 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_NONE
390 self._parameters["optmessages"] = 15
394 def _post_run(self,_oH=None):
396 if ("StoreSupplementaryCalculations" in self._parameters) and \
397 "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
398 for _A in self.StoredVariables["APosterioriCovariance"]:
399 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
400 self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
401 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
402 self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
403 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
404 _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
405 _C = numpy.dot(_EI, numpy.dot(_A, _EI))
406 self.StoredVariables["APosterioriCorrelations"].store( _C )
408 logging.debug("%s Nombre d'évaluation(s) de l'opérateur d'observation direct/tangent/adjoint.: %i/%i/%i", self._name, _oH["Direct"].nbcalls(0),_oH["Tangent"].nbcalls(0),_oH["Adjoint"].nbcalls(0))
409 logging.debug("%s Nombre d'appels au cache d'opérateur d'observation direct/tangent/adjoint..: %i/%i/%i", self._name, _oH["Direct"].nbcalls(3),_oH["Tangent"].nbcalls(3),_oH["Adjoint"].nbcalls(3))
410 logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
411 logging.debug("%s Terminé", self._name)
414 def get(self, key=None):
416 Renvoie l'une des variables stockées identifiée par la clé, ou le
417 dictionnaire de l'ensemble des variables disponibles en l'absence de
418 clé. Ce sont directement les variables sous forme objet qui sont
419 renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
420 des classes de persistance.
423 return self.StoredVariables[key]
425 return self.StoredVariables
427 def __contains__(self, key=None):
428 "D.__contains__(k) -> True if D has a key k, else False"
429 return key in self.StoredVariables
432 "D.keys() -> list of D's keys"
433 return self.StoredVariables.keys()
435 def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
437 Doit implémenter l'opération élémentaire de calcul d'assimilation sous
438 sa forme mathématique la plus naturelle possible.
440 raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
442 def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
444 Permet de définir dans l'algorithme des paramètres requis et leurs
445 caractéristiques par défaut.
448 raise ValueError("A name is mandatory to define a required parameter.")
450 self.__required_parameters[name] = {
452 "typecast" : typecast,
458 logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
460 def getRequiredParameters(self, noDetails=True):
462 Renvoie la liste des noms de paramètres requis ou directement le
463 dictionnaire des paramètres requis.
466 ks = list(self.__required_parameters.keys())
470 return self.__required_parameters
472 def setParameterValue(self, name=None, value=None):
474 Renvoie la valeur d'un paramètre requis de manière contrôlée
476 default = self.__required_parameters[name]["default"]
477 typecast = self.__required_parameters[name]["typecast"]
478 minval = self.__required_parameters[name]["minval"]
479 maxval = self.__required_parameters[name]["maxval"]
480 listval = self.__required_parameters[name]["listval"]
482 if value is None and default is None:
484 elif value is None and default is not None:
485 if typecast is None: __val = default
486 else: __val = typecast( default )
488 if typecast is None: __val = value
489 else: __val = typecast( value )
491 if minval is not None and (numpy.array(__val, float) < minval).any():
492 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
493 if maxval is not None and (numpy.array(__val, float) > maxval).any():
494 raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
495 if listval is not None:
496 if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
499 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
500 elif __val not in listval:
501 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
504 def __setParameters(self, fromDico={}):
506 Permet de stocker les paramètres reçus dans le dictionnaire interne.
508 self._parameters.update( fromDico )
509 for k in self.__required_parameters.keys():
510 if k in fromDico.keys():
511 self._parameters[k] = self.setParameterValue(k,fromDico[k])
513 self._parameters[k] = self.setParameterValue(k)
514 logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
516 # ==============================================================================
517 class Diagnostic(object):
519 Classe générale d'interface de type diagnostic
521 Ce template s'utilise de la manière suivante : il sert de classe "patron" en
522 même temps que l'une des classes de persistance, comme "OneScalar" par
525 Une classe élémentaire de diagnostic doit implémenter ses deux méthodes, la
526 méthode "_formula" pour écrire explicitement et proprement la formule pour
527 l'écriture mathématique du calcul du diagnostic (méthode interne non
528 publique), et "calculate" pour activer la précédente tout en ayant vérifié
529 et préparé les données, et pour stocker les résultats à chaque pas (méthode
530 externe d'activation).
532 def __init__(self, name = "", parameters = {}):
534 self.name = str(name)
535 self.parameters = dict( parameters )
537 def _formula(self, *args):
539 Doit implémenter l'opération élémentaire de diagnostic sous sa forme
540 mathématique la plus naturelle possible.
542 raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
544 def calculate(self, *args):
546 Active la formule de calcul avec les arguments correctement rangés
548 raise NotImplementedError("Diagnostic activation method has not been implemented!")
550 # ==============================================================================
551 class ParameterDictionary(object):
553 Classe générale d'interface de type dictionnaire de paramètres
556 name = "GenericParamDict",
563 self.__name = str(name)
567 if asScript is not None:
568 __Algo = ImportFromScript(asScript).getvalue( "Algorithm" )
569 __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
574 if __Algo is not None:
575 self.__A = str(__Algo)
576 if __Dict is not None:
577 self.__D = dict(__Dict)
580 "x.__repr__() <==> repr(x)"
581 return repr(self.__A)+"\n"+repr(self.__D)
584 "x.__str__() <==> str(x)"
585 return str(self.__A)+"\n"+str(self.__D)
587 # ==============================================================================
588 class DataObserver(object):
590 Classe générale d'interface de type observer
593 name = "GenericObserver",
600 self.__name = str(name)
605 if onVariable is None:
606 raise ValueError("setting an observer has to be done over a variable name, not over None.")
608 self.__V = str(onVariable)
611 self.__I = str(onVariable)
613 self.__I = str(withInfo)
615 if asString is not None:
616 __FunctionText = asString
617 elif (asTemplate is not None) and (asTemplate in ObserverTemplates):
618 __FunctionText = ObserverTemplates[asTemplate]
619 elif asScript is not None:
620 __FunctionText = ImportFromScript(asScript).getstring()
623 __Function = ObserverF(__FunctionText)
625 self.__O = __Function.getfunc()
627 return {self.__V:(self.__O, self.__I)}
630 "x.__repr__() <==> repr(x)"
631 return repr(self.__V)+"\n"+repr(self.__O)
634 "x.__str__() <==> str(x)"
635 return str(self.__V)+"\n"+str(self.__O)
637 # ==============================================================================
640 Classe générale d'interface de type état
643 name = "GenericVector",
645 asPersistentVector = None,
651 Permet de définir un vecteur :
652 - asVector : entrée des données, comme un vecteur compatible avec le
653 constructeur de numpy.matrix, ou "True" si entrée par script.
654 - asPersistentVector : entrée des données, comme une série de vecteurs
655 compatible avec le constructeur de numpy.matrix, ou comme un objet de
656 type Persistence, ou "True" si entrée par script.
657 - asScript : si un script valide est donné contenant une variable
658 nommée "name", la variable est de type "asVector" (par défaut) ou
659 "asPersistentVector" selon que l'une de ces variables est placée à
662 self.__name = str(name)
663 self.__check = bool(toBeChecked)
667 self.__is_vector = False
668 self.__is_series = False
670 if asScript is not None:
671 __Vector, __Series = None, None
672 if asPersistentVector:
673 __Series = ImportFromScript(asScript).getvalue( self.__name )
675 __Vector = ImportFromScript(asScript).getvalue( self.__name )
677 __Vector, __Series = asVector, asPersistentVector
679 if __Vector is not None:
680 self.__is_vector = True
681 self.__V = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
682 self.shape = self.__V.shape
683 self.size = self.__V.size
684 elif __Series is not None:
685 self.__is_series = True
686 if type(__Series) in [tuple, list, numpy.ndarray, numpy.matrix]:
687 self.__V = Persistence.OneVector(self.__name, basetype=numpy.matrix)
688 for member in __Series:
689 self.__V.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
690 import sys ; sys.stdout.flush()
693 if type(self.__V.shape) in (tuple, list):
694 self.shape = self.__V.shape
696 self.shape = self.__V.shape()
697 self.size = self.shape[0] * self.shape[1]
699 raise ValueError("The %s object is improperly defined, it requires at minima either a vector, a list/tuple of vectors or a persistent object. Please check your vector input."%self.__name)
701 if Scheduler is not None:
705 "Vérification du type interne"
706 return self.__is_vector
709 "Vérification du type interne"
710 return self.__is_series
713 "x.__repr__() <==> repr(x)"
714 return repr(self.__V)
717 "x.__str__() <==> str(x)"
720 # ==============================================================================
721 class Covariance(object):
723 Classe générale d'interface de type covariance
726 name = "GenericCovariance",
728 asEyeByScalar = None,
729 asEyeByVector = None,
735 Permet de définir une covariance :
736 - asCovariance : entrée des données, comme une matrice compatible avec
737 le constructeur de numpy.matrix
738 - asEyeByScalar : entrée des données comme un seul scalaire de variance,
739 multiplicatif d'une matrice de corrélation identité, aucune matrice
740 n'étant donc explicitement à donner
741 - asEyeByVector : entrée des données comme un seul vecteur de variance,
742 à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
743 n'étant donc explicitement à donner
744 - asCovObject : entrée des données comme un objet python, qui a les
745 methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
746 "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
747 "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
748 - toBeChecked : booléen indiquant si le caractère SDP de la matrice
749 pleine doit être vérifié
751 self.__name = str(name)
752 self.__check = bool(toBeChecked)
755 self.__is_scalar = False
756 self.__is_vector = False
757 self.__is_matrix = False
758 self.__is_object = False
760 if asScript is not None:
761 __Matrix, __Scalar, __Vector, __Object = None, None, None, None
763 __Scalar = _ImportFromScript(asScript).getvalue( "BackgroundError" )
765 __Vector = _ImportFromScript(asScript).getvalue( "BackgroundError" )
767 __Object = _ImportFromScript(asScript).getvalue( "BackgroundError" )
769 __Matrix = _ImportFromScript(asScript).getvalue( "BackgroundError" )
771 __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject
773 if __Scalar is not None:
774 if numpy.matrix(__Scalar).size != 1:
775 raise ValueError(' The diagonal multiplier given to define a sparse matrix is not a unique scalar value.\n Its actual measured size is %i. Please check your scalar input.'%numpy.matrix(__Scalar).size)
776 self.__is_scalar = True
777 self.__C = numpy.abs( float(__Scalar) )
780 elif __Vector is not None:
781 self.__is_vector = True
782 self.__C = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(__Vector, float ) ) ) )
783 self.shape = (self.__C.size,self.__C.size)
784 self.size = self.__C.size**2
785 elif __Matrix is not None:
786 self.__is_matrix = True
787 self.__C = numpy.matrix( __Matrix, float )
788 self.shape = self.__C.shape
789 self.size = self.__C.size
790 elif __Object is not None:
791 self.__is_object = True
793 for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
794 if not hasattr(self.__C,at):
795 raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
796 if hasattr(self.__C,"shape"):
797 self.shape = self.__C.shape
800 if hasattr(self.__C,"size"):
801 self.size = self.__C.size
806 # raise ValueError("The %s covariance matrix has to be specified either as a matrix, a vector for its diagonal or a scalar multiplying an identity matrix."%self.__name)
810 def __validate(self):
812 if self.ismatrix() and min(self.shape) != max(self.shape):
813 raise ValueError("The given matrix for %s is not a square one, its shape is %s. Please check your matrix input."%(self.__name,self.shape))
814 if self.isobject() and min(self.shape) != max(self.shape):
815 raise ValueError("The matrix given for \"%s\" is not a square one, its shape is %s. Please check your object input."%(self.__name,self.shape))
816 if self.isscalar() and self.__C <= 0:
817 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
818 if self.isvector() and (self.__C <= 0).any():
819 raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
820 if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
822 L = numpy.linalg.cholesky( self.__C )
824 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
825 if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
827 L = self.__C.cholesky()
829 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
832 "Vérification du type interne"
833 return self.__is_scalar
836 "Vérification du type interne"
837 return self.__is_vector
840 "Vérification du type interne"
841 return self.__is_matrix
844 "Vérification du type interne"
845 return self.__is_object
850 return Covariance(self.__name+"I", asCovariance = self.__C.I )
851 elif self.isvector():
852 return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
853 elif self.isscalar():
854 return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
855 elif self.isobject():
856 return Covariance(self.__name+"I", asCovObject = self.__C.getI() )
858 return None # Indispensable
863 return Covariance(self.__name+"T", asCovariance = self.__C.T )
864 elif self.isvector():
865 return Covariance(self.__name+"T", asEyeByVector = self.__C )
866 elif self.isscalar():
867 return Covariance(self.__name+"T", asEyeByScalar = self.__C )
868 elif self.isobject():
869 return Covariance(self.__name+"T", asCovObject = self.__C.getT() )
872 "Décomposition de Cholesky"
874 return Covariance(self.__name+"C", asCovariance = numpy.linalg.cholesky(self.__C) )
875 elif self.isvector():
876 return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
877 elif self.isscalar():
878 return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
879 elif self.isobject() and hasattr(self.__C,"cholesky"):
880 return Covariance(self.__name+"C", asCovObject = self.__C.cholesky() )
883 "Inversion de la décomposition de Cholesky"
885 return Covariance(self.__name+"H", asCovariance = numpy.linalg.cholesky(self.__C).I )
886 elif self.isvector():
887 return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
888 elif self.isscalar():
889 return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
890 elif self.isobject() and hasattr(self.__C,"choleskyI"):
891 return Covariance(self.__name+"H", asCovObject = self.__C.choleskyI() )
893 def diag(self, msize=None):
894 "Diagonale de la matrice"
896 return numpy.diag(self.__C)
897 elif self.isvector():
899 elif self.isscalar():
901 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
903 return self.__C * numpy.ones(int(msize))
904 elif self.isobject():
905 return self.__C.diag()
907 def asfullmatrix(self, msize=None):
911 elif self.isvector():
912 return numpy.matrix( numpy.diag(self.__C), float )
913 elif self.isscalar():
915 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
917 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
918 elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
919 return self.__C.asfullmatrix()
921 def trace(self, msize=None):
922 "Trace de la matrice"
924 return numpy.trace(self.__C)
925 elif self.isvector():
926 return float(numpy.sum(self.__C))
927 elif self.isscalar():
929 raise ValueError("the size of the %s covariance matrix has to be given in case of definition as a scalar over the diagonal."%(self.__name,))
931 return self.__C * int(msize)
932 elif self.isobject():
933 return self.__C.trace()
936 "x.__repr__() <==> repr(x)"
937 return repr(self.__C)
940 "x.__str__() <==> str(x)"
943 def __add__(self, other):
944 "x.__add__(y) <==> x+y"
945 if self.ismatrix() or self.isobject():
946 return self.__C + numpy.asmatrix(other)
947 elif self.isvector() or self.isscalar():
948 _A = numpy.asarray(other)
949 _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
950 return numpy.asmatrix(_A)
952 def __radd__(self, other):
953 "x.__radd__(y) <==> y+x"
954 raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
956 def __sub__(self, other):
957 "x.__sub__(y) <==> x-y"
958 if self.ismatrix() or self.isobject():
959 return self.__C - numpy.asmatrix(other)
960 elif self.isvector() or self.isscalar():
961 _A = numpy.asarray(other)
962 _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
963 return numpy.asmatrix(_A)
965 def __rsub__(self, other):
966 "x.__rsub__(y) <==> y-x"
967 raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
970 "x.__neg__() <==> -x"
973 def __mul__(self, other):
974 "x.__mul__(y) <==> x*y"
975 if self.ismatrix() and isinstance(other,numpy.matrix):
976 return self.__C * other
977 elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
978 or isinstance(other,list) \
979 or isinstance(other,tuple)):
980 if numpy.ravel(other).size == self.shape[1]: # Vecteur
981 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
982 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
983 return self.__C * numpy.asmatrix(other)
985 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
986 elif self.isvector() and (isinstance(other,numpy.matrix) \
987 or isinstance(other,numpy.ndarray) \
988 or isinstance(other,list) \
989 or isinstance(other,tuple)):
990 if numpy.ravel(other).size == self.shape[1]: # Vecteur
991 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
992 elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
993 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
995 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
996 elif self.isscalar() and isinstance(other,numpy.matrix):
997 return self.__C * other
998 elif self.isscalar() and (isinstance(other,numpy.ndarray) \
999 or isinstance(other,list) \
1000 or isinstance(other,tuple)):
1001 if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
1002 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
1004 return self.__C * numpy.asmatrix(other)
1005 elif self.isobject():
1006 return self.__C.__mul__(other)
1008 raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
1010 def __rmul__(self, other):
1011 "x.__rmul__(y) <==> y*x"
1012 if self.ismatrix() and isinstance(other,numpy.matrix):
1013 return other * self.__C
1014 elif self.isvector() and isinstance(other,numpy.matrix):
1015 if numpy.ravel(other).size == self.shape[0]: # Vecteur
1016 return numpy.asmatrix(numpy.ravel(other) * self.__C)
1017 elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
1018 return numpy.asmatrix(numpy.array(other) * self.__C)
1020 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
1021 elif self.isscalar() and isinstance(other,numpy.matrix):
1022 return other * self.__C
1023 elif self.isobject():
1024 return self.__C.__rmul__(other)
1026 raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
1029 "x.__len__() <==> len(x)"
1030 return self.shape[0]
1032 # ==============================================================================
1033 class ObserverF(object):
1035 Creation d'une fonction d'observateur a partir de son texte
1037 def __init__(self, corps=""):
1038 self.__corps = corps
1039 def func(self,var,info):
1040 "Fonction d'observation"
1043 "Restitution du pointeur de fonction dans l'objet"
1046 # ==============================================================================
1047 class ImportFromScript(object):
1049 Obtention d'une variable nommee depuis un fichier script importe
1051 def __init__(self, __filename=None):
1052 "Verifie l'existence et importe le script"
1053 __filename = __filename.rstrip(".py")
1054 if __filename is None:
1055 raise ValueError("The name of the file containing the variable to be imported has to be specified.")
1056 if not os.path.isfile(str(__filename)+".py"):
1057 raise ValueError("The file containing the variable to be imported doesn't seem to exist. The given file name is:\n \"%s\""%__filename)
1058 self.__scriptfile = __import__(__filename, globals(), locals(), [])
1059 self.__scriptstring = open(__filename+".py",'r').read()
1060 def getvalue(self, __varname=None, __synonym=None ):
1061 "Renvoie la variable demandee"
1062 if __varname is None:
1063 raise ValueError("The name of the variable to be imported has to be specified.")
1064 if not hasattr(self.__scriptfile, __varname):
1065 if __synonym is None:
1066 raise ValueError("The imported script file doesn't contain the specified variable \"%s\"."%__varname)
1067 elif not hasattr(self.__scriptfile, __synonym):
1068 raise ValueError("The imported script file doesn't contain the specified variable \"%s\"."%__synonym)
1070 return getattr(self.__scriptfile, __synonym)
1072 return getattr(self.__scriptfile, __varname)
1073 def getstring(self):
1074 "Renvoie le script complet"
1075 return self.__scriptstring
1077 # ==============================================================================
1078 def CostFunction3D(_x,
1079 _Hm = None, # Pour simuler Hm(x) : HO["Direct"].appliedTo
1080 _HmX = None, # Simulation déjà faite de Hm(x)
1081 _arg = None, # Arguments supplementaires pour Hm, sous la forme d'un tuple
1086 _SIV = False, # A résorber pour la 8.0
1087 _SSC = [], # self._parameters["StoreSupplementaryCalculations"]
1088 _nPS = 0, # nbPreviousSteps
1089 _QM = "DA", # QualityMeasure
1090 _SSV = {}, # Entrée et/ou sortie : self.StoredVariables
1091 _fRt = False, # Restitue ou pas la sortie étendue
1092 _sSc = True, # Stocke ou pas les SSC
1095 Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
1096 et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
1097 DFO, QuantileRegression
1103 for k in ["CostFunctionJ",
1109 "SimulatedObservationAtCurrentOptimum",
1110 "SimulatedObservationAtCurrentState",
1114 if hasattr(_SSV[k],"store"):
1115 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
1117 _X = numpy.asmatrix(numpy.ravel( _x )).T
1118 if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
1119 _SSV["CurrentState"].append( _X )
1121 if _HmX is not None:
1125 raise ValueError("COSTFUNCTION3D Operator has to be defined.")
1129 _HX = _Hm( _X, *_arg )
1130 _HX = numpy.asmatrix(numpy.ravel( _HX )).T
1132 if "SimulatedObservationAtCurrentState" in _SSC or \
1133 "SimulatedObservationAtCurrentOptimum" in _SSC:
1134 _SSV["SimulatedObservationAtCurrentState"].append( _HX )
1136 if numpy.any(numpy.isnan(_HX)):
1137 Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
1139 _Y = numpy.asmatrix(numpy.ravel( _Y )).T
1140 if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
1141 if _BI is None or _RI is None:
1142 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
1143 _Xb = numpy.asmatrix(numpy.ravel( _Xb )).T
1144 Jb = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
1145 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1146 elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
1148 raise ValueError("Observation error covariance matrix has to be properly defined!")
1150 Jo = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1151 elif _QM in ["LeastSquares", "LS", "L2"]:
1153 Jo = 0.5 * (_Y - _HX).T * (_Y - _HX)
1154 elif _QM in ["AbsoluteValue", "L1"]:
1156 Jo = numpy.sum( numpy.abs(_Y - _HX) )
1157 elif _QM in ["MaximumError", "ME"]:
1159 Jo = numpy.max( numpy.abs(_Y - _HX) )
1160 elif _QM in ["QR", "Null"]:
1164 raise ValueError("Unknown asked quality measure!")
1166 J = float( Jb ) + float( Jo )
1169 _SSV["CostFunctionJb"].append( Jb )
1170 _SSV["CostFunctionJo"].append( Jo )
1171 _SSV["CostFunctionJ" ].append( J )
1173 if "IndexOfOptimum" in _SSC or \
1174 "CurrentOptimum" in _SSC or \
1175 "SimulatedObservationAtCurrentOptimum" in _SSC:
1176 IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
1177 if "IndexOfOptimum" in _SSC:
1178 _SSV["IndexOfOptimum"].append( IndexMin )
1179 if "CurrentOptimum" in _SSC:
1180 _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
1181 if "SimulatedObservationAtCurrentOptimum" in _SSC:
1182 _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
1187 if _QM in ["QR"]: # Pour le QuantileRegression
1192 # ==============================================================================
1193 if __name__ == "__main__":
1194 print('\n AUTODIAGNOSTIC \n')