Compatibility correction for multiple numpy versions (REX [#25041])

dX correction for error with numpy 1.16.4 on DB09 (re)build

diff --git a/src/daComposant/daAlgorithms/AdjointTest.py b/src/daComposant/daAlgorithms/AdjointTest.py
index c9bbcecf710ae405447e1f7f998c6973081d970d..db8c6b34316759d26c95a5057200f1ce7bcdd606 100644 (file)
--- a/src/daComposant/daAlgorithms/AdjointTest.py
+++ b/src/daComposant/daAlgorithms/AdjointTest.py
@@ -21,7 +21,7 @@
 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
 
 import numpy
-from daCore import BasicObjects, PlatformInfo
+from daCore import BasicObjects, NumericObjects, PlatformInfo
 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
 
 # ==============================================================================
@@ -92,32 +92,26 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         Ht = HO["Tangent"].appliedInXTo
         Ha = HO["Adjoint"].appliedInXTo
         #
-        # ----------
         Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
         Perturbations.reverse()
         #
-        X       = numpy.ravel( Xb ).reshape((-1,1))
-        NormeX  = numpy.linalg.norm( X )
+        Xn       = numpy.ravel( Xb ).reshape((-1,1))
+        NormeX  = numpy.linalg.norm( Xn )
         if Y is None:
-            Y = numpy.ravel( Hm( X ) ).reshape((-1,1))
-        Y = numpy.ravel( Y ).reshape((-1,1))
-        NormeY = numpy.linalg.norm( Y )
+            Yn = numpy.ravel( Hm( Xn ) ).reshape((-1,1))
+        else:
+            Yn = numpy.ravel( Y ).reshape((-1,1))
+        NormeY = numpy.linalg.norm( Yn )
         if self._toStore("CurrentState"):
-            self.StoredVariables["CurrentState"].store( X )
+            self.StoredVariables["CurrentState"].store( Xn )
         if self._toStore("SimulatedObservationAtCurrentState"):
-            self.StoredVariables["SimulatedObservationAtCurrentState"].store( Y )
+            self.StoredVariables["SimulatedObservationAtCurrentState"].store( Yn )
         #
-        if len(self._parameters["InitialDirection"]) == 0:
-            dX0 = []
-            for v in X:
-                if abs(v) > 1.e-8:
-                    dX0.append( numpy.random.normal(0.,abs(v)) )
-                else:
-                    dX0.append( numpy.random.normal(0.,X.mean()) )
-        else:
-            dX0 = self._parameters["InitialDirection"]
-        #
-        dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.ravel( dX0 )
+        dX0 = NumericObjects.SetInitialDirection(
+            self._parameters["InitialDirection"],
+            self._parameters["AmplitudeOfInitialDirection"],
+            Xn,
+            )
         #
         # Entete des resultats
         # --------------------
@@ -156,10 +150,10 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             dX          = amplitude * dX0
             NormedX     = numpy.linalg.norm( dX )
             #
-            TangentFXdX = numpy.ravel( Ht( (X,dX) ) )
-            AdjointFXY  = numpy.ravel( Ha( (X,Y)  ) )
+            TangentFXdX = numpy.ravel( Ht( (Xn,dX) ) )
+            AdjointFXY  = numpy.ravel( Ha( (Xn,Yn)  ) )
             #
-            Residu = abs(float(numpy.dot( TangentFXdX, Y ) - numpy.dot( dX, AdjointFXY )))
+            Residu = abs(float(numpy.dot( TangentFXdX, Yn ) - numpy.dot( dX, AdjointFXY )))
             #
             msg = "  %2i  %5.0e   %9.3e   %9.3e   %9.3e   |  %9.3e"%(i,amplitude,NormeX,NormeY,NormedX,Residu)
             msgs += "\n" + __marge + msg

diff --git a/src/daComposant/daAlgorithms/GradientTest.py b/src/daComposant/daAlgorithms/GradientTest.py
index eb18c8fd5a99aea35e96be7a5540b88a56b9157d..0a410d04ed81bc45d178137bbbedf9dd3b628504 100644 (file)
--- a/src/daComposant/daAlgorithms/GradientTest.py
+++ b/src/daComposant/daAlgorithms/GradientTest.py
@@ -21,7 +21,7 @@
 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
 
 import math, numpy
-from daCore import BasicObjects, PlatformInfo
+from daCore import BasicObjects, NumericObjects, PlatformInfo
 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
 
 # ==============================================================================
@@ -125,25 +125,20 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         FX      = numpy.ravel( Hm( X ) ).reshape((-1,1))
         NormeX  = numpy.linalg.norm( X )
         NormeFX = numpy.linalg.norm( FX )
+        if NormeFX < mpr: NormeFX = mpr
         if self._toStore("CurrentState"):
             self.StoredVariables["CurrentState"].store( X )
         if self._toStore("SimulatedObservationAtCurrentState"):
             self.StoredVariables["SimulatedObservationAtCurrentState"].store( FX )
         #
-        if len(self._parameters["InitialDirection"]) == 0:
-            dX0 = []
-            for v in X:
-                if abs(v) > 1.e-8:
-                    dX0.append( numpy.random.normal(0.,abs(v)) )
-                else:
-                    dX0.append( numpy.random.normal(0.,X.mean()) )
-        else:
-            dX0 = numpy.ravel( self._parameters["InitialDirection"] )
-        #
-        dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.ravel( dX0 ).reshape((-1,1))
+        dX0 = NumericObjects.SetInitialDirection(
+            self._parameters["InitialDirection"],
+            self._parameters["AmplitudeOfInitialDirection"],
+            X,
+            )
         #
         if self._parameters["ResiduFormula"] in ["Taylor", "TaylorOnNorm"]:
-            dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
+            dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0.reshape((-1,1))
             GradFxdX = Ht( (X, dX1) )
             GradFxdX = numpy.ravel( GradFxdX ).reshape((-1,1))
             GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
@@ -236,7 +231,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         NormesdFXGdX = []
         #
         for i,amplitude in enumerate(Perturbations):
-            dX      = amplitude * dX0
+            dX      = amplitude * dX0.reshape((-1,1))
             #
             FX_plus_dX = Hm( X + dX )
             FX_plus_dX = numpy.ravel( FX_plus_dX ).reshape((-1,1))

diff --git a/src/daComposant/daAlgorithms/LinearityTest.py b/src/daComposant/daAlgorithms/LinearityTest.py
index fda2bc0ed98527162da1983691d09b15b6b794c5..6ae2c7fc652180058473bf0486561b86a83febe8 100644 (file)
--- a/src/daComposant/daAlgorithms/LinearityTest.py
+++ b/src/daComposant/daAlgorithms/LinearityTest.py
@@ -21,7 +21,7 @@
 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
 
 import math, numpy
-from daCore import BasicObjects, PlatformInfo
+from daCore import BasicObjects, NumericObjects, PlatformInfo
 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
 
 # ==============================================================================
@@ -99,44 +99,29 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
             import math
             return math.sqrt( ((numpy.ravel(V2) - numpy.ravel(V1))**2).sum() / float(numpy.ravel(V1).size) )
         #
-        # Operateurs
-        # ----------
         Hm = HO["Direct"].appliedTo
         if self._parameters["ResiduFormula"] in ["Taylor", "NominalTaylor", "NominalTaylorRMS"]:
             Ht = HO["Tangent"].appliedInXTo
         #
-        # Construction des perturbations
-        # ------------------------------
         Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
         Perturbations.reverse()
         #
-        # Calcul du point courant
-        # -----------------------
         Xn      = numpy.ravel(     Xb   ).reshape((-1,1))
         FX      = numpy.ravel( Hm( Xn ) ).reshape((-1,1))
         NormeX  = numpy.linalg.norm( Xn )
         NormeFX = numpy.linalg.norm( FX )
+        if NormeFX < mpr: NormeFX = mpr
         if self._toStore("CurrentState"):
             self.StoredVariables["CurrentState"].store( numpy.ravel(Xn) )
         if self._toStore("SimulatedObservationAtCurrentState"):
             self.StoredVariables["SimulatedObservationAtCurrentState"].store( numpy.ravel(FX) )
         #
-        # Fabrication de la direction de l'increment dX
-        # ---------------------------------------------
-        if len(self._parameters["InitialDirection"]) == 0:
-            dX0 = []
-            for v in Xn:
-                if abs(v) > 1.e-8:
-                    dX0.append( numpy.random.normal(0.,abs(v)) )
-                else:
-                    dX0.append( numpy.random.normal(0.,Xn.mean()) )
-        else:
-            dX0 = numpy.ravel( self._parameters["InitialDirection"] )
-        #
-        dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.ravel( dX0 ).reshape((-1,1))
+        dX0 = NumericObjects.SetInitialDirection(
+            self._parameters["InitialDirection"],
+            self._parameters["AmplitudeOfInitialDirection"],
+            Xn,
+            )
         #
-        # Calcul du gradient au point courant X pour l'increment dX
-        # ---------------------------------------------------------
         if self._parameters["ResiduFormula"] in ["Taylor", "NominalTaylor", "NominalTaylorRMS"]:
             dX1      = float(self._parameters["AmplitudeOfTangentPerturbation"]) * dX0
             GradFxdX = Ht( (Xn, dX1) )
@@ -252,7 +237,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         # Boucle sur les perturbations
         # ----------------------------
         for i,amplitude in enumerate(Perturbations):
-            dX      = amplitude * dX0
+            dX      = amplitude * dX0.reshape((-1,1))
             #
             if self._parameters["ResiduFormula"] == "CenteredDL":
                 if self._toStore("CurrentState"):

diff --git a/src/daComposant/daAlgorithms/TangentTest.py b/src/daComposant/daAlgorithms/TangentTest.py
index cd0565f0774127280387a3d09cc32db778e2098e..e9ecad0854cf36160e92487cb6ea3a17be827d0f 100644 (file)
--- a/src/daComposant/daAlgorithms/TangentTest.py
+++ b/src/daComposant/daAlgorithms/TangentTest.py
@@ -21,7 +21,7 @@
 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
 
 import numpy
-from daCore import BasicObjects, PlatformInfo
+from daCore import BasicObjects, NumericObjects, PlatformInfo
 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
 
 # ==============================================================================
@@ -114,19 +114,11 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         if self._toStore("SimulatedObservationAtCurrentState"):
             self.StoredVariables["SimulatedObservationAtCurrentState"].store( FX )
         #
-        # Fabrication de la direction de l'increment dX
-        # ---------------------------------------------
-        if len(self._parameters["InitialDirection"]) == 0:
-            dX0 = []
-            for v in Xn:
-                if abs(v) > 1.e-8:
-                    dX0.append( numpy.random.normal(0.,abs(v)) )
-                else:
-                    dX0.append( numpy.random.normal(0.,Xn.mean()) )
-        else:
-            dX0 = numpy.ravel( self._parameters["InitialDirection"] )
-        #
-        dX0 = float(self._parameters["AmplitudeOfInitialDirection"]) * numpy.ravel( dX0 ).reshape((-1,1))
+        dX0 = NumericObjects.SetInitialDirection(
+            self._parameters["InitialDirection"],
+            self._parameters["AmplitudeOfInitialDirection"],
+            Xn,
+            )
         #
         # Calcul du gradient au point courant X pour l'increment dX
         # qui est le tangent en X multiplie par dX
@@ -136,6 +128,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         GradFxdX = numpy.ravel( GradFxdX ).reshape((-1,1))
         GradFxdX = float(1./self._parameters["AmplitudeOfTangentPerturbation"]) * GradFxdX
         NormeGX  = numpy.linalg.norm( GradFxdX )
+        if NormeGX < mpr: NormeGX = mpr
         #
         # Entete des resultats
         # --------------------
@@ -184,7 +177,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         # Boucle sur les perturbations
         # ----------------------------
         for i,amplitude in enumerate(Perturbations):
-            dX      = amplitude * dX0
+            dX      = amplitude * dX0.reshape((-1,1))
             #
             if self._parameters["ResiduFormula"] == "Taylor":
                 FX_plus_dX  = numpy.ravel( Hm( Xn + dX ) ).reshape((-1,1))

diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py
index b6be0670dc6adf7c3843e68a6f7569973417e824..8f1a22649bf08bf9047b5e3648a2f018b11f0774 100644 (file)
--- a/src/daComposant/daCore/NumericObjects.py
+++ b/src/daComposant/daCore/NumericObjects.py
@@ -476,6 +476,34 @@ class FDApproximation(object):
             if self.__mfEnabled: return [_HaY,]
             else:                return _HaY
 
+# ==============================================================================
+def SetInitialDirection( __Direction = [], __Amplitude = 1., __Position = None ):
+    "Établit ou élabore une direction avec une amplitude"
+    #
+    if len(__Direction) == 0 and __Position is None:
+        raise ValueError("If initial direction is void, current position has to be given")
+    if abs(float(__Amplitude)) < mpr:
+        raise ValueError("Amplitude of perturbation can not be zero")
+    #
+    if len(__Direction) > 0:
+        __dX0 = numpy.asarray(__Direction)
+    else:
+        __dX0 = []
+        __X0 = numpy.ravel(numpy.asarray(__Position))
+        __mX0 = numpy.mean( __X0, dtype=mfp )
+        if abs(__mX0) < 2*mpr: __mX0 = 1. # Évite le problème de position nulle
+        for v in __X0:
+            if abs(v) > 1.e-8:
+                __dX0.append( numpy.random.normal(0.,abs(v)) )
+            else:
+                __dX0.append( numpy.random.normal(0.,__mX0) )
+    #
+    __dX0 = numpy.asarray(__dX0,float) # Évite le problème d'array de taille 1
+    __dX0 = numpy.ravel( __dX0 )       # Redresse les vecteurs
+    __dX0 = float(__Amplitude) * __dX0
+    #
+    return __dX0
+
 # ==============================================================================
 def EnsembleOfCenteredPerturbations( __bgCenter, __bgCovariance, __nbMembers ):
     "Génération d'un ensemble de taille __nbMembers-1 d'états aléatoires centrés"