src/daComposant/daAlgorithms/AdjointTest.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2023 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 import numpy
  24 from daCore import BasicObjects, NumericObjects, PlatformInfo
  25 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
  26 mfp = PlatformInfo.PlatformInfo().MaximumPrecision()
  27
  28 # ==============================================================================
  29 class ElementaryAlgorithm(BasicObjects.Algorithm):
  30     def __init__(self):
  31         BasicObjects.Algorithm.__init__(self, "ADJOINTTEST")
  32         self.defineRequiredParameter(
  33             name     = "ResiduFormula",
  34             default  = "ScalarProduct",
  35             typecast = str,
  36             message  = "Formule de résidu utilisée",
  37             listval  = ["ScalarProduct"],
  38             )
  39         self.defineRequiredParameter(
  40             name     = "AmplitudeOfInitialDirection",
  41             default  = 1.,
  42             typecast = float,
  43             message  = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal",
  44             )
  45         self.defineRequiredParameter(
  46             name     = "EpsilonMinimumExponent",
  47             default  = -8,
  48             typecast = int,
  49             message  = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément",
  50             minval   = -20,
  51             maxval   = 0,
  52             )
  53         self.defineRequiredParameter(
  54             name     = "InitialDirection",
  55             default  = [],
  56             typecast = list,
  57             message  = "Direction initiale de la dérivée directionnelle autour du point nominal",
  58             )
  59         self.defineRequiredParameter(
  60             name     = "NumberOfPrintedDigits",
  61             default  = 5,
  62             typecast = int,
  63             message  = "Nombre de chiffres affichés pour les impressions de réels",
  64             minval   = 0,
  65             )
  66         self.defineRequiredParameter(
  67             name     = "ResultTitle",
  68             default  = "",
  69             typecast = str,
  70             message  = "Titre du tableau et de la figure",
  71             )
  72         self.defineRequiredParameter(
  73             name     = "SetSeed",
  74             typecast = numpy.random.seed,
  75             message  = "Graine fixée pour le générateur aléatoire",
  76             )
  77         self.defineRequiredParameter(
  78             name     = "StoreSupplementaryCalculations",
  79             default  = [],
  80             typecast = tuple,
  81             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
  82             listval  = [
  83                 "CurrentState",
  84                 "Residu",
  85                 "SimulatedObservationAtCurrentState",
  86                 ]
  87             )
  88         self.requireInputArguments(
  89             mandatory= ("Xb", "HO" ),
  90             optional = ("Y", ),
  91             )
  92         self.setAttributes(tags=(
  93             "Checking",
  94             ))
  95
  96     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
  97         self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
  98         #
  99         Hm = HO["Direct"].appliedTo
 100         Ht = HO["Tangent"].appliedInXTo
 101         Ha = HO["Adjoint"].appliedInXTo
 102         #
 103         X0      = numpy.ravel( Xb ).reshape((-1,1))
 104         #
 105         # ----------
 106         __p = self._parameters["NumberOfPrintedDigits"]
 107         #
 108         __marge = 5*u" "
 109         __flech = 3*"="+"> "
 110         msgs  = ("\n") # 1
 111         if len(self._parameters["ResultTitle"]) > 0:
 112             __rt = str(self._parameters["ResultTitle"])
 113             msgs += (__marge + "====" + "="*len(__rt) + "====\n")
 114             msgs += (__marge + "    " + __rt + "\n")
 115             msgs += (__marge + "====" + "="*len(__rt) + "====\n")
 116         else:
 117             msgs += (__marge + "%s\n"%self._name)
 118             msgs += (__marge + "%s\n"%("="*len(self._name),))
 119         #
 120         msgs += ("\n")
 121         msgs += (__marge + "This test allows to analyze the quality of an adjoint operator associated\n")
 122         msgs += (__marge + "to some given direct operator F, applied to one single vector argument x.\n")
 123         msgs += (__marge + "If the adjoint operator is approximated and not given, the test measures\n")
 124         msgs += (__marge + "the quality of the automatic approximation, around an input checking point X.\n")
 125         msgs += ("\n")
 126         msgs += (__flech + "Information before launching:\n")
 127         msgs += (__marge + "-----------------------------\n")
 128         msgs += ("\n")
 129         msgs += (__marge + "Characteristics of input vector X, internally converted:\n")
 130         msgs += (__marge + "  Type...............: %s\n")%type( X0 )
 131         msgs += (__marge + "  Length of vector...: %i\n")%max(numpy.ravel( X0 ).shape)
 132         msgs += (__marge + "  Minimum value......: %."+str(__p)+"e\n")%numpy.min(  X0 )
 133         msgs += (__marge + "  Maximum value......: %."+str(__p)+"e\n")%numpy.max(  X0 )
 134         msgs += (__marge + "  Mean of vector.....: %."+str(__p)+"e\n")%numpy.mean( X0, dtype=mfp )
 135         msgs += (__marge + "  Standard error.....: %."+str(__p)+"e\n")%numpy.std(  X0, dtype=mfp )
 136         msgs += (__marge + "  L2 norm of vector..: %."+str(__p)+"e\n")%numpy.linalg.norm( X0 )
 137         msgs += ("\n")
 138         msgs += (__marge + "%s\n\n"%("-"*75,))
 139         msgs += (__flech + "Numerical quality indicators:\n")
 140         msgs += (__marge + "-----------------------------\n")
 141         msgs += ("\n")
 142         #
 143         if self._parameters["ResiduFormula"] == "ScalarProduct":
 144             msgs += (__marge + "Using the \"%s\" formula, one observes the residue R which is the\n"%self._parameters["ResiduFormula"])
 145             msgs += (__marge + "difference of two scalar products:\n")
 146             msgs += ("\n")
 147             msgs += (__marge + "    R(Alpha) = | < TangentF_X(dX) , Y > - < dX , AdjointF_X(Y) > |\n")
 148             msgs += ("\n")
 149             msgs += (__marge + "which must remain constantly equal to zero to the accuracy of the calculation.\n")
 150             msgs += (__marge + "One takes dX0 = Normal(0,X) and dX = Alpha*dX0, where F is the calculation\n")
 151             msgs += (__marge + "operator. If it is given, Y must be in the image of F. If it is not given,\n")
 152             msgs += (__marge + "one takes Y = F(X).\n")
 153             #
 154             __entete = str.rstrip("  i   Alpha  "     + \
 155                 str.center("||X||",2+__p+7)  + \
 156                 str.center("||Y||",2+__p+7)  + \
 157                 str.center("||dX||",2+__p+7) + \
 158                 str.center("R(Alpha)",2+__p+7))
 159             __nbtirets = len(__entete) + 2
 160             #
 161         msgs += ("\n")
 162         msgs += (__marge + "(Remark: numbers that are (about) under %.0e represent 0 to machine precision)\n"%mpr)
 163         print(msgs) # 1
 164         #
 165         Perturbations = [ 10**i for i in range(self._parameters["EpsilonMinimumExponent"],1) ]
 166         Perturbations.reverse()
 167         #
 168         NormeX  = numpy.linalg.norm( X0 )
 169         if Y is None:
 170             Yn = numpy.ravel( Hm( X0 ) ).reshape((-1,1))
 171         else:
 172             Yn = numpy.ravel( Y ).reshape((-1,1))
 173         NormeY = numpy.linalg.norm( Yn )
 174         if self._toStore("CurrentState"):
 175             self.StoredVariables["CurrentState"].store( X0 )
 176         if self._toStore("SimulatedObservationAtCurrentState"):
 177             self.StoredVariables["SimulatedObservationAtCurrentState"].store( Yn )
 178         #
 179         dX0 = NumericObjects.SetInitialDirection(
 180             self._parameters["InitialDirection"],
 181             self._parameters["AmplitudeOfInitialDirection"],
 182             X0,
 183             )
 184         #
 185         # Boucle sur les perturbations
 186         # ----------------------------
 187         msgs  = ("") # 2
 188         msgs += "\n" + __marge + "-"*__nbtirets
 189         msgs += "\n" + __marge + __entete
 190         msgs += "\n" + __marge + "-"*__nbtirets
 191         msgs += ("\n")
 192         __pf = "  %"+str(__p+7)+"."+str(__p)+"e"
 193         __ms = "  %2i  %5.0e"+(__pf*4)+"\n"
 194         for i,amplitude in enumerate(Perturbations):
 195             dX          = amplitude * dX0
 196             NormedX     = numpy.linalg.norm( dX )
 197             #
 198             if self._parameters["ResiduFormula"] == "ScalarProduct":
 199                 TangentFXdX = numpy.ravel( Ht( (X0,dX) ) )
 200                 AdjointFXY  = numpy.ravel( Ha( (X0,Yn)  ) )
 201                 #
 202                 Residu = abs(float(numpy.dot( TangentFXdX, Yn ) - numpy.dot( dX, AdjointFXY )))
 203                 #
 204                 self.StoredVariables["Residu"].store( Residu )
 205                 ttsep = __ms%(i,amplitude,NormeX,NormeY,NormedX,Residu)
 206                 msgs += __marge + ttsep
 207         #
 208         msgs += (__marge + "-"*__nbtirets + "\n\n")
 209         msgs += (__marge + "End of the \"%s\" verification by the \"%s\" formula.\n\n"%(self._name,self._parameters["ResiduFormula"]))
 210         msgs += (__marge + "%s\n"%("-"*75,))
 211         print(msgs) # 2
 212         #
 213         self._post_run(HO)
 214         return 0
 215
 216 # ==============================================================================
 217 if __name__ == "__main__":
 218     print('\n AUTODIAGNOSTIC\n')