Salome HOME
Updating test with unittest and extended user case definition
[modules/adao.git] / src / daComposant / daCore / BasicObjects.py
1 # -*- coding: utf-8 -*-
2 #
3 # Copyright (C) 2008-2019 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 __author__ = "Jean-Philippe ARGAUD"
27 __all__ = []
28
29 import os
30 import sys
31 import logging
32 import copy
33 import numpy
34 from functools import partial
35 from daCore import Persistence
36 from daCore import PlatformInfo
37 from daCore import Interfaces
38 from daCore import Templates
39 from daCore.Interfaces import ImportFromScript, ImportFromFile
40
41 # ==============================================================================
42 class CacheManager(object):
43     """
44     Classe générale de gestion d'un cache de calculs
45     """
46     def __init__(self,
47                  toleranceInRedundancy = 1.e-18,
48                  lenghtOfRedundancy    = -1,
49                 ):
50         """
51         Les caractéristiques de tolérance peuvent être modifées à la création.
52         """
53         self.__tolerBP  = float(toleranceInRedundancy)
54         self.__lenghtOR = int(lenghtOfRedundancy)
55         self.__initlnOR = self.__lenghtOR
56         self.clearCache()
57
58     def clearCache(self):
59         "Vide le cache"
60         self.__listOPCV = [] # Operator Previous Calculated Points, Results, Point Norms
61         # logging.debug("CM Tolerance de determination des doublons : %.2e", self.__tolerBP)
62
63     def wasCalculatedIn(self, xValue ): #, info="" ):
64         "Vérifie l'existence d'un calcul correspondant à la valeur"
65         __alc = False
66         __HxV = None
67         for i in range(min(len(self.__listOPCV),self.__lenghtOR)-1,-1,-1):
68             if not hasattr(xValue, 'size') or (xValue.size != self.__listOPCV[i][0].size):
69                 # 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)
70                 continue
71             if numpy.linalg.norm(numpy.ravel(xValue) - self.__listOPCV[i][0]) < self.__tolerBP * self.__listOPCV[i][2]:
72                 __alc  = True
73                 __HxV = self.__listOPCV[i][1]
74                 # logging.debug("CM Cas%s déja calculé, portant le numéro %i", info, i)
75                 break
76         return __alc, __HxV
77
78     def storeValueInX(self, xValue, HxValue ):
79         "Stocke un calcul correspondant à la valeur"
80         if self.__lenghtOR < 0:
81             self.__lenghtOR = 2 * xValue.size + 2
82             self.__initlnOR = self.__lenghtOR
83         while len(self.__listOPCV) > self.__lenghtOR:
84             # logging.debug("CM Réduction de la liste des cas à %i éléments par suppression du premier", self.__lenghtOR)
85             self.__listOPCV.pop(0)
86         self.__listOPCV.append( (
87             copy.copy(numpy.ravel(xValue)),
88             copy.copy(HxValue),
89             numpy.linalg.norm(xValue),
90             ) )
91
92     def disable(self):
93         "Inactive le cache"
94         self.__initlnOR = self.__lenghtOR
95         self.__lenghtOR = 0
96
97     def enable(self):
98         "Active le cache"
99         self.__lenghtOR = self.__initlnOR
100
101 # ==============================================================================
102 class Operator(object):
103     """
104     Classe générale d'interface de type opérateur simple
105     """
106     NbCallsAsMatrix = 0
107     NbCallsAsMethod = 0
108     NbCallsOfCached = 0
109     CM = CacheManager()
110     #
111     def __init__(self,
112         fromMethod           = None,
113         fromMatrix           = None,
114         avoidingRedundancy   = True,
115         inputAsMultiFunction = False,
116         extraArguments       = None,
117         ):
118         """
119         On construit un objet de ce type en fournissant, à l'aide de l'un des
120         deux mots-clé, soit une fonction ou un multi-fonction python, soit une
121         matrice.
122         Arguments :
123         - fromMethod : argument de type fonction Python
124         - fromMatrix : argument adapté au constructeur numpy.matrix
125         - avoidingRedundancy : booléen évitant (ou pas) les calculs redondants
126         - inputAsMultiFunction : booléen indiquant une fonction explicitement
127           définie (ou pas) en multi-fonction
128         - extraArguments : arguments supplémentaires passés à la fonction de
129           base et ses dérivées (tuple ou dictionnaire)
130         """
131         self.__NbCallsAsMatrix, self.__NbCallsAsMethod, self.__NbCallsOfCached = 0, 0, 0
132         self.__AvoidRC   = bool( avoidingRedundancy )
133         self.__inputAsMF = bool( inputAsMultiFunction )
134         self.__extraArgs = extraArguments
135         if   fromMethod is not None and self.__inputAsMF:
136             self.__Method = fromMethod # logtimer(fromMethod)
137             self.__Matrix = None
138             self.__Type   = "Method"
139         elif fromMethod is not None and not self.__inputAsMF:
140             self.__Method = partial( MultiFonction, _sFunction=fromMethod)
141             self.__Matrix = None
142             self.__Type   = "Method"
143         elif fromMatrix is not None:
144             self.__Method = None
145             self.__Matrix = numpy.matrix( fromMatrix, numpy.float )
146             self.__Type   = "Matrix"
147         else:
148             self.__Method = None
149             self.__Matrix = None
150             self.__Type   = None
151
152     def disableAvoidingRedundancy(self):
153         "Inactive le cache"
154         Operator.CM.disable()
155
156     def enableAvoidingRedundancy(self):
157         "Active le cache"
158         if self.__AvoidRC:
159             Operator.CM.enable()
160         else:
161             Operator.CM.disable()
162
163     def isType(self):
164         "Renvoie le type"
165         return self.__Type
166
167     def appliedTo(self, xValue, HValue = None, argsAsSerie = False):
168         """
169         Permet de restituer le résultat de l'application de l'opérateur à une
170         série d'arguments xValue. Cette méthode se contente d'appliquer, chaque
171         argument devant a priori être du bon type.
172         Arguments :
173         - les arguments par série sont :
174             - xValue : argument adapté pour appliquer l'opérateur
175             - HValue : valeur précalculée de l'opérateur en ce point
176         - argsAsSerie : indique si les arguments sont une mono ou multi-valeur
177         """
178         if argsAsSerie:
179             _xValue = xValue
180             _HValue = HValue
181         else:
182             _xValue = (xValue,)
183             if HValue is not None:
184                 _HValue = (HValue,)
185             else:
186                 _HValue = HValue
187         PlatformInfo.isIterable( _xValue, True, " in Operator.appliedTo" )
188         #
189         if _HValue is not None:
190             assert len(_xValue) == len(_HValue), "Incompatible number of elements in xValue and HValue"
191             HxValue = []
192             for i in range(len(_HValue)):
193                 HxValue.append( numpy.asmatrix( numpy.ravel( _HValue[i] ) ).T )
194                 if self.__AvoidRC:
195                     Operator.CM.storeValueInX(_xValue[i],HxValue[-1])
196         else:
197             HxValue = []
198             _xserie = []
199             _hindex = []
200             for i, xv in enumerate(_xValue):
201                 if self.__AvoidRC:
202                     __alreadyCalculated, __HxV = Operator.CM.wasCalculatedIn(xv)
203                 else:
204                     __alreadyCalculated = False
205                 #
206                 if __alreadyCalculated:
207                     self.__addOneCacheCall()
208                     _hv = __HxV
209                 else:
210                     if self.__Matrix is not None:
211                         self.__addOneMatrixCall()
212                         _hv = self.__Matrix * xv
213                     else:
214                         self.__addOneMethodCall()
215                         _xserie.append( xv )
216                         _hindex.append(  i )
217                         _hv = None
218                 HxValue.append( _hv )
219             #
220             if len(_xserie)>0 and self.__Matrix is None:
221                 if self.__extraArgs is None:
222                     _hserie = self.__Method( _xserie ) # Calcul MF
223                 else:
224                     _hserie = self.__Method( _xserie, self.__extraArgs ) # Calcul MF
225                 if not hasattr(_hserie, "pop"):
226                     raise TypeError("The user input multi-function doesn't seem to return sequence results, behaving like a mono-function. It has to be checked.")
227                 for i in _hindex:
228                     _xv = _xserie.pop(0)
229                     _hv = _hserie.pop(0)
230                     HxValue[i] = _hv
231                     if self.__AvoidRC:
232                         Operator.CM.storeValueInX(_xv,_hv)
233         #
234         if argsAsSerie: return HxValue
235         else:           return HxValue[-1]
236
237     def appliedControledFormTo(self, paires, argsAsSerie = False ):
238         """
239         Permet de restituer le résultat de l'application de l'opérateur à des
240         paires (xValue, uValue). Cette méthode se contente d'appliquer, son
241         argument devant a priori être du bon type. Si la uValue est None,
242         on suppose que l'opérateur ne s'applique qu'à xValue.
243         Arguments :
244         - paires : les arguments par paire sont :
245             - xValue : argument X adapté pour appliquer l'opérateur
246             - uValue : argument U adapté pour appliquer l'opérateur
247         - argsAsSerie : indique si l'argument est une mono ou multi-valeur
248         """
249         if argsAsSerie: _xuValue = paires
250         else:           _xuValue = (paires,)
251         PlatformInfo.isIterable( _xuValue, True, " in Operator.appliedControledFormTo" )
252         #
253         if self.__Matrix is not None:
254             HxValue = []
255             for paire in _xuValue:
256                 _xValue, _uValue = paire
257                 self.__addOneMatrixCall()
258                 HxValue.append( self.__Matrix * _xValue )
259         else:
260             HxValue = []
261             for paire in _xuValue:
262                 _xuValue = []
263                 _xValue, _uValue = paire
264                 if _uValue is not None:
265                     _xuValue.append( paire )
266                 else:
267                     _xuValue.append( _xValue )
268             self.__addOneMethodCall( len(_xuValue) )
269             if self.__extraArgs is None:
270                 HxValue = self.__Method( _xuValue ) # Calcul MF
271             else:
272                 HxValue = self.__Method( _xuValue, self.__extraArgs ) # Calcul MF
273         #
274         if argsAsSerie: return HxValue
275         else:           return HxValue[-1]
276
277     def appliedInXTo(self, paires, argsAsSerie = False ):
278         """
279         Permet de restituer le résultat de l'application de l'opérateur à une
280         série d'arguments xValue, sachant que l'opérateur est valable en
281         xNominal. Cette méthode se contente d'appliquer, son argument devant a
282         priori être du bon type. Si l'opérateur est linéaire car c'est une
283         matrice, alors il est valable en tout point nominal et xNominal peut
284         être quelconque. Il n'y a qu'une seule paire par défaut, et argsAsSerie
285         permet d'indiquer que l'argument est multi-paires.
286         Arguments :
287         - paires : les arguments par paire sont :
288             - xNominal : série d'arguments permettant de donner le point où
289               l'opérateur est construit pour être ensuite appliqué
290             - xValue : série d'arguments adaptés pour appliquer l'opérateur
291         - argsAsSerie : indique si l'argument est une mono ou multi-valeur
292         """
293         if argsAsSerie: _nxValue = paires
294         else:           _nxValue = (paires,)
295         PlatformInfo.isIterable( _nxValue, True, " in Operator.appliedInXTo" )
296         #
297         if self.__Matrix is not None:
298             HxValue = []
299             for paire in _nxValue:
300                 _xNominal, _xValue = paire
301                 self.__addOneMatrixCall()
302                 HxValue.append( self.__Matrix * _xValue )
303         else:
304             self.__addOneMethodCall( len(_nxValue) )
305             if self.__extraArgs is None:
306                 HxValue = self.__Method( _nxValue ) # Calcul MF
307             else:
308                 HxValue = self.__Method( _nxValue, self.__extraArgs ) # Calcul MF
309         #
310         if argsAsSerie: return HxValue
311         else:           return HxValue[-1]
312
313     def asMatrix(self, ValueForMethodForm = "UnknownVoidValue", argsAsSerie = False):
314         """
315         Permet de renvoyer l'opérateur sous la forme d'une matrice
316         """
317         if self.__Matrix is not None:
318             self.__addOneMatrixCall()
319             mValue = [self.__Matrix,]
320         elif ValueForMethodForm is not "UnknownVoidValue": # Ne pas utiliser "None"
321             mValue = []
322             if argsAsSerie:
323                 self.__addOneMethodCall( len(ValueForMethodForm) )
324                 for _vfmf in ValueForMethodForm:
325                     mValue.append( numpy.matrix( self.__Method(((_vfmf, None),)) ) )
326             else:
327                 self.__addOneMethodCall()
328                 mValue = self.__Method(((ValueForMethodForm, None),))
329         else:
330             raise ValueError("Matrix form of the operator defined as a function/method requires to give an operating point.")
331         #
332         if argsAsSerie: return mValue
333         else:           return mValue[-1]
334
335     def shape(self):
336         """
337         Renvoie la taille sous forme numpy si l'opérateur est disponible sous
338         la forme d'une matrice
339         """
340         if self.__Matrix is not None:
341             return self.__Matrix.shape
342         else:
343             raise ValueError("Matrix form of the operator is not available, nor the shape")
344
345     def nbcalls(self, which=None):
346         """
347         Renvoie les nombres d'évaluations de l'opérateur
348         """
349         __nbcalls = (
350             self.__NbCallsAsMatrix+self.__NbCallsAsMethod,
351             self.__NbCallsAsMatrix,
352             self.__NbCallsAsMethod,
353             self.__NbCallsOfCached,
354             Operator.NbCallsAsMatrix+Operator.NbCallsAsMethod,
355             Operator.NbCallsAsMatrix,
356             Operator.NbCallsAsMethod,
357             Operator.NbCallsOfCached,
358             )
359         if which is None: return __nbcalls
360         else:             return __nbcalls[which]
361
362     def __addOneMatrixCall(self):
363         "Comptabilise un appel"
364         self.__NbCallsAsMatrix   += 1 # Decompte local
365         Operator.NbCallsAsMatrix += 1 # Decompte global
366
367     def __addOneMethodCall(self, nb = 1):
368         "Comptabilise un appel"
369         self.__NbCallsAsMethod   += nb # Decompte local
370         Operator.NbCallsAsMethod += nb # Decompte global
371
372     def __addOneCacheCall(self):
373         "Comptabilise un appel"
374         self.__NbCallsOfCached   += 1 # Decompte local
375         Operator.NbCallsOfCached += 1 # Decompte global
376
377 # ==============================================================================
378 class FullOperator(object):
379     """
380     Classe générale d'interface de type opérateur complet
381     (Direct, Linéaire Tangent, Adjoint)
382     """
383     def __init__(self,
384                  name             = "GenericFullOperator",
385                  asMatrix         = None,
386                  asOneFunction    = None, # 1 Fonction
387                  asThreeFunctions = None, # 3 Fonctions in a dictionary
388                  asScript         = None, # 1 or 3 Fonction(s) by script
389                  asDict           = None, # Parameters
390                  appliedInX       = None,
391                  extraArguments   = None,
392                  avoidRC          = True,
393                  inputAsMF        = False,# Fonction(s) as Multi-Functions
394                  scheduledBy      = None,
395                  toBeChecked      = False,
396                  ):
397         ""
398         self.__name      = str(name)
399         self.__check     = bool(toBeChecked)
400         self.__extraArgs = extraArguments
401         #
402         self.__FO        = {}
403         #
404         __Parameters = {}
405         if (asDict is not None) and isinstance(asDict, dict):
406             __Parameters.update( asDict )
407             if "DifferentialIncrement" in asDict:
408                 __Parameters["withIncrement"]  = asDict["DifferentialIncrement"]
409             if "CenteredFiniteDifference" in asDict:
410                 __Parameters["withCenteredDF"] = asDict["CenteredFiniteDifference"]
411             if "EnableMultiProcessing" in asDict:
412                 __Parameters["withmpEnabled"]  = asDict["EnableMultiProcessing"]
413             if "NumberOfProcesses" in asDict:
414                 __Parameters["withmpWorkers"]  = asDict["NumberOfProcesses"]
415         #
416         if asScript is not None:
417             __Matrix, __Function = None, None
418             if asMatrix:
419                 __Matrix = ImportFromScript(asScript).getvalue( self.__name )
420             elif asOneFunction:
421                 __Function = { "Direct":ImportFromScript(asScript).getvalue( "DirectOperator" ) }
422                 __Function.update({"useApproximatedDerivatives":True})
423                 __Function.update(__Parameters)
424             elif asThreeFunctions:
425                 __Function = {
426                     "Direct" :ImportFromScript(asScript).getvalue( "DirectOperator" ),
427                     "Tangent":ImportFromScript(asScript).getvalue( "TangentOperator" ),
428                     "Adjoint":ImportFromScript(asScript).getvalue( "AdjointOperator" ),
429                     }
430                 __Function.update(__Parameters)
431         else:
432             __Matrix = asMatrix
433             if asOneFunction is not None:
434                 if isinstance(asOneFunction, dict) and "Direct" in asOneFunction:
435                     if asOneFunction["Direct"] is not None:
436                         __Function = asOneFunction
437                     else:
438                         raise ValueError("The function has to be given in a dictionnary which have 1 key (\"Direct\")")
439                 else:
440                     __Function = { "Direct":asOneFunction }
441                 __Function.update({"useApproximatedDerivatives":True})
442                 __Function.update(__Parameters)
443             elif asThreeFunctions is not None:
444                 if isinstance(asThreeFunctions, dict) and \
445                    ("Tangent" in asThreeFunctions) and (asThreeFunctions["Tangent"] is not None) and \
446                    ("Adjoint" in asThreeFunctions) and (asThreeFunctions["Adjoint"] is not None) and \
447                    (("useApproximatedDerivatives" not in asThreeFunctions) or not bool(asThreeFunctions["useApproximatedDerivatives"])):
448                     __Function = asThreeFunctions
449                 elif isinstance(asThreeFunctions, dict) and \
450                    ("Direct" in asThreeFunctions) and (asThreeFunctions["Direct"] is not None):
451                     __Function = asThreeFunctions
452                     __Function.update({"useApproximatedDerivatives":True})
453                 else:
454                     raise ValueError("The functions has to be given in a dictionnary which have either 1 key (\"Direct\") or 3 keys (\"Direct\" (optionnal), \"Tangent\" and \"Adjoint\")")
455                 if "Direct"  not in asThreeFunctions:
456                     __Function["Direct"] = asThreeFunctions["Tangent"]
457                 __Function.update(__Parameters)
458             else:
459                 __Function = None
460         #
461         # if sys.version_info[0] < 3 and isinstance(__Function, dict):
462         #     for k in ("Direct", "Tangent", "Adjoint"):
463         #         if k in __Function and hasattr(__Function[k],"__class__"):
464         #             if type(__Function[k]) is type(self.__init__):
465         #                 raise TypeError("can't use a class method (%s) as a function for the \"%s\" operator. Use a real function instead."%(type(__Function[k]),k))
466         #
467         if   appliedInX is not None and isinstance(appliedInX, dict):
468             __appliedInX = appliedInX
469         elif appliedInX is not None:
470             __appliedInX = {"HXb":appliedInX}
471         else:
472             __appliedInX = None
473         #
474         if scheduledBy is not None:
475             self.__T = scheduledBy
476         #
477         if isinstance(__Function, dict) and \
478                 ("useApproximatedDerivatives" in __Function) and bool(__Function["useApproximatedDerivatives"]) and \
479                 ("Direct" in __Function) and (__Function["Direct"] is not None):
480             if "withCenteredDF"            not in __Function: __Function["withCenteredDF"]            = False
481             if "withIncrement"             not in __Function: __Function["withIncrement"]             = 0.01
482             if "withdX"                    not in __Function: __Function["withdX"]                    = None
483             if "withAvoidingRedundancy"    not in __Function: __Function["withAvoidingRedundancy"]    = avoidRC
484             if "withToleranceInRedundancy" not in __Function: __Function["withToleranceInRedundancy"] = 1.e-18
485             if "withLenghtOfRedundancy"    not in __Function: __Function["withLenghtOfRedundancy"]    = -1
486             if "withmpEnabled"             not in __Function: __Function["withmpEnabled"]             = False
487             if "withmpWorkers"             not in __Function: __Function["withmpWorkers"]             = None
488             if "withmfEnabled"             not in __Function: __Function["withmfEnabled"]             = inputAsMF
489             from daNumerics.ApproximatedDerivatives import FDApproximation
490             FDA = FDApproximation(
491                 Function              = __Function["Direct"],
492                 centeredDF            = __Function["withCenteredDF"],
493                 increment             = __Function["withIncrement"],
494                 dX                    = __Function["withdX"],
495                 avoidingRedundancy    = __Function["withAvoidingRedundancy"],
496                 toleranceInRedundancy = __Function["withToleranceInRedundancy"],
497                 lenghtOfRedundancy    = __Function["withLenghtOfRedundancy"],
498                 mpEnabled             = __Function["withmpEnabled"],
499                 mpWorkers             = __Function["withmpWorkers"],
500                 mfEnabled             = __Function["withmfEnabled"],
501                 )
502             self.__FO["Direct"]  = Operator( fromMethod = FDA.DirectOperator,  avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
503             self.__FO["Tangent"] = Operator( fromMethod = FDA.TangentOperator, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
504             self.__FO["Adjoint"] = Operator( fromMethod = FDA.AdjointOperator, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
505         elif isinstance(__Function, dict) and \
506                 ("Direct" in __Function) and ("Tangent" in __Function) and ("Adjoint" in __Function) and \
507                 (__Function["Direct"] is not None) and (__Function["Tangent"] is not None) and (__Function["Adjoint"] is not None):
508             self.__FO["Direct"]  = Operator( fromMethod = __Function["Direct"],  avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
509             self.__FO["Tangent"] = Operator( fromMethod = __Function["Tangent"], avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
510             self.__FO["Adjoint"] = Operator( fromMethod = __Function["Adjoint"], avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF, extraArguments = self.__extraArgs )
511         elif asMatrix is not None:
512             __matrice = numpy.matrix( __Matrix, numpy.float )
513             self.__FO["Direct"]  = Operator( fromMatrix = __matrice,   avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF )
514             self.__FO["Tangent"] = Operator( fromMatrix = __matrice,   avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF )
515             self.__FO["Adjoint"] = Operator( fromMatrix = __matrice.T, avoidingRedundancy = avoidRC, inputAsMultiFunction = inputAsMF )
516             del __matrice
517         else:
518             raise ValueError("Improperly defined observation operator, it requires at minima either a matrix, a Direct for approximate derivatives or a Tangent/Adjoint pair.")
519         #
520         if __appliedInX is not None:
521             self.__FO["AppliedInX"] = {}
522             for key in list(__appliedInX.keys()):
523                 if type( __appliedInX[key] ) is type( numpy.matrix([]) ):
524                     # Pour le cas où l'on a une vraie matrice
525                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key].A1, numpy.float ).T
526                 elif type( __appliedInX[key] ) is type( numpy.array([]) ) and len(__appliedInX[key].shape) > 1:
527                     # Pour le cas où l'on a un vecteur représenté en array avec 2 dimensions
528                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key].reshape(len(__appliedInX[key]),), numpy.float ).T
529                 else:
530                     self.__FO["AppliedInX"][key] = numpy.matrix( __appliedInX[key],    numpy.float ).T
531         else:
532             self.__FO["AppliedInX"] = None
533
534     def getO(self):
535         return self.__FO
536
537     def __repr__(self):
538         "x.__repr__() <==> repr(x)"
539         return repr(self.__V)
540
541     def __str__(self):
542         "x.__str__() <==> str(x)"
543         return str(self.__V)
544
545 # ==============================================================================
546 class Algorithm(object):
547     """
548     Classe générale d'interface de type algorithme
549
550     Elle donne un cadre pour l'écriture d'une classe élémentaire d'algorithme
551     d'assimilation, en fournissant un container (dictionnaire) de variables
552     persistantes initialisées, et des méthodes d'accès à ces variables stockées.
553
554     Une classe élémentaire d'algorithme doit implémenter la méthode "run".
555     """
556     def __init__(self, name):
557         """
558         L'initialisation présente permet de fabriquer des variables de stockage
559         disponibles de manière générique dans les algorithmes élémentaires. Ces
560         variables de stockage sont ensuite conservées dans un dictionnaire
561         interne à l'objet, mais auquel on accède par la méthode "get".
562
563         Les variables prévues sont :
564             - APosterioriCorrelations : matrice de corrélations de la matrice A
565             - APosterioriCovariance : matrice de covariances a posteriori : A
566             - APosterioriStandardDeviations : vecteur des écart-types de la matrice A
567             - APosterioriVariances : vecteur des variances de la matrice A
568             - Analysis : vecteur d'analyse : Xa
569             - BMA : Background moins Analysis : Xa - Xb
570             - CostFunctionJ  : fonction-coût globale, somme des deux parties suivantes Jb et Jo
571             - CostFunctionJAtCurrentOptimum : fonction-coût globale à l'état optimal courant lors d'itérations
572             - CostFunctionJb : partie ébauche ou background de la fonction-coût : Jb
573             - CostFunctionJbAtCurrentOptimum : partie ébauche à l'état optimal courant lors d'itérations
574             - CostFunctionJo : partie observations de la fonction-coût : Jo
575             - CostFunctionJoAtCurrentOptimum : partie observations à l'état optimal courant lors d'itérations
576             - CurrentOptimum : état optimal courant lors d'itérations
577             - CurrentState : état courant lors d'itérations
578             - GradientOfCostFunctionJ  : gradient de la fonction-coût globale
579             - GradientOfCostFunctionJb : gradient de la partie ébauche de la fonction-coût
580             - GradientOfCostFunctionJo : gradient de la partie observations de la fonction-coût
581             - IndexOfOptimum : index de l'état optimal courant lors d'itérations
582             - Innovation : l'innovation : d = Y - H(X)
583             - InnovationAtCurrentState : l'innovation à l'état courant : dn = Y - H(Xn)
584             - JacobianMatrixAtCurrentState : matrice jacobienne à l'état courant
585             - JacobianMatrixAtOptimum : matrice jacobienne à l'optimum
586             - KalmanGainAtOptimum : gain de Kalman à l'optimum
587             - MahalanobisConsistency : indicateur de consistance des covariances
588             - OMA : Observation moins Analyse : Y - Xa
589             - OMB : Observation moins Background : Y - Xb
590             - PredictedState : état prédit courant lors d'itérations
591             - Residu : dans le cas des algorithmes de vérification
592             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
593             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
594             - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
595             - SimulatedObservationAtCurrentOptimum : l'état observé H(X) à l'état optimal courant
596             - SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
597             - SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
598             - SimulationQuantiles : états observés H(X) pour les quantiles demandés
599         On peut rajouter des variables à stocker dans l'initialisation de
600         l'algorithme élémentaire qui va hériter de cette classe
601         """
602         logging.debug("%s Initialisation", str(name))
603         self._m = PlatformInfo.SystemUsage()
604         #
605         self._name = str( name )
606         self._parameters = {"StoreSupplementaryCalculations":[]}
607         self.__required_parameters = {}
608         self.__required_inputs = {"RequiredInputValues":{"mandatory":(), "optional":()}}
609         #
610         self.StoredVariables = {}
611         self.StoredVariables["APosterioriCorrelations"]              = Persistence.OneMatrix(name = "APosterioriCorrelations")
612         self.StoredVariables["APosterioriCovariance"]                = Persistence.OneMatrix(name = "APosterioriCovariance")
613         self.StoredVariables["APosterioriStandardDeviations"]        = Persistence.OneVector(name = "APosterioriStandardDeviations")
614         self.StoredVariables["APosterioriVariances"]                 = Persistence.OneVector(name = "APosterioriVariances")
615         self.StoredVariables["Analysis"]                             = Persistence.OneVector(name = "Analysis")
616         self.StoredVariables["BMA"]                                  = Persistence.OneVector(name = "BMA")
617         self.StoredVariables["CostFunctionJ"]                        = Persistence.OneScalar(name = "CostFunctionJ")
618         self.StoredVariables["CostFunctionJAtCurrentOptimum"]        = Persistence.OneScalar(name = "CostFunctionJAtCurrentOptimum")
619         self.StoredVariables["CostFunctionJb"]                       = Persistence.OneScalar(name = "CostFunctionJb")
620         self.StoredVariables["CostFunctionJbAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
621         self.StoredVariables["CostFunctionJo"]                       = Persistence.OneScalar(name = "CostFunctionJo")
622         self.StoredVariables["CostFunctionJoAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
623         self.StoredVariables["CurrentOptimum"]                       = Persistence.OneVector(name = "CurrentOptimum")
624         self.StoredVariables["CurrentState"]                         = Persistence.OneVector(name = "CurrentState")
625         self.StoredVariables["GradientOfCostFunctionJ"]              = Persistence.OneVector(name = "GradientOfCostFunctionJ")
626         self.StoredVariables["GradientOfCostFunctionJb"]             = Persistence.OneVector(name = "GradientOfCostFunctionJb")
627         self.StoredVariables["GradientOfCostFunctionJo"]             = Persistence.OneVector(name = "GradientOfCostFunctionJo")
628         self.StoredVariables["IndexOfOptimum"]                       = Persistence.OneIndex(name = "IndexOfOptimum")
629         self.StoredVariables["Innovation"]                           = Persistence.OneVector(name = "Innovation")
630         self.StoredVariables["InnovationAtCurrentState"]             = Persistence.OneVector(name = "InnovationAtCurrentState")
631         self.StoredVariables["JacobianMatrixAtCurrentState"]         = Persistence.OneMatrix(name = "JacobianMatrixAtCurrentState")
632         self.StoredVariables["JacobianMatrixAtOptimum"]              = Persistence.OneMatrix(name = "JacobianMatrixAtOptimum")
633         self.StoredVariables["KalmanGainAtOptimum"]                  = Persistence.OneMatrix(name = "KalmanGainAtOptimum")
634         self.StoredVariables["MahalanobisConsistency"]               = Persistence.OneScalar(name = "MahalanobisConsistency")
635         self.StoredVariables["OMA"]                                  = Persistence.OneVector(name = "OMA")
636         self.StoredVariables["OMB"]                                  = Persistence.OneVector(name = "OMB")
637         self.StoredVariables["PredictedState"]                       = Persistence.OneVector(name = "PredictedState")
638         self.StoredVariables["Residu"]                               = Persistence.OneScalar(name = "Residu")
639         self.StoredVariables["SigmaBck2"]                            = Persistence.OneScalar(name = "SigmaBck2")
640         self.StoredVariables["SigmaObs2"]                            = Persistence.OneScalar(name = "SigmaObs2")
641         self.StoredVariables["SimulatedObservationAtBackground"]     = Persistence.OneVector(name = "SimulatedObservationAtBackground")
642         self.StoredVariables["SimulatedObservationAtCurrentOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentOptimum")
643         self.StoredVariables["SimulatedObservationAtCurrentState"]   = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
644         self.StoredVariables["SimulatedObservationAtOptimum"]        = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
645         self.StoredVariables["SimulationQuantiles"]                  = Persistence.OneMatrix(name = "SimulationQuantiles")
646
647     def _pre_run(self, Parameters, Xb=None, Y=None, R=None, B=None, Q=None ):
648         "Pré-calcul"
649         logging.debug("%s Lancement", self._name)
650         logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
651         #
652         # Mise a jour de self._parameters avec Parameters
653         self.__setParameters(Parameters)
654         #
655         # Corrections et complements
656         def __test_vvalue(argument, variable, argname):
657             if argument is None:
658                 if variable in self.__required_inputs["RequiredInputValues"]["mandatory"]:
659                     raise ValueError("%s %s vector %s has to be properly defined!"%(self._name,argname,variable))
660                 elif variable in self.__required_inputs["RequiredInputValues"]["optional"]:
661                     logging.debug("%s %s vector %s is not set, but is optional."%(self._name,argname,variable))
662                 else:
663                     logging.debug("%s %s vector %s is not set, but is not required."%(self._name,argname,variable))
664             else:
665                 logging.debug("%s %s vector %s is set, and its size is %i."%(self._name,argname,variable,numpy.array(argument).size))
666             return 0
667         __test_vvalue( Xb, "Xb", "Background or initial state" )
668         __test_vvalue( Y,  "Y",  "Observation" )
669         #
670         def __test_cvalue(argument, variable, argname):
671             if argument is None:
672                 if variable in self.__required_inputs["RequiredInputValues"]["mandatory"]:
673                     raise ValueError("%s %s error covariance matrix %s has to be properly defined!"%(self._name,argname,variable))
674                 elif variable in self.__required_inputs["RequiredInputValues"]["optional"]:
675                     logging.debug("%s %s error covariance matrix %s is not set, but is optional."%(self._name,argname,variable))
676                 else:
677                     logging.debug("%s %s error covariance matrix %s is not set, but is not required."%(self._name,argname,variable))
678             else:
679                 logging.debug("%s %s error covariance matrix %s is set."%(self._name,argname,variable))
680             return 0
681         __test_cvalue( R, "R", "Observation" )
682         __test_cvalue( B, "B", "Background" )
683         __test_cvalue( Q, "Q", "Evolution" )
684         #
685         if ("Bounds" in self._parameters) and isinstance(self._parameters["Bounds"], (list, tuple)) and (len(self._parameters["Bounds"]) > 0):
686             logging.debug("%s Prise en compte des bornes effectuee"%(self._name,))
687         else:
688             self._parameters["Bounds"] = None
689         #
690         if logging.getLogger().level < logging.WARNING:
691             self._parameters["optiprint"], self._parameters["optdisp"] = 1, 1
692             if PlatformInfo.has_scipy:
693                 import scipy.optimize
694                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_ALL
695             else:
696                 self._parameters["optmessages"] = 15
697         else:
698             self._parameters["optiprint"], self._parameters["optdisp"] = -1, 0
699             if PlatformInfo.has_scipy:
700                 import scipy.optimize
701                 self._parameters["optmessages"] = scipy.optimize.tnc.MSG_NONE
702             else:
703                 self._parameters["optmessages"] = 15
704         #
705         return 0
706
707     def _post_run(self,_oH=None):
708         "Post-calcul"
709         if ("StoreSupplementaryCalculations" in self._parameters) and \
710             "APosterioriCovariance" in self._parameters["StoreSupplementaryCalculations"]:
711             for _A in self.StoredVariables["APosterioriCovariance"]:
712                 if "APosterioriVariances" in self._parameters["StoreSupplementaryCalculations"]:
713                     self.StoredVariables["APosterioriVariances"].store( numpy.diag(_A) )
714                 if "APosterioriStandardDeviations" in self._parameters["StoreSupplementaryCalculations"]:
715                     self.StoredVariables["APosterioriStandardDeviations"].store( numpy.sqrt(numpy.diag(_A)) )
716                 if "APosterioriCorrelations" in self._parameters["StoreSupplementaryCalculations"]:
717                     _EI = numpy.diag(1./numpy.sqrt(numpy.diag(_A)))
718                     _C = numpy.dot(_EI, numpy.dot(_A, _EI))
719                     self.StoredVariables["APosterioriCorrelations"].store( _C )
720         if _oH is not None:
721             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))
722             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))
723         logging.debug("%s Taille mémoire utilisée de %.0f Mio", self._name, self._m.getUsedMemory("Mio"))
724         logging.debug("%s Terminé", self._name)
725         return 0
726
727     def _toStore(self, key):
728         "True if in StoreSupplementaryCalculations, else False"
729         return key in self._parameters["StoreSupplementaryCalculations"]
730
731     def get(self, key=None):
732         """
733         Renvoie l'une des variables stockées identifiée par la clé, ou le
734         dictionnaire de l'ensemble des variables disponibles en l'absence de
735         clé. Ce sont directement les variables sous forme objet qui sont
736         renvoyées, donc les méthodes d'accès à l'objet individuel sont celles
737         des classes de persistance.
738         """
739         if key is not None:
740             return self.StoredVariables[key]
741         else:
742             return self.StoredVariables
743
744     def __contains__(self, key=None):
745         "D.__contains__(k) -> True if D has a key k, else False"
746         return key in self.StoredVariables
747
748     def keys(self):
749         "D.keys() -> list of D's keys"
750         if hasattr(self, "StoredVariables"):
751             return self.StoredVariables.keys()
752         else:
753             return []
754
755     def pop(self, k, d):
756         "D.pop(k[,d]) -> v, remove specified key and return the corresponding value"
757         if hasattr(self, "StoredVariables"):
758             return self.StoredVariables.pop(k, d)
759         else:
760             try:
761                 msg = "'%s'"%k
762             except:
763                 raise TypeError("pop expected at least 1 arguments, got 0")
764             "If key is not found, d is returned if given, otherwise KeyError is raised"
765             try:
766                 return d
767             except:
768                 raise KeyError(msg)
769
770     def run(self, Xb=None, Y=None, H=None, M=None, R=None, B=None, Q=None, Parameters=None):
771         """
772         Doit implémenter l'opération élémentaire de calcul d'assimilation sous
773         sa forme mathématique la plus naturelle possible.
774         """
775         raise NotImplementedError("Mathematical assimilation calculation has not been implemented!")
776
777     def defineRequiredParameter(self, name = None, default = None, typecast = None, message = None, minval = None, maxval = None, listval = None):
778         """
779         Permet de définir dans l'algorithme des paramètres requis et leurs
780         caractéristiques par défaut.
781         """
782         if name is None:
783             raise ValueError("A name is mandatory to define a required parameter.")
784         #
785         self.__required_parameters[name] = {
786             "default"  : default,
787             "typecast" : typecast,
788             "minval"   : minval,
789             "maxval"   : maxval,
790             "listval"  : listval,
791             "message"  : message,
792             }
793         logging.debug("%s %s (valeur par défaut = %s)", self._name, message, self.setParameterValue(name))
794
795     def getRequiredParameters(self, noDetails=True):
796         """
797         Renvoie la liste des noms de paramètres requis ou directement le
798         dictionnaire des paramètres requis.
799         """
800         if noDetails:
801             return sorted(self.__required_parameters.keys())
802         else:
803             return self.__required_parameters
804
805     def setParameterValue(self, name=None, value=None):
806         """
807         Renvoie la valeur d'un paramètre requis de manière contrôlée
808         """
809         default  = self.__required_parameters[name]["default"]
810         typecast = self.__required_parameters[name]["typecast"]
811         minval   = self.__required_parameters[name]["minval"]
812         maxval   = self.__required_parameters[name]["maxval"]
813         listval  = self.__required_parameters[name]["listval"]
814         #
815         if value is None and default is None:
816             __val = None
817         elif value is None and default is not None:
818             if typecast is None: __val = default
819             else:                __val = typecast( default )
820         else:
821             if typecast is None: __val = value
822             else:                __val = typecast( value )
823         #
824         if minval is not None and (numpy.array(__val, float) < minval).any():
825             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be less than %s."%(name, __val, minval))
826         if maxval is not None and (numpy.array(__val, float) > maxval).any():
827             raise ValueError("The parameter named \"%s\" of value \"%s\" can not be greater than %s."%(name, __val, maxval))
828         if listval is not None:
829             if typecast is list or typecast is tuple or isinstance(__val,list) or isinstance(__val,tuple):
830                 for v in __val:
831                     if v not in listval:
832                         raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%(v, name, listval))
833             elif __val not in listval:
834                 raise ValueError("The value \"%s\" of the parameter named \"%s\" is not allowed, it has to be in the list %s."%( __val, name,listval))
835         return __val
836
837     def requireInputArguments(self, mandatory=(), optional=()):
838         """
839         Permet d'imposer des arguments requises en entrée
840         """
841         self.__required_inputs["RequiredInputValues"]["mandatory"] = tuple( mandatory )
842         self.__required_inputs["RequiredInputValues"]["optional"]  = tuple( optional )
843
844     def __setParameters(self, fromDico={}):
845         """
846         Permet de stocker les paramètres reçus dans le dictionnaire interne.
847         """
848         self._parameters.update( fromDico )
849         for k in self.__required_parameters.keys():
850             if k in fromDico.keys():
851                 self._parameters[k] = self.setParameterValue(k,fromDico[k])
852             else:
853                 self._parameters[k] = self.setParameterValue(k)
854             logging.debug("%s %s : %s", self._name, self.__required_parameters[k]["message"], self._parameters[k])
855
856 # ==============================================================================
857 class AlgorithmAndParameters(object):
858     """
859     Classe générale d'interface d'action pour l'algorithme et ses paramètres
860     """
861     def __init__(self,
862                  name               = "GenericAlgorithm",
863                  asAlgorithm        = None,
864                  asDict             = None,
865                  asScript           = None,
866                 ):
867         """
868         """
869         self.__name       = str(name)
870         self.__A          = None
871         self.__P          = {}
872         #
873         self.__algorithm         = {}
874         self.__algorithmFile     = None
875         self.__algorithmName     = None
876         #
877         self.updateParameters( asDict, asScript )
878         #
879         if asAlgorithm is None and asScript is not None:
880             __Algo = ImportFromScript(asScript).getvalue( "Algorithm" )
881         else:
882             __Algo = asAlgorithm
883         #
884         if __Algo is not None:
885             self.__A = str(__Algo)
886             self.__P.update( {"Algorithm":self.__A} )
887         #
888         self.__setAlgorithm( self.__A )
889
890     def updateParameters(self,
891                  asDict     = None,
892                  asScript   = None,
893                 ):
894         "Mise a jour des parametres"
895         if asDict is None and asScript is not None:
896             __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
897         else:
898             __Dict = asDict
899         #
900         if __Dict is not None:
901             self.__P.update( dict(__Dict) )
902
903     def executePythonScheme(self, asDictAO = None):
904         "Permet de lancer le calcul d'assimilation"
905         Operator.CM.clearCache()
906         #
907         if not isinstance(asDictAO, dict):
908             raise ValueError("The objects for algorithm calculation have to be given together as a dictionnary, and they are not")
909         if   hasattr(asDictAO["Background"],"getO"):        self.__Xb = asDictAO["Background"].getO()
910         elif hasattr(asDictAO["CheckingPoint"],"getO"):     self.__Xb = asDictAO["CheckingPoint"].getO()
911         else:                                               self.__Xb = None
912         if hasattr(asDictAO["Observation"],"getO"):         self.__Y  = asDictAO["Observation"].getO()
913         else:                                               self.__Y  = asDictAO["Observation"]
914         if hasattr(asDictAO["ControlInput"],"getO"):        self.__U  = asDictAO["ControlInput"].getO()
915         else:                                               self.__U  = asDictAO["ControlInput"]
916         if hasattr(asDictAO["ObservationOperator"],"getO"): self.__HO = asDictAO["ObservationOperator"].getO()
917         else:                                               self.__HO = asDictAO["ObservationOperator"]
918         if hasattr(asDictAO["EvolutionModel"],"getO"):      self.__EM = asDictAO["EvolutionModel"].getO()
919         else:                                               self.__EM = asDictAO["EvolutionModel"]
920         if hasattr(asDictAO["ControlModel"],"getO"):        self.__CM = asDictAO["ControlModel"].getO()
921         else:                                               self.__CM = asDictAO["ControlModel"]
922         self.__B = asDictAO["BackgroundError"]
923         self.__R = asDictAO["ObservationError"]
924         self.__Q = asDictAO["EvolutionError"]
925         #
926         self.__shape_validate()
927         #
928         self.__algorithm.run(
929             Xb         = self.__Xb,
930             Y          = self.__Y,
931             U          = self.__U,
932             HO         = self.__HO,
933             EM         = self.__EM,
934             CM         = self.__CM,
935             R          = self.__R,
936             B          = self.__B,
937             Q          = self.__Q,
938             Parameters = self.__P,
939             )
940         return 0
941
942     def executeYACSScheme(self, FileName=None):
943         "Permet de lancer le calcul d'assimilation"
944         if FileName is None or not os.path.exists(FileName):
945             raise ValueError("a YACS file name has to be given for YACS execution.\n")
946         else:
947             __file    = os.path.abspath(FileName)
948             logging.debug("The YACS file name is \"%s\"."%__file)
949         if not PlatformInfo.has_salome or \
950             not PlatformInfo.has_yacs or \
951             not PlatformInfo.has_adao:
952             raise ImportError("\n\n"+\
953                 "Unable to get SALOME, YACS or ADAO environnement variables.\n"+\
954                 "Please load the right environnement before trying to use it.\n")
955         #
956         import pilot
957         import SALOMERuntime
958         import loader
959         SALOMERuntime.RuntimeSALOME_setRuntime()
960
961         r = pilot.getRuntime()
962         xmlLoader = loader.YACSLoader()
963         xmlLoader.registerProcCataLoader()
964         try:
965             catalogAd = r.loadCatalog("proc", __file)
966             r.addCatalog(catalogAd)
967         except:
968             pass
969
970         try:
971             p = xmlLoader.load(__file)
972         except IOError as ex:
973             print("The YACS XML schema file can not be loaded: %s"%(ex,))
974
975         logger = p.getLogger("parser")
976         if not logger.isEmpty():
977             print("The imported YACS XML schema has errors on parsing:")
978             print(logger.getStr())
979
980         if not p.isValid():
981             print("The YACS XML schema is not valid and will not be executed:")
982             print(p.getErrorReport())
983
984         info=pilot.LinkInfo(pilot.LinkInfo.ALL_DONT_STOP)
985         p.checkConsistency(info)
986         if info.areWarningsOrErrors():
987             print("The YACS XML schema is not coherent and will not be executed:")
988             print(info.getGlobalRepr())
989
990         e = pilot.ExecutorSwig()
991         e.RunW(p)
992         if p.getEffectiveState() != pilot.DONE:
993             print(p.getErrorReport())
994         #
995         return 0
996
997     def get(self, key = None):
998         "Vérifie l'existence d'une clé de variable ou de paramètres"
999         if key in self.__algorithm:
1000             return self.__algorithm.get( key )
1001         elif key in self.__P:
1002             return self.__P[key]
1003         else:
1004             return self.__P
1005
1006     def pop(self, k, d):
1007         "Necessaire pour le pickling"
1008         return self.__algorithm.pop(k, d)
1009
1010     def getAlgorithmRequiredParameters(self, noDetails=True):
1011         "Renvoie la liste des paramètres requis selon l'algorithme"
1012         return self.__algorithm.getRequiredParameters(noDetails)
1013
1014     def setObserver(self, __V, __O, __I, __S):
1015         if self.__algorithm is None \
1016             or isinstance(self.__algorithm, dict) \
1017             or not hasattr(self.__algorithm,"StoredVariables"):
1018             raise ValueError("No observer can be build before choosing an algorithm.")
1019         if __V not in self.__algorithm:
1020             raise ValueError("An observer requires to be set on a variable named %s which does not exist."%__V)
1021         else:
1022             self.__algorithm.StoredVariables[ __V ].setDataObserver(
1023                     Scheduler      = __S,
1024                     HookFunction   = __O,
1025                     HookParameters = __I,
1026                     )
1027
1028     def removeObserver(self, __V, __O, __A = False):
1029         if self.__algorithm is None \
1030             or isinstance(self.__algorithm, dict) \
1031             or not hasattr(self.__algorithm,"StoredVariables"):
1032             raise ValueError("No observer can be removed before choosing an algorithm.")
1033         if __V not in self.__algorithm:
1034             raise ValueError("An observer requires to be removed on a variable named %s which does not exist."%__V)
1035         else:
1036             return self.__algorithm.StoredVariables[ __V ].removeDataObserver(
1037                     HookFunction   = __O,
1038                     AllObservers   = __A,
1039                     )
1040
1041     def hasObserver(self, __V):
1042         if self.__algorithm is None \
1043             or isinstance(self.__algorithm, dict) \
1044             or not hasattr(self.__algorithm,"StoredVariables"):
1045             return False
1046         if __V not in self.__algorithm:
1047             return False
1048         return self.__algorithm.StoredVariables[ __V ].hasDataObserver()
1049
1050     def keys(self):
1051         return list(self.__algorithm.keys()) + list(self.__P.keys())
1052
1053     def __contains__(self, key=None):
1054         "D.__contains__(k) -> True if D has a key k, else False"
1055         return key in self.__algorithm or key in self.__P
1056
1057     def __repr__(self):
1058         "x.__repr__() <==> repr(x)"
1059         return repr(self.__A)+", "+repr(self.__P)
1060
1061     def __str__(self):
1062         "x.__str__() <==> str(x)"
1063         return str(self.__A)+", "+str(self.__P)
1064
1065     def __setAlgorithm(self, choice = None ):
1066         """
1067         Permet de sélectionner l'algorithme à utiliser pour mener à bien l'étude
1068         d'assimilation. L'argument est un champ caractère se rapportant au nom
1069         d'un algorithme réalisant l'opération sur les arguments fixes.
1070         """
1071         if choice is None:
1072             raise ValueError("Error: algorithm choice has to be given")
1073         if self.__algorithmName is not None:
1074             raise ValueError("Error: algorithm choice has already been done as \"%s\", it can't be changed."%self.__algorithmName)
1075         daDirectory = "daAlgorithms"
1076         #
1077         # Recherche explicitement le fichier complet
1078         # ------------------------------------------
1079         module_path = None
1080         for directory in sys.path:
1081             if os.path.isfile(os.path.join(directory, daDirectory, str(choice)+'.py')):
1082                 module_path = os.path.abspath(os.path.join(directory, daDirectory))
1083         if module_path is None:
1084             raise ImportError("No algorithm module named \"%s\" has been found in the search path.\n             The search path is %s"%(choice, sys.path))
1085         #
1086         # Importe le fichier complet comme un module
1087         # ------------------------------------------
1088         try:
1089             sys_path_tmp = sys.path ; sys.path.insert(0,module_path)
1090             self.__algorithmFile = __import__(str(choice), globals(), locals(), [])
1091             if not hasattr(self.__algorithmFile, "ElementaryAlgorithm"):
1092                 raise ImportError("this module does not define a valid elementary algorithm.")
1093             self.__algorithmName = str(choice)
1094             sys.path = sys_path_tmp ; del sys_path_tmp
1095         except ImportError as e:
1096             raise ImportError("The module named \"%s\" was found, but is incorrect at the import stage.\n             The import error message is: %s"%(choice,e))
1097         #
1098         # Instancie un objet du type élémentaire du fichier
1099         # -------------------------------------------------
1100         self.__algorithm = self.__algorithmFile.ElementaryAlgorithm()
1101         return 0
1102
1103     def __shape_validate(self):
1104         """
1105         Validation de la correspondance correcte des tailles des variables et
1106         des matrices s'il y en a.
1107         """
1108         if self.__Xb is None:                      __Xb_shape = (0,)
1109         elif hasattr(self.__Xb,"size"):            __Xb_shape = (self.__Xb.size,)
1110         elif hasattr(self.__Xb,"shape"):
1111             if isinstance(self.__Xb.shape, tuple): __Xb_shape = self.__Xb.shape
1112             else:                                  __Xb_shape = self.__Xb.shape()
1113         else: raise TypeError("The background (Xb) has no attribute of shape: problem !")
1114         #
1115         if self.__Y is None:                       __Y_shape = (0,)
1116         elif hasattr(self.__Y,"size"):             __Y_shape = (self.__Y.size,)
1117         elif hasattr(self.__Y,"shape"):
1118             if isinstance(self.__Y.shape, tuple):  __Y_shape = self.__Y.shape
1119             else:                                  __Y_shape = self.__Y.shape()
1120         else: raise TypeError("The observation (Y) has no attribute of shape: problem !")
1121         #
1122         if self.__U is None:                       __U_shape = (0,)
1123         elif hasattr(self.__U,"size"):             __U_shape = (self.__U.size,)
1124         elif hasattr(self.__U,"shape"):
1125             if isinstance(self.__U.shape, tuple):  __U_shape = self.__U.shape
1126             else:                                  __U_shape = self.__U.shape()
1127         else: raise TypeError("The control (U) has no attribute of shape: problem !")
1128         #
1129         if self.__B is None:                       __B_shape = (0,0)
1130         elif hasattr(self.__B,"shape"):
1131             if isinstance(self.__B.shape, tuple):  __B_shape = self.__B.shape
1132             else:                                  __B_shape = self.__B.shape()
1133         else: raise TypeError("The a priori errors covariance matrix (B) has no attribute of shape: problem !")
1134         #
1135         if self.__R is None:                       __R_shape = (0,0)
1136         elif hasattr(self.__R,"shape"):
1137             if isinstance(self.__R.shape, tuple):  __R_shape = self.__R.shape
1138             else:                                  __R_shape = self.__R.shape()
1139         else: raise TypeError("The observation errors covariance matrix (R) has no attribute of shape: problem !")
1140         #
1141         if self.__Q is None:                       __Q_shape = (0,0)
1142         elif hasattr(self.__Q,"shape"):
1143             if isinstance(self.__Q.shape, tuple):  __Q_shape = self.__Q.shape
1144             else:                                  __Q_shape = self.__Q.shape()
1145         else: raise TypeError("The evolution errors covariance matrix (Q) has no attribute of shape: problem !")
1146         #
1147         if len(self.__HO) == 0:                              __HO_shape = (0,0)
1148         elif isinstance(self.__HO, dict):                    __HO_shape = (0,0)
1149         elif hasattr(self.__HO["Direct"],"shape"):
1150             if isinstance(self.__HO["Direct"].shape, tuple): __HO_shape = self.__HO["Direct"].shape
1151             else:                                            __HO_shape = self.__HO["Direct"].shape()
1152         else: raise TypeError("The observation operator (H) has no attribute of shape: problem !")
1153         #
1154         if len(self.__EM) == 0:                              __EM_shape = (0,0)
1155         elif isinstance(self.__EM, dict):                    __EM_shape = (0,0)
1156         elif hasattr(self.__EM["Direct"],"shape"):
1157             if isinstance(self.__EM["Direct"].shape, tuple): __EM_shape = self.__EM["Direct"].shape
1158             else:                                            __EM_shape = self.__EM["Direct"].shape()
1159         else: raise TypeError("The evolution model (EM) has no attribute of shape: problem !")
1160         #
1161         if len(self.__CM) == 0:                              __CM_shape = (0,0)
1162         elif isinstance(self.__CM, dict):                    __CM_shape = (0,0)
1163         elif hasattr(self.__CM["Direct"],"shape"):
1164             if isinstance(self.__CM["Direct"].shape, tuple): __CM_shape = self.__CM["Direct"].shape
1165             else:                                            __CM_shape = self.__CM["Direct"].shape()
1166         else: raise TypeError("The control model (CM) has no attribute of shape: problem !")
1167         #
1168         # Vérification des conditions
1169         # ---------------------------
1170         if not( len(__Xb_shape) == 1 or min(__Xb_shape) == 1 ):
1171             raise ValueError("Shape characteristic of background (Xb) is incorrect: \"%s\"."%(__Xb_shape,))
1172         if not( len(__Y_shape) == 1 or min(__Y_shape) == 1 ):
1173             raise ValueError("Shape characteristic of observation (Y) is incorrect: \"%s\"."%(__Y_shape,))
1174         #
1175         if not( min(__B_shape) == max(__B_shape) ):
1176             raise ValueError("Shape characteristic of a priori errors covariance matrix (B) is incorrect: \"%s\"."%(__B_shape,))
1177         if not( min(__R_shape) == max(__R_shape) ):
1178             raise ValueError("Shape characteristic of observation errors covariance matrix (R) is incorrect: \"%s\"."%(__R_shape,))
1179         if not( min(__Q_shape) == max(__Q_shape) ):
1180             raise ValueError("Shape characteristic of evolution errors covariance matrix (Q) is incorrect: \"%s\"."%(__Q_shape,))
1181         if not( min(__EM_shape) == max(__EM_shape) ):
1182             raise ValueError("Shape characteristic of evolution operator (EM) is incorrect: \"%s\"."%(__EM_shape,))
1183         #
1184         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and not( __HO_shape[1] == max(__Xb_shape) ):
1185             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and state (X) \"%s\" are incompatible."%(__HO_shape,__Xb_shape))
1186         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and not( __HO_shape[0] == max(__Y_shape) ):
1187             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation (Y) \"%s\" are incompatible."%(__HO_shape,__Y_shape))
1188         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and len(self.__B) > 0 and not( __HO_shape[1] == __B_shape[0] ):
1189             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and a priori errors covariance matrix (B) \"%s\" are incompatible."%(__HO_shape,__B_shape))
1190         if len(self.__HO) > 0 and not isinstance(self.__HO, dict) and len(self.__R) > 0 and not( __HO_shape[0] == __R_shape[1] ):
1191             raise ValueError("Shape characteristic of observation operator (H) \"%s\" and observation errors covariance matrix (R) \"%s\" are incompatible."%(__HO_shape,__R_shape))
1192         #
1193         if self.__B is not None and len(self.__B) > 0 and not( __B_shape[1] == max(__Xb_shape) ):
1194             if self.__algorithmName in ["EnsembleBlue",]:
1195                 asPersistentVector = self.__Xb.reshape((-1,min(__B_shape)))
1196                 self.__Xb = Persistence.OneVector("Background", basetype=numpy.matrix)
1197                 for member in asPersistentVector:
1198                     self.__Xb.store( numpy.matrix( numpy.ravel(member), numpy.float ).T )
1199                 __Xb_shape = min(__B_shape)
1200             else:
1201                 raise ValueError("Shape characteristic of a priori errors covariance matrix (B) \"%s\" and background (Xb) \"%s\" are incompatible."%(__B_shape,__Xb_shape))
1202         #
1203         if self.__R is not None and len(self.__R) > 0 and not( __R_shape[1] == max(__Y_shape) ):
1204             raise ValueError("Shape characteristic of observation errors covariance matrix (R) \"%s\" and observation (Y) \"%s\" are incompatible."%(__R_shape,__Y_shape))
1205         #
1206         if self.__EM is not None and len(self.__EM) > 0 and not isinstance(self.__EM, dict) and not( __EM_shape[1] == max(__Xb_shape) ):
1207             raise ValueError("Shape characteristic of evolution model (EM) \"%s\" and state (X) \"%s\" are incompatible."%(__EM_shape,__Xb_shape))
1208         #
1209         if self.__CM is not None and len(self.__CM) > 0 and not isinstance(self.__CM, dict) and not( __CM_shape[1] == max(__U_shape) ):
1210             raise ValueError("Shape characteristic of control model (CM) \"%s\" and control (U) \"%s\" are incompatible."%(__CM_shape,__U_shape))
1211         #
1212         if ("Bounds" in self.__P) \
1213             and (isinstance(self.__P["Bounds"], list) or isinstance(self.__P["Bounds"], tuple)) \
1214             and (len(self.__P["Bounds"]) != max(__Xb_shape)):
1215             raise ValueError("The number \"%s\" of bound pairs for the state (X) components is different of the size \"%s\" of the state itself." \
1216                 %(len(self.__P["Bounds"]),max(__Xb_shape)))
1217         #
1218         return 1
1219
1220 # ==============================================================================
1221 class RegulationAndParameters(object):
1222     """
1223     Classe générale d'interface d'action pour la régulation et ses paramètres
1224     """
1225     def __init__(self,
1226                  name               = "GenericRegulation",
1227                  asAlgorithm        = None,
1228                  asDict             = None,
1229                  asScript           = None,
1230                 ):
1231         """
1232         """
1233         self.__name       = str(name)
1234         self.__P          = {}
1235         #
1236         if asAlgorithm is None and asScript is not None:
1237             __Algo = ImportFromScript(asScript).getvalue( "Algorithm" )
1238         else:
1239             __Algo = asAlgorithm
1240         #
1241         if asDict is None and asScript is not None:
1242             __Dict = ImportFromScript(asScript).getvalue( self.__name, "Parameters" )
1243         else:
1244             __Dict = asDict
1245         #
1246         if __Dict is not None:
1247             self.__P.update( dict(__Dict) )
1248         #
1249         if __Algo is not None:
1250             self.__P.update( {"Algorithm":self.__A} )
1251
1252     def get(self, key = None):
1253         "Vérifie l'existence d'une clé de variable ou de paramètres"
1254         if key in self.__P:
1255             return self.__P[key]
1256         else:
1257             return self.__P
1258
1259 # ==============================================================================
1260 class DataObserver(object):
1261     """
1262     Classe générale d'interface de type observer
1263     """
1264     def __init__(self,
1265                  name        = "GenericObserver",
1266                  onVariable  = None,
1267                  asTemplate  = None,
1268                  asString    = None,
1269                  asScript    = None,
1270                  asObsObject = None,
1271                  withInfo    = None,
1272                  scheduledBy = None,
1273                  withAlgo    = None,
1274                 ):
1275         """
1276         """
1277         self.__name       = str(name)
1278         self.__V          = None
1279         self.__O          = None
1280         self.__I          = None
1281         #
1282         if onVariable is None:
1283             raise ValueError("setting an observer has to be done over a variable name or a list of variable names, not over None.")
1284         elif type(onVariable) in (tuple, list):
1285             self.__V = tuple(map( str, onVariable ))
1286             if withInfo is None:
1287                 self.__I = self.__V
1288             else:
1289                 self.__I = (str(withInfo),)*len(self.__V)
1290         elif isinstance(onVariable, str):
1291             self.__V = (onVariable,)
1292             if withInfo is None:
1293                 self.__I = (onVariable,)
1294             else:
1295                 self.__I = (str(withInfo),)
1296         else:
1297             raise ValueError("setting an observer has to be done over a variable name or a list of variable names.")
1298         #
1299         if asString is not None:
1300             __FunctionText = asString
1301         elif (asTemplate is not None) and (asTemplate in Templates.ObserverTemplates):
1302             __FunctionText = Templates.ObserverTemplates[asTemplate]
1303         elif asScript is not None:
1304             __FunctionText = ImportFromScript(asScript).getstring()
1305         else:
1306             __FunctionText = ""
1307         __Function = ObserverF(__FunctionText)
1308         #
1309         if asObsObject is not None:
1310             self.__O = asObsObject
1311         else:
1312             self.__O = __Function.getfunc()
1313         #
1314         for k in range(len(self.__V)):
1315             ename = self.__V[k]
1316             einfo = self.__I[k]
1317             if ename not in withAlgo:
1318                 raise ValueError("An observer is asked to be set on a variable named %s which does not exist."%ename)
1319             else:
1320                 withAlgo.setObserver(ename, self.__O, einfo, scheduledBy)
1321
1322     def __repr__(self):
1323         "x.__repr__() <==> repr(x)"
1324         return repr(self.__V)+"\n"+repr(self.__O)
1325
1326     def __str__(self):
1327         "x.__str__() <==> str(x)"
1328         return str(self.__V)+"\n"+str(self.__O)
1329
1330 # ==============================================================================
1331 class State(object):
1332     """
1333     Classe générale d'interface de type état
1334     """
1335     def __init__(self,
1336                  name               = "GenericVector",
1337                  asVector           = None,
1338                  asPersistentVector = None,
1339                  asScript           = None,
1340                  asDataFile         = None,
1341                  colNames           = None,
1342                  colMajor           = False,
1343                  scheduledBy        = None,
1344                  toBeChecked        = False,
1345                 ):
1346         """
1347         Permet de définir un vecteur :
1348         - asVector : entrée des données, comme un vecteur compatible avec le
1349           constructeur de numpy.matrix, ou "True" si entrée par script.
1350         - asPersistentVector : entrée des données, comme une série de vecteurs
1351           compatible avec le constructeur de numpy.matrix, ou comme un objet de
1352           type Persistence, ou "True" si entrée par script.
1353         - asScript : si un script valide est donné contenant une variable
1354           nommée "name", la variable est de type "asVector" (par défaut) ou
1355           "asPersistentVector" selon que l'une de ces variables est placée à
1356           "True".
1357         - asDataFile : si un ou plusieurs fichiers valides sont donnés
1358           contenant des valeurs en colonnes, elles-mêmes nommées "colNames"
1359           (s'il n'y a pas de nom de colonne indiquée, on cherche une colonne
1360           nommée "name"), on récupère les colonnes et on les range ligne après
1361           ligne (colMajor=False) ou colonne après colonne (colMajor=True). La
1362           variable résultante est de type "asVector" (par défaut) ou
1363           "asPersistentVector" selon que l'une de ces variables est placée à
1364           "True".
1365         """
1366         self.__name       = str(name)
1367         self.__check      = bool(toBeChecked)
1368         #
1369         self.__V          = None
1370         self.__T          = None
1371         self.__is_vector  = False
1372         self.__is_series  = False
1373         #
1374         if asScript is not None:
1375             __Vector, __Series = None, None
1376             if asPersistentVector:
1377                 __Series = ImportFromScript(asScript).getvalue( self.__name )
1378             else:
1379                 __Vector = ImportFromScript(asScript).getvalue( self.__name )
1380         elif asDataFile is not None:
1381             __Vector, __Series = None, None
1382             if asPersistentVector:
1383                 if colNames is not None:
1384                     __Series = ImportFromFile(asDataFile).getvalue( colNames )[1]
1385                 else:
1386                     __Series = ImportFromFile(asDataFile).getvalue( [self.__name,] )[1]
1387                 if bool(colMajor) and not ImportFromFile(asDataFile).getformat() == "application/numpy.npz":
1388                     __Series = numpy.transpose(__Series)
1389                 elif not bool(colMajor) and ImportFromFile(asDataFile).getformat() == "application/numpy.npz":
1390                     __Series = numpy.transpose(__Series)
1391             else:
1392                 if colNames is not None:
1393                     __Vector = ImportFromFile(asDataFile).getvalue( colNames )[1]
1394                 else:
1395                     __Vector = ImportFromFile(asDataFile).getvalue( [self.__name,] )[1]
1396                 if bool(colMajor):
1397                     __Vector = numpy.ravel(__Vector, order = "F")
1398                 else:
1399                     __Vector = numpy.ravel(__Vector, order = "C")
1400         else:
1401             __Vector, __Series = asVector, asPersistentVector
1402         #
1403         if __Vector is not None:
1404             self.__is_vector = True
1405             self.__V         = numpy.matrix( numpy.asmatrix(__Vector).A1, numpy.float ).T
1406             self.shape       = self.__V.shape
1407             self.size        = self.__V.size
1408         elif __Series is not None:
1409             self.__is_series  = True
1410             if isinstance(__Series, (tuple, list, numpy.ndarray, numpy.matrix, str)):
1411                 self.__V = Persistence.OneVector(self.__name, basetype=numpy.matrix)
1412                 if isinstance(__Series, str): __Series = eval(__Series)
1413                 for member in __Series:
1414                     self.__V.store( numpy.matrix( numpy.asmatrix(member).A1, numpy.float ).T )
1415             else:
1416                 self.__V = __Series
1417             if isinstance(self.__V.shape, (tuple, list)):
1418                 self.shape       = self.__V.shape
1419             else:
1420                 self.shape       = self.__V.shape()
1421             if len(self.shape) == 1:
1422                 self.shape       = (self.shape[0],1)
1423             self.size        = self.shape[0] * self.shape[1]
1424         else:
1425             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)
1426         #
1427         if scheduledBy is not None:
1428             self.__T = scheduledBy
1429
1430     def getO(self, withScheduler=False):
1431         if withScheduler:
1432             return self.__V, self.__T
1433         elif self.__T is None:
1434             return self.__V
1435         else:
1436             return self.__V
1437
1438     def isvector(self):
1439         "Vérification du type interne"
1440         return self.__is_vector
1441
1442     def isseries(self):
1443         "Vérification du type interne"
1444         return self.__is_series
1445
1446     def __repr__(self):
1447         "x.__repr__() <==> repr(x)"
1448         return repr(self.__V)
1449
1450     def __str__(self):
1451         "x.__str__() <==> str(x)"
1452         return str(self.__V)
1453
1454 # ==============================================================================
1455 class Covariance(object):
1456     """
1457     Classe générale d'interface de type covariance
1458     """
1459     def __init__(self,
1460                  name          = "GenericCovariance",
1461                  asCovariance  = None,
1462                  asEyeByScalar = None,
1463                  asEyeByVector = None,
1464                  asCovObject   = None,
1465                  asScript      = None,
1466                  toBeChecked   = False,
1467                 ):
1468         """
1469         Permet de définir une covariance :
1470         - asCovariance : entrée des données, comme une matrice compatible avec
1471           le constructeur de numpy.matrix
1472         - asEyeByScalar : entrée des données comme un seul scalaire de variance,
1473           multiplicatif d'une matrice de corrélation identité, aucune matrice
1474           n'étant donc explicitement à donner
1475         - asEyeByVector : entrée des données comme un seul vecteur de variance,
1476           à mettre sur la diagonale d'une matrice de corrélation, aucune matrice
1477           n'étant donc explicitement à donner
1478         - asCovObject : entrée des données comme un objet python, qui a les
1479           methodes obligatoires "getT", "getI", "diag", "trace", "__add__",
1480           "__sub__", "__neg__", "__mul__", "__rmul__" et facultatives "shape",
1481           "size", "cholesky", "choleskyI", "asfullmatrix", "__repr__", "__str__"
1482         - toBeChecked : booléen indiquant si le caractère SDP de la matrice
1483           pleine doit être vérifié
1484         """
1485         self.__name       = str(name)
1486         self.__check      = bool(toBeChecked)
1487         #
1488         self.__C          = None
1489         self.__is_scalar  = False
1490         self.__is_vector  = False
1491         self.__is_matrix  = False
1492         self.__is_object  = False
1493         #
1494         if asScript is not None:
1495             __Matrix, __Scalar, __Vector, __Object = None, None, None, None
1496             if asEyeByScalar:
1497                 __Scalar = ImportFromScript(asScript).getvalue( self.__name )
1498             elif asEyeByVector:
1499                 __Vector = ImportFromScript(asScript).getvalue( self.__name )
1500             elif asCovObject:
1501                 __Object = ImportFromScript(asScript).getvalue( self.__name )
1502             else:
1503                 __Matrix = ImportFromScript(asScript).getvalue( self.__name )
1504         else:
1505             __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject
1506         #
1507         if __Scalar is not None:
1508             if numpy.matrix(__Scalar).size != 1:
1509                 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)
1510             self.__is_scalar = True
1511             self.__C         = numpy.abs( float(__Scalar) )
1512             self.shape       = (0,0)
1513             self.size        = 0
1514         elif __Vector is not None:
1515             self.__is_vector = True
1516             self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(__Vector, float ) ) ) )
1517             self.shape       = (self.__C.size,self.__C.size)
1518             self.size        = self.__C.size**2
1519         elif __Matrix is not None:
1520             self.__is_matrix = True
1521             self.__C         = numpy.matrix( __Matrix, float )
1522             self.shape       = self.__C.shape
1523             self.size        = self.__C.size
1524         elif __Object is not None:
1525             self.__is_object = True
1526             self.__C         = __Object
1527             for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
1528                 if not hasattr(self.__C,at):
1529                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
1530             if hasattr(self.__C,"shape"):
1531                 self.shape       = self.__C.shape
1532             else:
1533                 self.shape       = (0,0)
1534             if hasattr(self.__C,"size"):
1535                 self.size        = self.__C.size
1536             else:
1537                 self.size        = 0
1538         else:
1539             pass
1540             # 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)
1541         #
1542         self.__validate()
1543
1544     def __validate(self):
1545         "Validation"
1546         if self.__C is None:
1547             raise UnboundLocalError("%s covariance matrix value has not been set!"%(self.__name,))
1548         if self.ismatrix() and min(self.shape) != max(self.shape):
1549             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))
1550         if self.isobject() and min(self.shape) != max(self.shape):
1551             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))
1552         if self.isscalar() and self.__C <= 0:
1553             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your scalar input %s."%(self.__name,self.__C))
1554         if self.isvector() and (self.__C <= 0).any():
1555             raise ValueError("The \"%s\" covariance matrix is not positive-definite. Please check your vector input."%(self.__name,))
1556         if self.ismatrix() and (self.__check or logging.getLogger().level < logging.WARNING):
1557             try:
1558                 L = numpy.linalg.cholesky( self.__C )
1559             except:
1560                 raise ValueError("The %s covariance matrix is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
1561         if self.isobject() and (self.__check or logging.getLogger().level < logging.WARNING):
1562             try:
1563                 L = self.__C.cholesky()
1564             except:
1565                 raise ValueError("The %s covariance object is not symmetric positive-definite. Please check your matrix input."%(self.__name,))
1566
1567     def isscalar(self):
1568         "Vérification du type interne"
1569         return self.__is_scalar
1570
1571     def isvector(self):
1572         "Vérification du type interne"
1573         return self.__is_vector
1574
1575     def ismatrix(self):
1576         "Vérification du type interne"
1577         return self.__is_matrix
1578
1579     def isobject(self):
1580         "Vérification du type interne"
1581         return self.__is_object
1582
1583     def getI(self):
1584         "Inversion"
1585         if   self.ismatrix():
1586             return Covariance(self.__name+"I", asCovariance  = self.__C.I )
1587         elif self.isvector():
1588             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
1589         elif self.isscalar():
1590             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
1591         elif self.isobject():
1592             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
1593         else:
1594             return None # Indispensable
1595
1596     def getT(self):
1597         "Transposition"
1598         if   self.ismatrix():
1599             return Covariance(self.__name+"T", asCovariance  = self.__C.T )
1600         elif self.isvector():
1601             return Covariance(self.__name+"T", asEyeByVector = self.__C )
1602         elif self.isscalar():
1603             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
1604         elif self.isobject():
1605             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
1606
1607     def cholesky(self):
1608         "Décomposition de Cholesky"
1609         if   self.ismatrix():
1610             return Covariance(self.__name+"C", asCovariance  = numpy.linalg.cholesky(self.__C) )
1611         elif self.isvector():
1612             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
1613         elif self.isscalar():
1614             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
1615         elif self.isobject() and hasattr(self.__C,"cholesky"):
1616             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
1617
1618     def choleskyI(self):
1619         "Inversion de la décomposition de Cholesky"
1620         if   self.ismatrix():
1621             return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
1622         elif self.isvector():
1623             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
1624         elif self.isscalar():
1625             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
1626         elif self.isobject() and hasattr(self.__C,"choleskyI"):
1627             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
1628
1629     def diag(self, msize=None):
1630         "Diagonale de la matrice"
1631         if   self.ismatrix():
1632             return numpy.diag(self.__C)
1633         elif self.isvector():
1634             return self.__C
1635         elif self.isscalar():
1636             if msize is None:
1637                 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,))
1638             else:
1639                 return self.__C * numpy.ones(int(msize))
1640         elif self.isobject():
1641             return self.__C.diag()
1642
1643     def asfullmatrix(self, msize=None):
1644         "Matrice pleine"
1645         if   self.ismatrix():
1646             return self.__C
1647         elif self.isvector():
1648             return numpy.matrix( numpy.diag(self.__C), float )
1649         elif self.isscalar():
1650             if msize is None:
1651                 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,))
1652             else:
1653                 return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
1654         elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
1655             return self.__C.asfullmatrix()
1656
1657     def trace(self, msize=None):
1658         "Trace de la matrice"
1659         if   self.ismatrix():
1660             return numpy.trace(self.__C)
1661         elif self.isvector():
1662             return float(numpy.sum(self.__C))
1663         elif self.isscalar():
1664             if msize is None:
1665                 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,))
1666             else:
1667                 return self.__C * int(msize)
1668         elif self.isobject():
1669             return self.__C.trace()
1670
1671     def getO(self):
1672         return self
1673
1674     def __repr__(self):
1675         "x.__repr__() <==> repr(x)"
1676         return repr(self.__C)
1677
1678     def __str__(self):
1679         "x.__str__() <==> str(x)"
1680         return str(self.__C)
1681
1682     def __add__(self, other):
1683         "x.__add__(y) <==> x+y"
1684         if   self.ismatrix() or self.isobject():
1685             return self.__C + numpy.asmatrix(other)
1686         elif self.isvector() or self.isscalar():
1687             _A = numpy.asarray(other)
1688             _A.reshape(_A.size)[::_A.shape[1]+1] += self.__C
1689             return numpy.asmatrix(_A)
1690
1691     def __radd__(self, other):
1692         "x.__radd__(y) <==> y+x"
1693         raise NotImplementedError("%s covariance matrix __radd__ method not available for %s type!"%(self.__name,type(other)))
1694
1695     def __sub__(self, other):
1696         "x.__sub__(y) <==> x-y"
1697         if   self.ismatrix() or self.isobject():
1698             return self.__C - numpy.asmatrix(other)
1699         elif self.isvector() or self.isscalar():
1700             _A = numpy.asarray(other)
1701             _A.reshape(_A.size)[::_A.shape[1]+1] = self.__C - _A.reshape(_A.size)[::_A.shape[1]+1]
1702             return numpy.asmatrix(_A)
1703
1704     def __rsub__(self, other):
1705         "x.__rsub__(y) <==> y-x"
1706         raise NotImplementedError("%s covariance matrix __rsub__ method not available for %s type!"%(self.__name,type(other)))
1707
1708     def __neg__(self):
1709         "x.__neg__() <==> -x"
1710         return - self.__C
1711
1712     def __mul__(self, other):
1713         "x.__mul__(y) <==> x*y"
1714         if   self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
1715             return self.__C * other
1716         elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
1717             if numpy.ravel(other).size == self.shape[1]: # Vecteur
1718                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
1719             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
1720                 return self.__C * numpy.asmatrix(other)
1721             else:
1722                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asmatrix(other).shape,self.__name))
1723         elif self.isvector() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
1724             if numpy.ravel(other).size == self.shape[1]: # Vecteur
1725                 return numpy.asmatrix(self.__C * numpy.ravel(other)).T
1726             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
1727                 return numpy.asmatrix((self.__C * (numpy.asarray(other).transpose())).transpose())
1728             else:
1729                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
1730         elif self.isscalar() and isinstance(other,numpy.matrix):
1731             return self.__C * other
1732         elif self.isscalar() and isinstance(other, (list, numpy.ndarray, tuple)):
1733             if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
1734                 return self.__C * numpy.asmatrix(numpy.ravel(other)).T
1735             else:
1736                 return self.__C * numpy.asmatrix(other)
1737         elif self.isobject():
1738             return self.__C.__mul__(other)
1739         else:
1740             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
1741
1742     def __rmul__(self, other):
1743         "x.__rmul__(y) <==> y*x"
1744         if self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
1745             return other * self.__C
1746         elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
1747             if numpy.ravel(other).size == self.shape[1]: # Vecteur
1748                 return numpy.asmatrix(numpy.ravel(other)) * self.__C
1749             elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
1750                 return numpy.asmatrix(other) * self.__C
1751             else:
1752                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.asmatrix(other).shape,self.shape,self.__name))
1753         elif self.isvector() and isinstance(other,numpy.matrix):
1754             if numpy.ravel(other).size == self.shape[0]: # Vecteur
1755                 return numpy.asmatrix(numpy.ravel(other) * self.__C)
1756             elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
1757                 return numpy.asmatrix(numpy.array(other) * self.__C)
1758             else:
1759                 raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.ravel(other).shape,self.shape,self.__name))
1760         elif self.isscalar() and isinstance(other,numpy.matrix):
1761             return other * self.__C
1762         elif self.isobject():
1763             return self.__C.__rmul__(other)
1764         else:
1765             raise NotImplementedError("%s covariance matrix __rmul__ method not available for %s type!"%(self.__name,type(other)))
1766
1767     def __len__(self):
1768         "x.__len__() <==> len(x)"
1769         return self.shape[0]
1770
1771 # ==============================================================================
1772 class ObserverF(object):
1773     """
1774     Creation d'une fonction d'observateur a partir de son texte
1775     """
1776     def __init__(self, corps=""):
1777         self.__corps = corps
1778     def func(self,var,info):
1779         "Fonction d'observation"
1780         exec(self.__corps)
1781     def getfunc(self):
1782         "Restitution du pointeur de fonction dans l'objet"
1783         return self.func
1784
1785 # ==============================================================================
1786 class CaseLogger(object):
1787     """
1788     Conservation des commandes de creation d'un cas
1789     """
1790     def __init__(self, __name="", __objname="case", __addViewers=None, __addLoaders=None):
1791         self.__name     = str(__name)
1792         self.__objname  = str(__objname)
1793         self.__logSerie = []
1794         self.__switchoff = False
1795         self.__viewers = {
1796             "TUI" :Interfaces._TUIViewer,
1797             "SCD" :Interfaces._SCDViewer,
1798             "YACS":Interfaces._YACSViewer,
1799             }
1800         self.__loaders = {
1801             "TUI" :Interfaces._TUIViewer,
1802             "COM" :Interfaces._COMViewer,
1803             }
1804         if __addViewers is not None:
1805             self.__viewers.update(dict(__addViewers))
1806         if __addLoaders is not None:
1807             self.__loaders.update(dict(__addLoaders))
1808
1809     def register(self, __command=None, __keys=None, __local=None, __pre=None, __switchoff=False):
1810         "Enregistrement d'une commande individuelle"
1811         if __command is not None and __keys is not None and __local is not None and not self.__switchoff:
1812             if "self" in __keys: __keys.remove("self")
1813             self.__logSerie.append( (str(__command), __keys, __local, __pre, __switchoff) )
1814             if __switchoff:
1815                 self.__switchoff = True
1816         if not __switchoff:
1817             self.__switchoff = False
1818
1819     def dump(self, __filename=None, __format="TUI", __upa=""):
1820         "Restitution normalisée des commandes"
1821         if __format in self.__viewers:
1822             __formater = self.__viewers[__format](self.__name, self.__objname, self.__logSerie)
1823         else:
1824             raise ValueError("Dumping as \"%s\" is not available"%__format)
1825         return __formater.dump(__filename, __upa)
1826
1827     def load(self, __filename=None, __content=None, __object=None, __format="TUI"):
1828         "Chargement normalisé des commandes"
1829         if __format in self.__loaders:
1830             __formater = self.__loaders[__format]()
1831         else:
1832             raise ValueError("Loading as \"%s\" is not available"%__format)
1833         return __formater.load(__filename, __content, __object)
1834
1835 # ==============================================================================
1836 def MultiFonction( __xserie, _extraArguments = None, _sFunction = lambda x: x ):
1837     """
1838     Pour une liste ordonnée de vecteurs en entrée, renvoie en sortie la liste
1839     correspondante de valeurs de la fonction en argument
1840     """
1841     if not PlatformInfo.isIterable( __xserie ):
1842         raise TypeError("MultiFonction not iterable unkown input type: %s"%(type(__xserie),))
1843     #
1844     __multiHX = []
1845     if _extraArguments is None:
1846         for __xvalue in __xserie:
1847             __multiHX.append( _sFunction( __xvalue ) )
1848     elif _extraArguments is not None and isinstance(_extraArguments, (list, tuple, map)):
1849         for __xvalue in __xserie:
1850             __multiHX.append( _sFunction( __xvalue, *_extraArguments ) )
1851     elif _extraArguments is not None and isinstance(_extraArguments, dict):
1852         for __xvalue in __xserie:
1853             __multiHX.append( _sFunction( __xvalue, **_extraArguments ) )
1854     else:
1855         raise TypeError("MultiFonction extra arguments unkown input type: %s"%(type(_extraArguments),))
1856     #
1857     return __multiHX
1858
1859 # ==============================================================================
1860 def CostFunction3D(_x,
1861                    _Hm  = None,  # Pour simuler Hm(x) : HO["Direct"].appliedTo
1862                    _HmX = None,  # Simulation déjà faite de Hm(x)
1863                    _arg = None,  # Arguments supplementaires pour Hm, sous la forme d'un tuple
1864                    _BI  = None,
1865                    _RI  = None,
1866                    _Xb  = None,
1867                    _Y   = None,
1868                    _SIV = False, # A résorber pour la 8.0
1869                    _SSC = [],    # self._parameters["StoreSupplementaryCalculations"]
1870                    _nPS = 0,     # nbPreviousSteps
1871                    _QM  = "DA",  # QualityMeasure
1872                    _SSV = {},    # Entrée et/ou sortie : self.StoredVariables
1873                    _fRt = False, # Restitue ou pas la sortie étendue
1874                    _sSc = True,  # Stocke ou pas les SSC
1875                   ):
1876     """
1877     Fonction-coût générale utile pour les algorithmes statiques/3D : 3DVAR, BLUE
1878     et dérivés, Kalman et dérivés, LeastSquares, SamplingTest, PSO, SA, Tabu,
1879     DFO, QuantileRegression
1880     """
1881     if not _sSc:
1882         _SIV = False
1883         _SSC = {}
1884     else:
1885         for k in ["CostFunctionJ",
1886                   "CostFunctionJb",
1887                   "CostFunctionJo",
1888                   "CurrentOptimum",
1889                   "CurrentState",
1890                   "IndexOfOptimum",
1891                   "SimulatedObservationAtCurrentOptimum",
1892                   "SimulatedObservationAtCurrentState",
1893                  ]:
1894             if k not in _SSV:
1895                 _SSV[k] = []
1896             if hasattr(_SSV[k],"store"):
1897                 _SSV[k].append = _SSV[k].store # Pour utiliser "append" au lieu de "store"
1898     #
1899     _X  = numpy.asmatrix(numpy.ravel( _x )).T
1900     if _SIV or "CurrentState" in _SSC or "CurrentOptimum" in _SSC:
1901         _SSV["CurrentState"].append( _X )
1902     #
1903     if _HmX is not None:
1904         _HX = _HmX
1905     else:
1906         if _Hm is None:
1907             raise ValueError("COSTFUNCTION3D Operator has to be defined.")
1908         if _arg is None:
1909             _HX = _Hm( _X )
1910         else:
1911             _HX = _Hm( _X, *_arg )
1912     _HX = numpy.asmatrix(numpy.ravel( _HX )).T
1913     #
1914     if "SimulatedObservationAtCurrentState" in _SSC or \
1915        "SimulatedObservationAtCurrentOptimum" in _SSC:
1916         _SSV["SimulatedObservationAtCurrentState"].append( _HX )
1917     #
1918     if numpy.any(numpy.isnan(_HX)):
1919         Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
1920     else:
1921         _Y   = numpy.asmatrix(numpy.ravel( _Y )).T
1922         if _QM in ["AugmentedWeightedLeastSquares", "AWLS", "AugmentedPonderatedLeastSquares", "APLS", "DA"]:
1923             if _BI is None or _RI is None:
1924                 raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
1925             _Xb  = numpy.asmatrix(numpy.ravel( _Xb )).T
1926             Jb  = 0.5 * (_X - _Xb).T * _BI * (_X - _Xb)
1927             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1928         elif _QM in ["WeightedLeastSquares", "WLS", "PonderatedLeastSquares", "PLS"]:
1929             if _RI is None:
1930                 raise ValueError("Observation error covariance matrix has to be properly defined!")
1931             Jb  = 0.
1932             Jo  = 0.5 * (_Y - _HX).T * _RI * (_Y - _HX)
1933         elif _QM in ["LeastSquares", "LS", "L2"]:
1934             Jb  = 0.
1935             Jo  = 0.5 * (_Y - _HX).T * (_Y - _HX)
1936         elif _QM in ["AbsoluteValue", "L1"]:
1937             Jb  = 0.
1938             Jo  = numpy.sum( numpy.abs(_Y - _HX) )
1939         elif _QM in ["MaximumError", "ME"]:
1940             Jb  = 0.
1941             Jo  = numpy.max( numpy.abs(_Y - _HX) )
1942         elif _QM in ["QR", "Null"]:
1943             Jb  = 0.
1944             Jo  = 0.
1945         else:
1946             raise ValueError("Unknown asked quality measure!")
1947         #
1948         J   = float( Jb ) + float( Jo )
1949     #
1950     if _sSc:
1951         _SSV["CostFunctionJb"].append( Jb )
1952         _SSV["CostFunctionJo"].append( Jo )
1953         _SSV["CostFunctionJ" ].append( J )
1954     #
1955     if "IndexOfOptimum" in _SSC or \
1956        "CurrentOptimum" in _SSC or \
1957        "SimulatedObservationAtCurrentOptimum" in _SSC:
1958         IndexMin = numpy.argmin( _SSV["CostFunctionJ"][_nPS:] ) + _nPS
1959     if "IndexOfOptimum" in _SSC:
1960         _SSV["IndexOfOptimum"].append( IndexMin )
1961     if "CurrentOptimum" in _SSC:
1962         _SSV["CurrentOptimum"].append( _SSV["CurrentState"][IndexMin] )
1963     if "SimulatedObservationAtCurrentOptimum" in _SSC:
1964         _SSV["SimulatedObservationAtCurrentOptimum"].append( _SSV["SimulatedObservationAtCurrentState"][IndexMin] )
1965     #
1966     if _fRt:
1967         return _SSV
1968     else:
1969         if _QM in ["QR"]: # Pour le QuantileRegression
1970             return _HX
1971         else:
1972             return J
1973
1974 # ==============================================================================
1975 if __name__ == "__main__":
1976     print('\n AUTODIAGNOSTIC \n')