Salome HOME
Minor improvements and fixes for internal variables
[modules/adao.git] / src / daComposant / daCore / BasicObjects.py
index 340612437730232b165d9a618594406ef50a0981..7b8953ac6cb8a0c3d3d811d2e512f379403717f0 100644 (file)
@@ -524,6 +524,7 @@ class FullOperator(object):
                 centeredDF            = __Function["CenteredFiniteDifference"],
                 increment             = __Function["DifferentialIncrement"],
                 dX                    = __Function["withdX"],
+                extraArguments        = self.__extraArgs,
                 avoidingRedundancy    = __Function["withAvoidingRedundancy"],
                 toleranceInRedundancy = __Function["withToleranceInRedundancy"],
                 lenghtOfRedundancy    = __Function["withLenghtOfRedundancy"],
@@ -623,6 +624,7 @@ class Algorithm(object):
             - OMB : Observation moins Background : Y - Xb
             - ForecastState : état prédit courant lors d'itérations
             - Residu : dans le cas des algorithmes de vérification
+            - SampledStateForQuantiles : échantillons d'états pour l'estimation des quantiles
             - SigmaBck2 : indicateur de correction optimale des erreurs d'ébauche
             - SigmaObs2 : indicateur de correction optimale des erreurs d'observation
             - SimulatedObservationAtBackground : l'état observé H(Xb) à l'ébauche
@@ -660,7 +662,8 @@ class Algorithm(object):
         self.StoredVariables["CostFunctionJbAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJbAtCurrentOptimum")
         self.StoredVariables["CostFunctionJo"]                       = Persistence.OneScalar(name = "CostFunctionJo")
         self.StoredVariables["CostFunctionJoAtCurrentOptimum"]       = Persistence.OneScalar(name = "CostFunctionJoAtCurrentOptimum")
-        self.StoredVariables["CurrentIterationNumber"]               = Persistence.OneIndex(name = "CurrentIterationNumber")
+        self.StoredVariables["CurrentEnsembleState"]                 = Persistence.OneMatrix(name = "CurrentEnsembleState")
+        self.StoredVariables["CurrentIterationNumber"]               = Persistence.OneIndex(name  = "CurrentIterationNumber")
         self.StoredVariables["CurrentOptimum"]                       = Persistence.OneVector(name = "CurrentOptimum")
         self.StoredVariables["CurrentState"]                         = Persistence.OneVector(name = "CurrentState")
         self.StoredVariables["ForecastState"]                        = Persistence.OneVector(name = "ForecastState")
@@ -679,6 +682,7 @@ class Algorithm(object):
         self.StoredVariables["OMA"]                                  = Persistence.OneVector(name = "OMA")
         self.StoredVariables["OMB"]                                  = Persistence.OneVector(name = "OMB")
         self.StoredVariables["Residu"]                               = Persistence.OneScalar(name = "Residu")
+        self.StoredVariables["SampledStateForQuantiles"]             = Persistence.OneMatrix(name = "SampledStateForQuantiles")
         self.StoredVariables["SigmaBck2"]                            = Persistence.OneScalar(name = "SigmaBck2")
         self.StoredVariables["SigmaObs2"]                            = Persistence.OneScalar(name = "SigmaObs2")
         self.StoredVariables["SimulatedObservationAtBackground"]     = Persistence.OneVector(name = "SimulatedObservationAtBackground")
@@ -1756,15 +1760,20 @@ class Covariance(object):
             __Matrix, __Scalar, __Vector, __Object = asCovariance, asEyeByScalar, asEyeByVector, asCovObject
         #
         if __Scalar is not None:
-            if numpy.matrix(__Scalar).size != 1:
-                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)
+            if isinstance(__Scalar, str):
+                __Scalar = __Scalar.replace(";"," ").replace(","," ").split()
+                if len(__Scalar) > 0: __Scalar = __Scalar[0]
+            if numpy.array(__Scalar).size != 1:
+                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.array(__Scalar).size)
             self.__is_scalar = True
             self.__C         = numpy.abs( float(__Scalar) )
             self.shape       = (0,0)
             self.size        = 0
         elif __Vector is not None:
+            if isinstance(__Vector, str):
+                __Vector = __Vector.replace(";"," ").replace(","," ").split()
             self.__is_vector = True
-            self.__C         = numpy.abs( numpy.array( numpy.ravel( numpy.matrix(__Vector, float ) ) ) )
+            self.__C         = numpy.abs( numpy.array( numpy.ravel( __Vector ), dtype=float ) )
             self.shape       = (self.__C.size,self.__C.size)
             self.size        = self.__C.size**2
         elif __Matrix is not None:
@@ -1775,7 +1784,7 @@ class Covariance(object):
         elif __Object is not None:
             self.__is_object = True
             self.__C         = __Object
-            for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__mul__","__rmul__"):
+            for at in ("getT","getI","diag","trace","__add__","__sub__","__neg__","__matmul__","__mul__","__rmatmul__","__rmul__"):
                 if not hasattr(self.__C,at):
                     raise ValueError("The matrix given for %s as an object has no attribute \"%s\". Please check your object input."%(self.__name,at))
             if hasattr(self.__C,"shape"):
@@ -1834,12 +1843,12 @@ class Covariance(object):
     def getI(self):
         "Inversion"
         if   self.ismatrix():
-            return Covariance(self.__name+"I", asCovariance  = self.__C.I )
+            return Covariance(self.__name+"I", asCovariance  = numpy.linalg.inv(self.__C) )
         elif self.isvector():
             return Covariance(self.__name+"I", asEyeByVector = 1. / self.__C )
         elif self.isscalar():
             return Covariance(self.__name+"I", asEyeByScalar = 1. / self.__C )
-        elif self.isobject():
+        elif self.isobject() and hasattr(self.__C,"getI"):
             return Covariance(self.__name+"I", asCovObject   = self.__C.getI() )
         else:
             return None # Indispensable
@@ -1852,8 +1861,10 @@ class Covariance(object):
             return Covariance(self.__name+"T", asEyeByVector = self.__C )
         elif self.isscalar():
             return Covariance(self.__name+"T", asEyeByScalar = self.__C )
-        elif self.isobject():
+        elif self.isobject() and hasattr(self.__C,"getT"):
             return Covariance(self.__name+"T", asCovObject   = self.__C.getT() )
+        else:
+            raise AttributeError("the %s covariance matrix has no getT attribute."%(self.__name,))
 
     def cholesky(self):
         "Décomposition de Cholesky"
@@ -1865,41 +1876,49 @@ class Covariance(object):
             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
         elif self.isobject() and hasattr(self.__C,"cholesky"):
             return Covariance(self.__name+"C", asCovObject   = self.__C.cholesky() )
+        else:
+            raise AttributeError("the %s covariance matrix has no cholesky attribute."%(self.__name,))
 
     def choleskyI(self):
         "Inversion de la décomposition de Cholesky"
         if   self.ismatrix():
-            return Covariance(self.__name+"H", asCovariance  = numpy.linalg.cholesky(self.__C).I )
+            return Covariance(self.__name+"H", asCovariance  = numpy.linalg.inv(numpy.linalg.cholesky(self.__C)) )
         elif self.isvector():
             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
         elif self.isscalar():
             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
         elif self.isobject() and hasattr(self.__C,"choleskyI"):
             return Covariance(self.__name+"H", asCovObject   = self.__C.choleskyI() )
+        else:
+            raise AttributeError("the %s covariance matrix has no choleskyI attribute."%(self.__name,))
 
     def sqrtm(self):
         "Racine carrée matricielle"
         if   self.ismatrix():
             import scipy
-            return Covariance(self.__name+"C", asCovariance  = scipy.linalg.sqrtm(self.__C) )
+            return Covariance(self.__name+"C", asCovariance  = numpy.real(scipy.linalg.sqrtm(self.__C)) )
         elif self.isvector():
             return Covariance(self.__name+"C", asEyeByVector = numpy.sqrt( self.__C ) )
         elif self.isscalar():
             return Covariance(self.__name+"C", asEyeByScalar = numpy.sqrt( self.__C ) )
-        elif self.isobject() and hasattr(self.__C,"sqrt"):
-            return Covariance(self.__name+"C", asCovObject   = self.__C.sqrt() )
+        elif self.isobject() and hasattr(self.__C,"sqrtm"):
+            return Covariance(self.__name+"C", asCovObject   = self.__C.sqrtm() )
+        else:
+            raise AttributeError("the %s covariance matrix has no sqrtm attribute."%(self.__name,))
 
     def sqrtmI(self):
         "Inversion de la racine carrée matricielle"
         if   self.ismatrix():
             import scipy
-            return Covariance(self.__name+"H", asCovariance  = scipy.linalg.sqrtm(self.__C).I )
+            return Covariance(self.__name+"H", asCovariance  = numpy.linalg.inv(numpy.real(scipy.linalg.sqrtm(self.__C))) )
         elif self.isvector():
             return Covariance(self.__name+"H", asEyeByVector = 1.0 / numpy.sqrt( self.__C ) )
         elif self.isscalar():
             return Covariance(self.__name+"H", asEyeByScalar = 1.0 / numpy.sqrt( self.__C ) )
-        elif self.isobject() and hasattr(self.__C,"sqrtI"):
-            return Covariance(self.__name+"H", asCovObject   = self.__C.sqrtI() )
+        elif self.isobject() and hasattr(self.__C,"sqrtmI"):
+            return Covariance(self.__name+"H", asCovObject   = self.__C.sqrtmI() )
+        else:
+            raise AttributeError("the %s covariance matrix has no sqrtmI attribute."%(self.__name,))
 
     def diag(self, msize=None):
         "Diagonale de la matrice"
@@ -1912,22 +1931,10 @@ class Covariance(object):
                 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,))
             else:
                 return self.__C * numpy.ones(int(msize))
-        elif self.isobject():
+        elif self.isobject() and hasattr(self.__C,"diag"):
             return self.__C.diag()
-
-    def asfullmatrix(self, msize=None):
-        "Matrice pleine"
-        if   self.ismatrix():
-            return self.__C
-        elif self.isvector():
-            return numpy.matrix( numpy.diag(self.__C), float )
-        elif self.isscalar():
-            if msize is None:
-                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,))
-            else:
-                return numpy.matrix( self.__C * numpy.eye(int(msize)), float )
-        elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
-            return self.__C.asfullmatrix()
+        else:
+            raise AttributeError("the %s covariance matrix has no diag attribute."%(self.__name,))
 
     def trace(self, msize=None):
         "Trace de la matrice"
@@ -1942,6 +1949,28 @@ class Covariance(object):
                 return self.__C * int(msize)
         elif self.isobject():
             return self.__C.trace()
+        else:
+            raise AttributeError("the %s covariance matrix has no trace attribute."%(self.__name,))
+
+    def asfullmatrix(self, msize=None):
+        "Matrice pleine"
+        if   self.ismatrix():
+            return numpy.asarray(self.__C)
+        elif self.isvector():
+            return numpy.asarray( numpy.diag(self.__C), float )
+        elif self.isscalar():
+            if msize is None:
+                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,))
+            else:
+                return numpy.asarray( self.__C * numpy.eye(int(msize)), float )
+        elif self.isobject() and hasattr(self.__C,"asfullmatrix"):
+            return self.__C.asfullmatrix()
+        else:
+            raise AttributeError("the %s covariance matrix has no asfullmatrix attribute."%(self.__name,))
+
+    def assparsematrix(self):
+        "Valeur sparse"
+        return self.__C
 
     def getO(self):
         return self
@@ -1984,6 +2013,36 @@ class Covariance(object):
         "x.__neg__() <==> -x"
         return - self.__C
 
+    def __matmul__(self, other):
+        "x.__mul__(y) <==> x@y"
+        if   self.ismatrix() and isinstance(other, (int, float)):
+            return numpy.asarray(self.__C) * other
+        elif self.ismatrix() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
+            if numpy.ravel(other).size == self.shape[1]: # Vecteur
+                return numpy.ravel(self.__C @ numpy.ravel(other))
+            elif numpy.asarray(other).shape[0] == self.shape[1]: # Matrice
+                return numpy.asarray(self.__C) @ numpy.asarray(other)
+            else:
+                raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.asarray(other).shape,self.__name))
+        elif self.isvector() and isinstance(other, (list, numpy.matrix, numpy.ndarray, tuple)):
+            if numpy.ravel(other).size == self.shape[1]: # Vecteur
+                return numpy.ravel(self.__C) * numpy.ravel(other)
+            elif numpy.asarray(other).shape[0] == self.shape[1]: # Matrice
+                return numpy.ravel(self.__C).reshape((-1,1)) * numpy.asarray(other)
+            else:
+                raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(self.shape,numpy.ravel(other).shape,self.__name))
+        elif self.isscalar() and isinstance(other,numpy.matrix):
+            return numpy.asarray(self.__C * other)
+        elif self.isscalar() and isinstance(other, (list, numpy.ndarray, tuple)):
+            if len(numpy.asarray(other).shape) == 1 or numpy.asarray(other).shape[1] == 1 or numpy.asarray(other).shape[0] == 1:
+                return self.__C * numpy.ravel(other)
+            else:
+                return self.__C * numpy.asarray(other)
+        elif self.isobject():
+            return self.__C.__matmul__(other)
+        else:
+            raise NotImplementedError("%s covariance matrix __matmul__ method not available for %s type!"%(self.__name,type(other)))
+
     def __mul__(self, other):
         "x.__mul__(y) <==> x*y"
         if   self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
@@ -2014,6 +2073,31 @@ class Covariance(object):
         else:
             raise NotImplementedError("%s covariance matrix __mul__ method not available for %s type!"%(self.__name,type(other)))
 
+    def __rmatmul__(self, other):
+        "x.__rmul__(y) <==> y@x"
+        if self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
+            return other * self.__C
+        elif self.ismatrix() and isinstance(other, (list, numpy.ndarray, tuple)):
+            if numpy.ravel(other).size == self.shape[1]: # Vecteur
+                return numpy.asmatrix(numpy.ravel(other)) * self.__C
+            elif numpy.asmatrix(other).shape[0] == self.shape[1]: # Matrice
+                return numpy.asmatrix(other) * self.__C
+            else:
+                raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.asmatrix(other).shape,self.shape,self.__name))
+        elif self.isvector() and isinstance(other,numpy.matrix):
+            if numpy.ravel(other).size == self.shape[0]: # Vecteur
+                return numpy.asmatrix(numpy.ravel(other) * self.__C)
+            elif numpy.asmatrix(other).shape[1] == self.shape[0]: # Matrice
+                return numpy.asmatrix(numpy.array(other) * self.__C)
+            else:
+                raise ValueError("operands could not be broadcast together with shapes %s %s in %s matrix"%(numpy.ravel(other).shape,self.shape,self.__name))
+        elif self.isscalar() and isinstance(other,numpy.matrix):
+            return other * self.__C
+        elif self.isobject():
+            return self.__C.__rmatmul__(other)
+        else:
+            raise NotImplementedError("%s covariance matrix __rmatmul__ method not available for %s type!"%(self.__name,type(other)))
+
     def __rmul__(self, other):
         "x.__rmul__(y) <==> y*x"
         if self.ismatrix() and isinstance(other, (int, numpy.matrix, float)):
@@ -2139,14 +2223,7 @@ def MultiFonction(
     #
     # Calculs effectifs
     if __mpEnabled:
-        _jobs = []
-        if _extraArguments is None:
-            _jobs = __xserie
-        elif _extraArguments is not None and isinstance(_extraArguments, (list, tuple, map)):
-            for __xvalue in __xserie:
-                _jobs.append( [__xvalue, ] + list(_extraArguments) )
-        else:
-            raise TypeError("MultiFonction extra arguments unkown input type: %s"%(type(_extraArguments),))
+        _jobs = __xserie
         # logging.debug("MULTF Internal multiprocessing calculations begin : evaluation of %i point(s)"%(len(_jobs),))
         import multiprocessing
         with multiprocessing.Pool(__mpWorkers) as pool: