Salome HOME
5e9b4e481b8970dba565fc482e1fe875528bfa50
[modules/adao.git] / src / daComposant / daCore / BasicObjects.py
1 #-*-coding:iso-8859-1-*-
2 #
3 # Copyright (C) 2008-2017 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 """
24     Définit les outils généraux élémentaires.
25
26     Ce module est destiné à être appelée par AssimilationStudy.
27 """
28 __author__ = "Jean-Philippe ARGAUD"
29 __all__ = []
30
31 import logging, copy
32 import numpy
33 from daCore import Persistence
34 from daCore import PlatformInfo
35
36 # ==============================================================================
37 class CacheManager(object):
38     """
39     Classe générale de gestion d'un cache de calculs
40     """
41     def __init__(self,
42                  toleranceInRedundancy = 1.e-18,
43                  lenghtOfRedundancy    = -1,
44                 ):
45         """
46         Les caractéristiques de tolérance peuvent être modifées à la création.
47         """
48         self.__tolerBP  = float(toleranceInRedundancy)
49         self.__lenghtOR = int(lenghtOfRedundancy)
50         self.__initlnOR = self.__lenghtOR
51         self.clearCache()
52
53     def clearCache(self):
54         "Vide le cache"
55         self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
56         # logging.debug("CM Tolerance de determination des doublons : %.2e", self.__tolerBP)
57
58     def wasCalculatedIn(self, xValue ): #, info="" ):
59         "Vérifie l'existence d'un calcul correspondant à la valeur"
60         __alc = False
61         __HxV = None
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)
65                 continue
66             if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
67                 __alc  = True
68                 __HxV = self.__listOPCV[i][1]
69                 # logging.debug("CM Cas%s déja calculé, portant le numéro %i", info, i)
70                 break
71         return __alc, __HxV
72
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)),
83             copy.copy(HxValue),
84             numpy.linalg.norm(xValue),
85             ) )
86
87     def disable(self):
88         "Inactive le cache"
89         self.__initlnOR = self.__lenghtOR
90         self.__lenghtOR = 0
91
92     def enable(self):
93         "Active le cache"
94         self.__lenghtOR = self.__initlnOR
95
96 # ==============================================================================
97 class Operator(object):
98     """
99     Classe générale d'interface de type opérateur
100     """
101     NbCallsAsMatrix = 0
102     NbCallsAsMethod = 0
103     NbCallsOfCached = 0
104     CM = CacheManager()
105     #
106     def __init__(self, fromMethod=None, fromMatrix=None, avoidingRedundancy = True):
107         """
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.
110         Arguments :
111         - fromMethod : argument de type fonction Python
112         - fromMatrix : argument adapté au constructeur numpy.matrix
113         - avoidingRedundancy : évite ou pas les calculs redondants
114         """
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
119             self.__Matrix = None
120             self.__Type   = "Method"
121         elif fromMatrix is not None:
122             self.__Method = None
123             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
124             self.__Type   = "Matrix"
125         else:
126             self.__Method = None
127             self.__Matrix = None
128             self.__Type   = None
129
130     def disableAvoidingRedundancy(self):
131         "Inactive le cache"
132         Operator.CM.disable()
133
134     def enableAvoidingRedundancy(self):
135         "Active le cache"
136         if self.__AvoidRC:
137             Operator.CM.enable()
138         else:
139             Operator.CM.disable()
140
141     def isType(self):
142         "Renvoie le type"
143         return self.__Type
144
145     def appliedTo(self, xValue, HValue = None):
146         """
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.
150         Arguments :
151         - xValue : argument adapté pour appliquer l'opérateur
152         """
153         if HValue is not None:
154             HxValue = numpy.asmatrix( numpy.ravel( HValue ) ).T
155             if self.__AvoidRC:
156                 Operator.CM.storeValueInX(xValue,HxValue)
157         else:
158             if self.__AvoidRC:
159                 __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xValue)
160             else:
161                 __alreadyCalculated = False
162             #
163             if __alreadyCalculated:
164                 self.__addOneCacheCall()
165                 HxValue = __HxV
166             else:
167                 if self.__Matrix is not None:
168                     self.__addOneMatrixCall()
169                     HxValue = self.__Matrix * xValue
170                 else:
171                     self.__addOneMethodCall()
172                     HxValue = self.__Method( xValue )
173                 if self.__AvoidRC:
174                     Operator.CM.storeValueInX(xValue,HxValue)
175         #
176         return HxValue
177
178     def appliedControledFormTo(self, paire ):
179         """
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.
184         Arguments :
185         - xValue : argument X adapté pour appliquer l'opérateur
186         - uValue : argument U adapté pour appliquer l'opérateur
187         """
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) )
196         else:
197             self.__addOneMethodCall()
198             return self.__Method( xValue )
199
200     def appliedInXTo(self, paire ):
201         """
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
207         d'utiliser xNominal.
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
212         """
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
218         else:
219             self.__addOneMethodCall()
220             return self.__Method( (xNominal, xValue) )
221
222     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue"):
223         """
224         Permet de renvoyer l'opérateur sous la forme d'une matrice
225         """
226         if self.__Matrix is not None:
227             self.__addOneMatrixCall()
228             return self.__Matrix
229         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
230             self.__addOneMethodCall()
231             return numpy.matrix( self.__Method( (ValueForMethodForm, None) ) )
232         else:
233             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
234
235     def shape(self):
236         """
237         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
238         la forme d'une matrice
239         """
240         if self.__Matrix is not None:
241             return self.__Matrix.shape
242         else:
243             raise ValueError("Matrix form of the operator is not available, nor the shape")
244
245     def nbcalls(self, which=None):
246         """
247         Renvoie les nombres d'évaluations de l'opérateur
248         """
249         __nbcalls = (
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,
258             )
259         if which is None: return __nbcalls
260         else:             return __nbcalls[which]
261
262     def __addOneMatrixCall(self):
263         "Comptabilise un appel"
264         self.__NbCallsAsMatrix   += 1 # Decompte local
265         Operator.NbCallsAsMatrix += 1 # Decompte global
266
267     def __addOneMethodCall(self):
268         "Comptabilise un appel"
269         self.__NbCallsAsMethod   += 1 # Decompte local
270         Operator.NbCallsAsMethod += 1 # Decompte global
271
272     def __addOneCacheCall(self):
273         "Comptabilise un appel"
274         self.__NbCallsOfCached   += 1 # Decompte local
275         Operator.NbCallsOfCached += 1 # Decompte global
276
277 # ==============================================================================
278 class Algorithm(object):
279     """
280     Classe générale d'interface de type algorithme
281
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.
285
286     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
287     """
288     def __init__(self, name):
289         """
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".
294
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
322         """
323         logging.debug("%s Initialisation", str(name))
324         self._m = PlatformInfo.SystemUsage()
325         #
326         self._name = str( name )
327         self._parameters = {"StoreSupplementaryCalculations":[]}
328         self.__required_parameters = {}
329         self.StoredVariables = {}
330         #
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")
362
363     def _pre_run(self, Parameters ):
364         "Pré-calcul"
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"))
367         #
368         # Mise a jour de self._parameters avec Parameters
369         self.__setParameters(Parameters)
370         #
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,))
374         else:
375             self._parameters["Bounds"] = None
376         #
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
382             else:
383                 self._parameters["optmessages"] = 15
384         else:
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
389             else:
390                 self._parameters["optmessages"] = 15
391         #
392         return 0
393
394     def _post_run(self,_oH=None):
395         "Post-calcul"
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 )
407         if _oH is not None:
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)
412         return 0
413
414     def get(self, key=None):
415         """
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.
421         """
422         if key is not None:
423             return self.StoredVariables[key]
424         else:
425             return self.StoredVariables
426
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
430
431     def keys(self):
432         "D.keys() -> list of D's keys"
433         return self.StoredVariables.keys()
434
435     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
436         """
437         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
438         sa forme mathématique la plus naturelle possible.
439         """
440         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
441
442     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
443         """
444         Permet de définir dans l'algorithme des paramètres requis et leurs
445         caractéristiques par défaut.
446         """
447         if name is None:
448             raise ValueError("A name is mandatory to define a required parameter.")
449         #
450         self.__required_parameters[name] = {
451             "default"  : default,
452             "typecast" : typecast,
453             "minval"   : minval,
454             "maxval"   : maxval,
455             "listval"  : listval,
456             "message"  : message,
457             }
458         logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
459
460     def getRequiredParameters(self, noDetails=True):
461         """
462         Renvoie la liste des noms de paramètres requis ou directement le
463         dictionnaire des paramètres requis.
464         """
465         if noDetails:
466             ks = list(self.__required_parameters.keys())
467             ks.sort()
468             return ks
469         else:
470             return self.__required_parameters
471
472     def setParameterValue(self, name=None, value=None):
473         """
474         Renvoie la valeur d'un paramètre requis de manière contrôlée
475         """
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"]
481         #
482         if value is None and default is None:
483             __val = None
484         elif value is None and default is not None:
485             if typecast is None: __val = default
486             else:                __val = typecast( default )
487         else:
488             if typecast is None: __val = value
489             else:                __val = typecast( value )
490         #
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):
497                 for v in __val:
498                     if v not in listval:
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))
502         return __val
503
504     def __setParameters(self, fromDico={}):
505         """
506         Permet de stocker les paramètres reçus dans le dictionnaire interne.
507         """
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])
512             else:
513                 self._parameters[k] = self.setParameterValue(k)
514             logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
515
516 # ==============================================================================
517 class Diagnostic(object):
518     """
519     Classe générale d'interface de type diagnostic
520
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
523     exemple.
524
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).
531     """
532     def __init__(self, name = "", parameters = {}):
533         "Initialisation"
534         self.name       = str(name)
535         self.parameters = dict( parameters )
536
537     def _formula(self, *args):
538         """
539         Doit implémenter l'opération élémentaire de diagnostic sous sa forme
540         mathématique la plus naturelle possible.
541         """
542         raise NotImplementedError("Diagnostic mathematical formula has not been implemented!")
543
544     def calculate(self, *args):
545         """
546         Active la formule de calcul avec les arguments correctement rangés
547         """
548         raise NotImplementedError("Diagnostic activation method has not been implemented!")
549
550 # ==============================================================================
551 class Vector(object):
552     """
553     Classe générale d'interface de type vecteur
554     """
555     def __init__(self,
556                  name               = "GenericVector",
557                  asVector           = None,
558                  asPersistentVector = None,
559                  Scheduler          = None,
560                  toBeChecked        = False,
561                 ):
562         """
563         Permet de définir un vecteur :
564         - asVector : entrée des données, comme un vecteur compatible avec
565           le constructeur de numpy.matrix
566         - asPersistentVector : entrée des données, comme une série de
567           vecteurs compatible avec le constructeur de numpy.matrix, ou
568           comme un objet de type Persistence
569         """
570         self.__name       = str(name)
571         self.__check      = bool(toBeChecked)
572         #
573         self.__V          = None
574         self.__T          = None
575         self.__is_vector  = False
576         self.__is_series  = False
577         if asVector is not None:
578             self.__is_vector = True
579             self.__V         = numpy.matrix( numpy.ravel(numpy.matrix(asVector)), numpy.float ).T
580             self.shape       = self.__V.shape
581             self.size        = self.__V.size
582         elif asPersistentVector is not None:
583             self.__is_series  = True
584             if type(asPersistentVector) in [tuple, list, numpy.ndarray, numpy.matrix]:
585                 self.__V = Persistence.OneVector(self.__name, basetype=numpy.matrix)
586                 for member in asPersistentVector:
587                     self.__V.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
588                 import sys ; sys.stdout.flush()
589             else:
590                 self.__V = asPersistentVector
591             if type(self.__V.shape) in (tuple, list):
592                 self.shape       = self.__V.shape
593             else:
594                 self.shape       = self.__V.shape()
595             self.size        = self.shape[0] * self.shape[1]
596         else:
597             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)
598         #
599         if Scheduler is not None:
600             self.__T = Scheduler
601
602     def isvector(self):
603         "Vérification du type interne"
604         return self.__is_vector
605
606     def isseries(self):
607         "Vérification du type interne"
608         return self.__is_series
609
610     def __repr__(self):
611         "x.__repr__() <==> repr(x)"
612         return repr(self.__V)
613
614     def __str__(self):
615         "x.__str__() <==> str(x)"
616         return str(self.__V)
617
618 # ==============================================================================
619 class Covariance(object):
620     """
621     Classe générale d'interface de type covariance
622     """
623     def __init__(self,
624                  name          = "GenericCovariance",
625                  asCovariance  = None,
626                  asEyeByScalar = None,
627                  asEyeByVector = None,
628                  asCovObject   = None,
629                  toBeChecked   = False,
630                 ):
631         """
632         Permet de définir une covariance :
633         - asCovariance : entrée des données, comme une matrice compatible avec
634           le constructeur de numpy.matrix
635         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
636           multiplicatif d'une matrice de corrélation identité, aucune matrice
637           n'étant donc explicitement à donner
638         - asEyeByVector : entrée des données comme un seul vecteur de variance,
639           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
640           n'étant donc explicitement à donner
641         - asCovObject : entrée des données comme un objet python, qui a les
642           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
643           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
644           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
645         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
646           pleine doit être vérifié
647         """
648         self.__name       = str(name)
649         self.__check      = bool(toBeChecked)
650         #
651         self.__C          = None
652         self.__is_scalar  = False
653         self.__is_vector  = False
654         self.__is_matrix  = False
655         self.__is_object  = False
656         if asEyeByScalar is not None:
657             if numpy.matrix(asEyeByScalar).size != 1:
658                 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(asEyeByScalar).size)
659             self.__is_scalar = True
660             self.__C         = numpy.abs( float(asEyeByScalar) )
661             self.shape       = (0,0)
662             self.size        = 0
663         elif asEyeByVector is not None:
664             self.__is_vector = True
665             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(asEyeByVector, float ) ) ) )
666             self.shape       = (self.__C.size,self.__C.size)
667             self.size        = self.__C.size**2
668         elif asCovariance is not None:
669             self.__is_matrix = True
670             self.__C         = numpy.matrix( asCovariance, float )
671             self.shape       = self.__C.shape
672             self.size        = self.__C.size
673         elif asCovObject is not None:
674             self.__is_object = True
675             self.__C         = asCovObject
676             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
677                 if not hasattr(self.__C,at):
678                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
679             if hasattr(self.__C,"shape"):
680                 self.shape       = self.__C.shape
681             else:
682                 self.shape       = (0,0)
683             if hasattr(self.__C,"size"):
684                 self.size        = self.__C.size
685             else:
686                 self.size        = 0
687         else:
688             pass
689             # 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)
690         #
691         self.__validate()
692
693     def __validate(self):
694         "Validation"
695         if self.ismatrix() and min(self.shape) != max(self.shape):
696             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))
697         if self.isobject() and min(self.shape) != max(self.shape):
698             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))
699         if self.isscalar() and self.__C <= 0:
700             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
701         if self.isvector() and (self.__C <= 0).any():
702             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
703         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
704             try:
705                 L = numpy.linalg.cholesky( self.__C )
706             except:
707                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
708         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
709             try:
710                 L = self.__C.cholesky()
711             except:
712                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
713
714     def isscalar(self):
715         "Vérification du type interne"
716         return self.__is_scalar
717
718     def isvector(self):
719         "Vérification du type interne"
720         return self.__is_vector
721
722     def ismatrix(self):
723         "Vérification du type interne"
724         return self.__is_matrix
725
726     def isobject(self):
727         "Vérification du type interne"
728         return self.__is_object
729
730     def getI(self):
731         "Inversion"
732         if   self.ismatrix():
733             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
734         elif self.isvector():
735             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
736         elif self.isscalar():
737             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
738         elif self.isobject():
739             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
740         else:
741             return None # Indispensable
742
743     def getT(self):
744         "Transposition"
745         if   self.ismatrix():
746             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
747         elif self.isvector():
748             return Covariance(self.__name+"T", asEyeByVector = self.__C )
749         elif self.isscalar():
750             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
751         elif self.isobject():
752             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
753
754     def cholesky(self):
755         "Décomposition de Cholesky"
756         if   self.ismatrix():
757             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
758         elif self.isvector():
759             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
760         elif self.isscalar():
761             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
762         elif self.isobject() and hasattr(self.__C,"cholesky"):
763             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
764
765     def choleskyI(self):
766         "Inversion de la décomposition de Cholesky"
767         if   self.ismatrix():
768             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
769         elif self.isvector():
770             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
771         elif self.isscalar():
772             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
773         elif self.isobject() and hasattr(self.__C,"choleskyI"):
774             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
775
776     def diag(self, msize=None):
777         "Diagonale de la matrice"
778         if   self.ismatrix():
779             return numpy.diag(self.__C)
780         elif self.isvector():
781             return self.__C
782         elif self.isscalar():
783             if msize is None:
784                 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,))
785             else:
786                 return self.__C * numpy.ones(int(msize))
787         elif self.isobject():
788             return self.__C.diag()
789
790     def asfullmatrix(self, msize=None):
791         "Matrice pleine"
792         if   self.ismatrix():
793             return self.__C
794         elif self.isvector():
795             return numpy.matrix( numpy.diag(self.__C), float )
796         elif self.isscalar():
797             if msize is None:
798                 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,))
799             else:
800                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
801         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
802             return self.__C.asfullmatrix()
803
804     def trace(self, msize=None):
805         "Trace de la matrice"
806         if   self.ismatrix():
807             return numpy.trace(self.__C)
808         elif self.isvector():
809             return float(numpy.sum(self.__C))
810         elif self.isscalar():
811             if msize is None:
812                 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,))
813             else:
814                 return self.__C * int(msize)
815         elif self.isobject():
816             return self.__C.trace()
817
818     def __repr__(self):
819         "x.__repr__() <==> repr(x)"
820         return repr(self.__C)
821
822     def __str__(self):
823         "x.__str__() <==> str(x)"
824         return str(self.__C)
825
826     def __add__(self, other):
827         "x.__add__(y) <==> x+y"
828         if   self.ismatrix() or self.isobject():
829             return self.__C + numpy.asmatrix(other)
830         elif self.isvector() or self.isscalar():
831             _A = numpy.asarray(other)
832             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
833             return numpy.asmatrix(_A)
834
835     def __radd__(self, other):
836         "x.__radd__(y) <==> y+x"
837         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
838
839     def __sub__(self, other):
840         "x.__sub__(y) <==> x-y"
841         if   self.ismatrix() or self.isobject():
842             return self.__C - numpy.asmatrix(other)
843         elif self.isvector() or self.isscalar():
844             _A = numpy.asarray(other)
845             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
846             return numpy.asmatrix(_A)
847
848     def __rsub__(self, other):
849         "x.__rsub__(y) <==> y-x"
850         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
851
852     def __neg__(self):
853         "x.__neg__() <==> -x"
854         return - self.__C
855
856     def __mul__(self, other):
857         "x.__mul__(y) <==> x*y"
858         if   self.ismatrix() and isinstance(other,numpy.matrix):
859             return self.__C * other
860         elif self.ismatrix() and (isinstance(other,numpy.ndarray) \
861                                or isinstance(other,list) \
862                                or isinstance(other,tuple)):
863             if numpy.ravel(other).size == self.shape[1]: # Vecteur
864                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
865             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
866                 return self.__C * numpy.asmatrix(other)
867             else:
868                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
869         elif self.isvector() and (isinstance(other,numpy.matrix) \
870                                or isinstance(other,numpy.ndarray) \
871                                or isinstance(other,list) \
872                                or isinstance(other,tuple)):
873             if numpy.ravel(other).size == self.shape[1]: # Vecteur
874                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
875             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
876                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
877             else:
878                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
879         elif self.isscalar() and isinstance(other,numpy.matrix):
880             return self.__C * other
881         elif self.isscalar() and (isinstance(other,numpy.ndarray) \
882                                or isinstance(other,list) \
883                                or isinstance(other,tuple)):
884             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
885                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
886             else:
887                 return self.__C * numpy.asmatrix(other)
888         elif self.isobject():
889             return self.__C.__mul__(other)
890         else:
891             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
892
893     def __rmul__(self, other):
894         "x.__rmul__(y) <==> y*x"
895         if self.ismatrix() and isinstance(other,numpy.matrix):
896             return other * self.__C
897         elif self.isvector() and isinstance(other,numpy.matrix):
898             if numpy.ravel(other).size == self.shape[0]: # Vecteur
899                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
900             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
901                 return numpy.asmatrix(numpy.array(other) * self.__C)
902             else:
903                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
904         elif self.isscalar() and isinstance(other,numpy.matrix):
905             return other * self.__C
906         elif self.isobject():
907             return self.__C.__rmul__(other)
908         else:
909             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
910
911     def __len__(self):
912         "x.__len__() <==> len(x)"
913         return self.shape[0]
914
915 # ==============================================================================
916 def CostFunction3D(_x,
917                    _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
918                    _HmX = None,  # Simulation déjà faite de Hm(x)
919                    _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
920                    _BI  = None,
921                    _RI  = None,
922                    _Xb  = None,
923                    _Y   = None,
924                    _SIV = False, # A résorber pour la 8.0
925                    _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
926                    _nPS = 0,     # nbPreviousSteps
927                    _QM  = "DA",  # QualityMeasure
928                    _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
929                    _fRt = False, # Restitue ou pas la sortie étendue
930                    _sSc = True,  # Stocke ou pas les SSC
931                   ):
932     """
933     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
934     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
935     DFO, QuantileRegression
936     """
937     if not _sSc:
938         _SIV = False
939         _SSC = {}
940     else:
941         for k in ["CostFunctionJ",
942                   "CostFunctionJb",
943                   "CostFunctionJo",
944                   "CurrentOptimum",
945                   "CurrentState",
946                   "IndexOfOptimum",
947                   "SimulatedObservationAtCurrentOptimum",
948                   "SimulatedObservationAtCurrentState",
949                  ]:
950             if k not in _SSV:
951                 _SSV[k] = []
952             if hasattr(_SSV[k],"store"):
953                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
954     #
955     _X  = numpy.asmatrix(numpy.ravel( _x )).T
956     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
957         _SSV["CurrentState"].append( _X )
958     #
959     if _HmX is not None:
960         _HX = _HmX
961     else:
962         if _Hm is None:
963             raise ValueError("COSTFUNCTION3D Operator has to be defined.")
964         if _arg is None:
965             _HX = _Hm( _X )
966         else:
967             _HX = _Hm( _X, *_arg )
968     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
969     #
970     if "SimulatedObservationAtCurrentState" in _SSC or \
971        "SimulatedObservationAtCurrentOptimum" in _SSC:
972         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
973     #
974     if numpy.any(numpy.isnan(_HX)):
975         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
976     else:
977         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
978         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
979             if _BI is None or _RI is None:
980                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
981             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
982             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
983             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
984         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
985             if _RI is None:
986                 raise ValueError("Observation error covariance matrix has to be properly defined!")
987             Jb  = 0.
988             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
989         elif _QM in ["LeastSquares", "LS", "L2"]:
990             Jb  = 0.
991             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
992         elif _QM in ["AbsoluteValue", "L1"]:
993             Jb  = 0.
994             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
995         elif _QM in ["MaximumError", "ME"]:
996             Jb  = 0.
997             Jo  = numpy.max( numpy.abs(_Y - _HX) )
998         elif _QM in ["QR", "Null"]:
999             Jb  = 0.
1000             Jo  = 0.
1001         else:
1002             raise ValueError("Unknown asked quality measure!")
1003         #
1004         J   = float( Jb ) + float( Jo )
1005     #
1006     if _sSc:
1007         _SSV["CostFunctionJb"].append( Jb )
1008         _SSV["CostFunctionJo"].append( Jo )
1009         _SSV["CostFunctionJ" ].append( J )
1010     #
1011     if "IndexOfOptimum" in _SSC or \
1012        "CurrentOptimum" in _SSC or \
1013        "SimulatedObservationAtCurrentOptimum" in _SSC:
1014         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
1015     if "IndexOfOptimum" in _SSC:
1016         _SSV["IndexOfOptimum"].append( IndexMin )
1017     if "CurrentOptimum" in _SSC:
1018         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
1019     if "SimulatedObservationAtCurrentOptimum" in _SSC:
1020         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
1021     #
1022     if _fRt:
1023         return _SSV
1024     else:
1025         if _QM in ["QR"]: # Pour le QuantileRegression
1026             return _HX
1027         else:
1028             return J
1029
1030 # ==============================================================================
1031 if __name__ == "__main__":
1032     print('\n AUTODIAGNOSTIC \n')