From 98bd7f51f7339d6348492046668804040bfdf85a Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Wed, 23 Mar 2022 22:27:30 +0100 Subject: [PATCH] Compatibility correction for multiple numpy versions (REX [#25041]) dX correction for error with numpy 1.16.4 on DB09 (re)build --- src/daComposant/daAlgorithms/AdjointTest.py | 40 ++++++++----------- src/daComposant/daAlgorithms/GradientTest.py | 23 +++++------ src/daComposant/daAlgorithms/LinearityTest.py | 31 ++++---------- src/daComposant/daAlgorithms/TangentTest.py | 23 ++++------- src/daComposant/daCore/NumericObjects.py | 28 +++++++++++++ 5 files changed, 70 insertions(+), 75 deletions(-) diff --git a/src/daComposant/daAlgorithms/AdjointTest.py b/src/daComposant/daAlgorithms/AdjointTest.py index c9bbcec..db8c6b3 100644 --- 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 eb18c8f..0a410d0 100644 --- 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 fda2bc0..6ae2c7f 100644 --- 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 cd0565f..e9ecad0 100644 --- 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 b6be067..8f1a226 100644 --- 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" -- 2.30.2