From d2cc966dc54ab5bbee9d1ea3832042436c56519c Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Mon, 5 May 2014 08:42:30 +0200 Subject: [PATCH] Adding tangent test algorithm --- doc/en/reference.rst | 31 +++ doc/fr/reference.rst | 33 ++++ src/daComposant/daAlgorithms/TangentTest.py | 184 ++++++++++++++++++ .../daYacsSchemaCreator/infos_daComposant.py | 5 + 4 files changed, 253 insertions(+) create mode 100644 src/daComposant/daAlgorithms/TangentTest.py diff --git a/doc/en/reference.rst b/doc/en/reference.rst index 0861af8..0b89bc0 100644 --- a/doc/en/reference.rst +++ b/doc/en/reference.rst @@ -759,6 +759,7 @@ Optional and required commands for checking algorithms .. index:: single: GradientTest .. index:: single: LinearityTest .. index:: single: ObserverTest +.. index:: single: TangentTest .. index:: single: AlgorithmParameters .. index:: single: AmplitudeOfInitialDirection @@ -924,6 +925,36 @@ for an ADAO checking case`_. *"ObserverTest"*, and will not be used. The easiest way is to give "1" as a STRING for both, *"ObservationOperator"* having to be of type *Matrix*. +**"TangentTest"** + + *Required commands* + *"CheckingPoint", + "ObservationOperator"* + + AmplitudeOfInitialDirection + This key indicates the scaling of the initial perturbation build as a vector + used for the directional derivative around the nominal checking point. The + default is 1, that means no scaling. + + EpsilonMinimumExponent + This key indicates the minimal exponent value of the power of 10 coefficient + to be used to decrease the increment multiplier. The default is -8, and it + has to be between 0 and -20. For example, its default value leads to + calculate the residue of the scalar product formula with a fixed increment + multiplied from 1.e0 to 1.e-8. + + InitialDirection + This key indicates the vector direction used for the directional derivative + around the nominal checking point. It has to be a vector. If not specified, + this direction defaults to a random perturbation around zero of the same + vector size than the checking point. + + SetSeed + This key allow to give an integer in order to fix the seed of the random + generator used to generate the ensemble. A convenient value is for example + 1000. By default, the seed is left uninitialized, and so use the default + initialization from the computer. + Requirements for functions describing an operator ------------------------------------------------- diff --git a/doc/fr/reference.rst b/doc/fr/reference.rst index 0993e46..5fbaf73 100644 --- a/doc/fr/reference.rst +++ b/doc/fr/reference.rst @@ -796,6 +796,7 @@ Commandes optionnelles et requises pour les algorithmes de v .. index:: single: GradientTest .. index:: single: LinearityTest .. index:: single: ObserverTest +.. index:: single: TangentTest .. index:: single: AlgorithmParameters .. index:: single: AmplitudeOfInitialDirection @@ -969,6 +970,38 @@ commandes et mots-cl plus simple est de donner "1" comme un STRING pour les deux, l'*"ObservationOperator"* devant être de type *Matrix*. +**"TangentTest"** + + *Commandes obligatoires* + *"CheckingPoint", + "ObservationOperator"* + + AmplitudeOfInitialDirection + Cette clé indique la mise à l'échelle de la perturbation initiale construite + comme un vecteur utilisé pour la dérivée directionnelle autour du point + nominal de vérification. La valeur par défaut est de 1, ce qui signifie pas + de mise à l'échelle. + + EpsilonMinimumExponent + Cette clé indique la valeur de l'exposant minimal du coefficient en + puissance de 10 qui doit être utilisé pour faire décroître le multiplicateur + de l'incrément. La valeur par défaut est de -8, et elle doit être entre 0 et + -20. Par exemple, la valeur par défaut conduit à calculer le résidu de la + formule avec un incrément fixe multiplié par 1.e0 jusqu'à 1.e-8. + + InitialDirection + Cette clé indique la direction vectorielle utilisée pour la dérivée + directionnelle autour du point nominal de vérification. Cela doit être un + vecteur. Si elle n'est pas spécifiée, la direction par défaut est une + perturbation par défaut autour de zero de la même taille vectorielle que le + point de vérification. + + SetSeed + Cette clé permet de donner un nombre entier pour fixer la graine du + générateur aléatoire utilisé pour générer l'ensemble. Un valeur pratique est + par exemple 1000. Par défaut, la graine est laissée non initialisée, et elle + utilise ainsi l'initialisation par défaut de l'ordinateur. + Exigences pour les fonctions décrivant un opérateur --------------------------------------------------- diff --git a/src/daComposant/daAlgorithms/TangentTest.py b/src/daComposant/daAlgorithms/TangentTest.py new file mode 100644 index 0000000..1106e10 --- /dev/null +++ b/src/daComposant/daAlgorithms/TangentTest.py @@ -0,0 +1,184 @@ +#-*-coding:iso-8859-1-*- +# +# Copyright (C) 2008-2014 EDF R&D +# +# This library is free software; you can redistribute it and/or +# modify it under the terms of the GNU Lesser General Public +# License as published by the Free Software Foundation; either +# version 2.1 of the License. +# +# This library is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with this library; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +# +# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +# +# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D + +import logging +from daCore import BasicObjects +import numpy, math + +# ============================================================================== +class ElementaryAlgorithm(BasicObjects.Algorithm): + def __init__(self): + BasicObjects.Algorithm.__init__(self, "TANGENTTEST") + self.defineRequiredParameter( + name = "ResiduFormula", + default = "Taylor", + typecast = str, + message = "Formule de résidu utilisée", + listval = ["Taylor"], + ) + self.defineRequiredParameter( + name = "EpsilonMinimumExponent", + default = -8, + typecast = int, + message = "Exposant minimal en puissance de 10 pour le multiplicateur d'incrément", + minval = -20, + maxval = 0, + ) + self.defineRequiredParameter( + name = "InitialDirection", + default = [], + typecast = list, + message = "Direction initiale de la dérivée directionnelle autour du point nominal", + ) + self.defineRequiredParameter( + name = "AmplitudeOfInitialDirection", + default = 1., + typecast = float, + message = "Amplitude de la direction initiale de la dérivée directionnelle autour du point nominal", + ) + self.defineRequiredParameter( + name = "SetSeed", + typecast = numpy.random.seed, + message = "Graine fixée pour le générateur aléatoire", + ) + self.defineRequiredParameter( + name = "ResultTitle", + default = "", + typecast = str, + message = "Titre du tableau et de la figure", + ) + + def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None): + self._pre_run() + # + # Paramètres de pilotage + # ---------------------- + self.setParameters(Parameters) + # + # Opérateurs + # ---------- + Hm = HO["Direct"].appliedTo + Ht = HO["Tangent"].appliedInXTo + # + # Construction des perturbations + # ------------------------------ + Perturbations = [ 10**i for i in xrange(self._parameters["EpsilonMinimumExponent"],1) ] + Perturbations.reverse() + # + # Calcul du point courant + # ----------------------- + Xn = numpy.asmatrix(numpy.ravel( Xb )).T + FX = numpy.asmatrix(numpy.ravel( Hm( Xn ) )).T + NormeX = numpy.linalg.norm( Xn ) + NormeFX = numpy.linalg.norm( FX ) + # + # Fabrication de la direction de l'incrément dX + # ---------------------------------------------- + if len(self._parameters["InitialDirection"]) == 0: + dX0 = [] + for v in Xn.A1: + 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.matrix( dX0 ).T + # + # Calcul du gradient au point courant X pour l'incrément dX + # qui est le tangent en X multiplié par dX + # --------------------------------------------------------- + GradFxdX = Ht( (Xn, dX0) ) + GradFxdX = numpy.asmatrix(numpy.ravel( GradFxdX )).T + NormeGX = numpy.linalg.norm( GradFxdX ) + # + # Entete des resultats + # -------------------- + __marge = 12*" " + if self._parameters["ResiduFormula"] == "Taylor": + __entete = " i Alpha ||X|| ||F(X)|| | R(Alpha) |R-1|/Alpha " + __msgdoc = """ + On observe le résidu provenant du rapport d'incréments utilisant le + linéaire tangent : + + || F(X+Alpha*dX) - F(X) || + R(Alpha) = ----------------------------- + || Alpha * TangentF_X * dX || + + qui doit rester stable en 1+O(Alpha) jusqu'à ce que l'on atteigne la + précision du calcul. + + Lorsque |R-1|/Alpha est inférieur ou égal à une valeur stable + lorsque Alpha varie, le tangent est valide, jusqu'à ce que l'on + atteigne la précision du calcul. + + Si |R-1|/Alpha est très faible, le code F est vraisemblablement + linéaire ou quasi-linéaire, et le tangent est valide jusqu'à ce que + l'on atteigne la précision du calcul. + + On prend dX0 = Normal(0,X) et dX = Alpha*dX0. F est le code de calcul. + """ + # + if len(self._parameters["ResultTitle"]) > 0: + msgs = "\n" + msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n" + msgs += __marge + " " + self._parameters["ResultTitle"] + "\n" + msgs += __marge + "====" + "="*len(self._parameters["ResultTitle"]) + "====\n" + else: + msgs = "" + msgs += __msgdoc + # + __nbtirets = len(__entete) + msgs += "\n" + __marge + "-"*__nbtirets + msgs += "\n" + __marge + __entete + msgs += "\n" + __marge + "-"*__nbtirets + # + # Boucle sur les perturbations + # ---------------------------- + for i,amplitude in enumerate(Perturbations): + dX = amplitude * dX0 + # + if self._parameters["ResiduFormula"] == "Taylor": + FX_plus_dX = numpy.asmatrix(numpy.ravel( Hm( Xn + dX ) )).T + # + Residu = numpy.linalg.norm( FX_plus_dX - FX ) / (amplitude * NormeGX) + # + self.StoredVariables["CostFunctionJ"].store( Residu ) + msg = " %2i %5.0e %9.3e %9.3e | %11.5e %5.1e"%(i,amplitude,NormeX,NormeFX,Residu,abs(Residu-1.)/amplitude) + msgs += "\n" + __marge + msg + # + msgs += "\n" + __marge + "-"*__nbtirets + msgs += "\n" + # + # Sorties eventuelles + # ------------------- + print + print "Results of tangent check by \"%s\" formula:"%self._parameters["ResiduFormula"] + print msgs + # + self._post_run(HO) + return 0 + +# ============================================================================== +if __name__ == "__main__": + print '\n AUTODIAGNOSTIC \n' diff --git a/src/daSalome/daYacsSchemaCreator/infos_daComposant.py b/src/daSalome/daYacsSchemaCreator/infos_daComposant.py index 326da38..64431a0 100644 --- a/src/daSalome/daYacsSchemaCreator/infos_daComposant.py +++ b/src/daSalome/daYacsSchemaCreator/infos_daComposant.py @@ -76,6 +76,7 @@ CheckAlgos = [ "GradientTest", "AdjointTest", "ObserverTest", + "TangentTest", ] AlgoDataRequirements = {} @@ -152,6 +153,10 @@ AlgoDataRequirements["AdjointTest"] = [ AlgoDataRequirements["ObserverTest"] = [ "Observers", ] +AlgoDataRequirements["TangentTest"] = [ + "CheckingPoint", + "ObservationOperator", + ] AlgoType = {} AlgoType["3DVAR"] = "Optim" -- 2.39.2