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