src/daComposant/daCore/BasicObjects.py

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