From e5a70dbf1ff528a1428da3d7637ab02a26451597 Mon Sep 17 00:00:00 2001 From: Jean-Philippe ARGAUD Date: Thu, 8 Sep 2022 13:30:14 +0200 Subject: [PATCH] Code and documentation update --- doc/en/ref_algorithm_FunctionTest.rst | 13 +- doc/en/ref_algorithm_ParallelFunctionTest.rst | 4 +- doc/en/ref_observers_requirements.rst | 23 +- doc/en/scripts/simple_FunctionTest.rst | 28 +- .../scripts/simple_ParallelFunctionTest.rst | 39 +- doc/en/snippets/NoUnconditionalOutput.rst | 1 + doc/en/snippets/ShowElementarySummary.rst | 9 + doc/fr/ref_algorithm_AdjointTest.rst | 25 +- doc/fr/ref_algorithm_ExtendedBlue.rst | 5 +- doc/fr/ref_algorithm_FunctionTest.rst | 14 +- doc/fr/ref_algorithm_GradientTest.rst | 11 +- doc/fr/ref_algorithm_LinearityTest.rst | 16 +- doc/fr/ref_algorithm_ParallelFunctionTest.rst | 16 +- doc/fr/ref_algorithm_TangentTest.rst | 13 +- doc/fr/ref_observers_requirements.rst | 46 +- doc/fr/scripts/simple_FunctionTest.rst | 22 +- .../scripts/simple_ParallelFunctionTest.rst | 32 +- doc/fr/snippets/NoUnconditionalOutput.rst | 1 + doc/fr/snippets/ShowElementarySummary.rst | 9 + src/daComposant/daAlgorithms/Atoms/ecwnlls.py | 2 + .../daAlgorithms/Atoms/incr3dvar.py | 2 + .../daAlgorithms/Atoms/lbfgsb19hlt.py | 501 ++++++++++++++++++ src/daComposant/daAlgorithms/Atoms/mmqr.py | 6 +- .../daAlgorithms/Atoms/psas3dvar.py | 2 + .../daAlgorithms/Atoms/std3dvar.py | 2 + .../daAlgorithms/Atoms/std4dvar.py | 2 + .../daAlgorithms/Atoms/van3dvar.py | 2 + src/daComposant/daCore/Aidsm.py | 1 - src/daComposant/daCore/AssimilationStudy.py | 1 - src/daComposant/daCore/BasicObjects.py | 1 - src/daComposant/daCore/ExtendedLogging.py | 1 - src/daComposant/daCore/Interfaces.py | 1 - src/daComposant/daCore/NumericObjects.py | 1 - src/daComposant/daCore/Persistence.py | 1 - src/daComposant/daCore/PlatformInfo.py | 2 - test/test1002/Performances.py | 15 +- 36 files changed, 740 insertions(+), 130 deletions(-) create mode 100644 doc/en/snippets/NoUnconditionalOutput.rst create mode 100644 doc/en/snippets/ShowElementarySummary.rst create mode 100644 doc/fr/snippets/NoUnconditionalOutput.rst create mode 100644 doc/fr/snippets/ShowElementarySummary.rst create mode 100644 src/daComposant/daAlgorithms/Atoms/lbfgsb19hlt.py diff --git a/doc/en/ref_algorithm_FunctionTest.rst b/doc/en/ref_algorithm_FunctionTest.rst index ca4866a..71c6e15 100644 --- a/doc/en/ref_algorithm_FunctionTest.rst +++ b/doc/en/ref_algorithm_FunctionTest.rst @@ -30,11 +30,10 @@ Checking algorithm "*FunctionTest*" .. ------------------------------------ .. .. include:: snippets/Header2Algo01.rst -This algorithm allows to verify that an operator, in particular the -observation one, is working correctly and that its call is compatible -with its usage in ADAO algorithms. In practice, it allows to call one -or several times the operator, activating or not the "debug" mode -during execution. +This algorithm allows to verify that a given operator :math:`F`, in particular +the observation one, is working correctly and that its call is compatible with +its usage in ADAO algorithms. In practice, it allows to call one or several +times the operator, activating or not the "debug" mode during execution. Statistics on input and output vectors for each execution of operator are given, and an another global statistic is given at the end of the checking @@ -59,6 +58,8 @@ themselves beforehand with the intended test .. include:: snippets/SetDebug.rst +.. include:: snippets/ShowElementarySummary.rst + StoreSupplementaryCalculations .. index:: single: StoreSupplementaryCalculations @@ -82,7 +83,7 @@ StoreSupplementaryCalculations .. ------------------------------------ .. .. include:: snippets/Header2Algo04.rst -*None* +.. include:: snippets/NoUnconditionalOutput.rst .. ------------------------------------ .. .. include:: snippets/Header2Algo05.rst diff --git a/doc/en/ref_algorithm_ParallelFunctionTest.rst b/doc/en/ref_algorithm_ParallelFunctionTest.rst index df13223..dece3a8 100644 --- a/doc/en/ref_algorithm_ParallelFunctionTest.rst +++ b/doc/en/ref_algorithm_ParallelFunctionTest.rst @@ -59,6 +59,8 @@ themselves beforehand with the intended test .. include:: snippets/SetDebug.rst +.. include:: snippets/ShowElementarySummary.rst + StoreSupplementaryCalculations .. index:: single: StoreSupplementaryCalculations @@ -82,7 +84,7 @@ StoreSupplementaryCalculations .. ------------------------------------ .. .. include:: snippets/Header2Algo04.rst -*None* +.. include:: snippets/NoUnconditionalOutput.rst .. ------------------------------------ .. .. include:: snippets/Header2Algo05.rst diff --git a/doc/en/ref_observers_requirements.rst b/doc/en/ref_observers_requirements.rst index 59f5849..6804c05 100644 --- a/doc/en/ref_observers_requirements.rst +++ b/doc/en/ref_observers_requirements.rst @@ -33,7 +33,7 @@ Requirements for functions describing an "*observer*" Some special variables, internal to the optimization process and used inside calculation, can be monitored during an ADAO calculation. These variables can be printed, plotted, saved, etc. by the user. This can be done using some -"*observer*", sometimes also called "callback" on a variable. They are special +"*observer*", sometimes also called "callback", on a variable. They are special Python functions, each one associated with a given variable, as conceptually described in the following figure: @@ -61,10 +61,10 @@ figure: .. centered:: **Choosing its entry type for an "observer" function** -The "*observer*" function can be given as an explicit script (entry of type +An "*observer*" function can be given as an explicit script (entry of type "*String*"), as a script in an external file (entry of type "*Script*"), or by using a template or pattern (entry of type"*Template*"). The templates are -available by default in ADAO using the graphical interface EFICAS or the text +available by default in ADAO, using the graphical interface EFICAS or the text interface TUI, and are detailed in the following :ref:`section_ref_observers_templates`. These templates are simple scripts that can be tuned by the user, either in the integrated edition stage of the case @@ -92,16 +92,16 @@ To use directly this "*observer*" capability, the user must use or build a script that have on standard input (that is, in the naming space) the variables ``var`` and ``info``. The variable ``var`` is to be used as an object of list/tuple type, that contains the history of the variable of interest, indexed -by the iterating steps. Only the body of the "*observer*" function has to be -specified by the user, not the function call itself. +by the iterating and/or time steps. Only the body of the "*observer*" function +has to be specified by the user, not the Python ``def`` function call itself. As an example, here is a very simple script (similar to the "*ValuePrinter*" template), that can be used to print the value of the monitored variable:: print(" --->",info," Value =",var[-1]) -Stored as a Python file or as an explicit string, these script lines can be -associated to each variable found in the keyword "*SELECTION*" of the +Stored as a Python file or as an explicit string, this or these script lines +can be associated to each variable found in the keyword "*SELECTION*" of the "*Observers*" command of the ADAO case: "*Analysis*", "*CurrentState*", "*CostFunction*"... The current value of the variable will for example be printed at each step of the optimization or data assimilation algorithm. The @@ -117,6 +117,15 @@ execution of this "*observer*" is simply never activated. crash before being registered as an "*observer*" function. The debugging can otherwise be really difficult! +Some "*observer*" allow the creation of successive files or figures, which are +uniquely numbered and, if applicable, stored by default in the standard +``/tmp`` directory. In the case where this information needs to be modified (as +for example when the ``/tmp`` directory is a virtual or local non-permanent +folder, or when one wishes to have a numbering according to the iteration), the +user is encouraged to take inspiration from a model that is suitable for him +and to modify it by specifying differently this shared information. Then, the +modified function can be used in a "*String*" or "*Script*" input. + Hereinafter we give the identifier and the contents of all the available "*observer*" models. diff --git a/doc/en/scripts/simple_FunctionTest.rst b/doc/en/scripts/simple_FunctionTest.rst index 74a376b..e147d6f 100644 --- a/doc/en/scripts/simple_FunctionTest.rst +++ b/doc/en/scripts/simple_FunctionTest.rst @@ -1,16 +1,16 @@ -.. index:: single: FunctionTest (exemple) - -Cet exemple décrit le test du bon fonctionnement d'un opérateur et que son -appel se déroule de manière compatible avec son usage dans les algorithmes -d'ADAO. Les information nécessaires sont minimales, à savoir ici un opérateur -de type observation :math:`H` et un état :math:`\mathbf{x}^b` sur lequel le -tester (nommé "*CheckingPoint*" pour le test). - -Le test est répété un nombre paramétrable de fois, et une statistique finale -permet de vérifier rapidement le bon comportement de l'opérateur. Le diagnostic -le plus simple consiste à vérifier, à la fin, l'ordre de grandeur des valeurs -indiquées comme la moyenne des différences entre les sorties répétées et leur -moyenne ("*mean of the differences between the outputs Y and their mean Ym*"). -Pour un opérateur normal, ces valeurs doivent être proches du zéro numérique. +.. index:: single: FunctionTest (example) +This example describes the test of the correct operation of a given operator, +and that its call proceeds in a way compatible with its common use in the ADAO +algorithms. The required information are minimal, namely here an operator +:math:`F` (described for the test by the observation command +"*ObservationOperator*"), and a state :math:`\mathbf{x}^b` to test it on +(described for the test by the command "*CheckingPoint*"). +The test is repeated a configurable number of times, and a final statistic +makes it possible to quickly verify the operator's good behavior. The simplest +diagnostic consists in checking, at the very end of the display, the order of +magnitude of the values indicated as the mean of the differences between the +repeated outputs and their mean, under the part entitled "*Characteristics of +the mean of the differences between the outputs Y and their mean Ym*". For a +satisfactory operator, these values should be close to the numerical zero. diff --git a/doc/en/scripts/simple_ParallelFunctionTest.rst b/doc/en/scripts/simple_ParallelFunctionTest.rst index b068a99..eb902e9 100644 --- a/doc/en/scripts/simple_ParallelFunctionTest.rst +++ b/doc/en/scripts/simple_ParallelFunctionTest.rst @@ -1,14 +1,29 @@ .. index:: single: FunctionTest (example) -This example describes the test of the correct operation of an operator and -that its call is compatible with its use in the ADAO algorithms. The necessary -information are minimal, namely here an operator of type observation :math:`H` -and a state :math:`\mathbf{x}^b` on which to test it (named "*CheckingPoint*" -for the test). - -The test is repeated a customizable number of times, and a final statistic -allows to quickly check the good behavior of the operator. The simplest -diagnostic consists in checking, at the end, the order of magnitude of the -values indicated as the average of the differences between the repeated outputs -Y and their mean Ym*. For a typical operator, these values should be close to -the numerical zero. +This example describes the test of the correct operation of a given operator, +and that its call proceeds in a way compatible with its common use in parallel +in the ADAO algorithms. The required information are minimal, namely here an +operator :math:`F` (described for the test by the observation command +"*ObservationOperator*"), and a state :math:`\mathbf{x}^b` to test it on +(described for the test by the command "*CheckingPoint*"). + +The test is repeated a configurable number of times, and a final statistic +makes it possible to quickly verify the operator's good behavior. The simplest +diagnostic consists in checking, at the very end of the display, the order of +magnitude of the values indicated as the mean of the differences between the +repeated outputs and their mean, under the part entitled "*Characteristics of +the mean of the differences between the outputs Y and their mean Ym*". For a +satisfactory operator, these values should be close to the numerical zero. + +.. note:: + + .. index:: single: EnableMultiProcessingInEvaluation + + It can be useful to make sure that the evaluation of the operator is really + done in parallel, and for example that there is no forced use of a + parallelism acceleration, which would avoid a real parallel test. For this + purpose, it is recommended to systematically use the boolean special + parameter "*EnableMultiProcessingInEvaluation*", exclusively reserved for + this purpose, of the operator declaration command. The use of this + parameter is illustrated in this example. It should not be used in any + other case. diff --git a/doc/en/snippets/NoUnconditionalOutput.rst b/doc/en/snippets/NoUnconditionalOutput.rst new file mode 100644 index 0000000..6bd824c --- /dev/null +++ b/doc/en/snippets/NoUnconditionalOutput.rst @@ -0,0 +1 @@ + *None* diff --git a/doc/en/snippets/ShowElementarySummary.rst b/doc/en/snippets/ShowElementarySummary.rst new file mode 100644 index 0000000..8ffd5f5 --- /dev/null +++ b/doc/en/snippets/ShowElementarySummary.rst @@ -0,0 +1,9 @@ +.. index:: single: ShowElementarySummary + +ShowElementarySummary + *Boolean value*. This variable leads to the activation, or not, of the + calculation and display of a summary at each elementary evaluation of the + test. The default value is "True", the choices are "True" or "False". + + Example : + ``{"ShowElementarySummary":False}`` diff --git a/doc/fr/ref_algorithm_AdjointTest.rst b/doc/fr/ref_algorithm_AdjointTest.rst index 8402291..91f4865 100644 --- a/doc/fr/ref_algorithm_AdjointTest.rst +++ b/doc/fr/ref_algorithm_AdjointTest.rst @@ -30,24 +30,35 @@ Algorithme de vérification "*AdjointTest*" .. ------------------------------------ .. .. include:: snippets/Header2Algo01.rst -Cet algorithme permet de vérifier la qualité de l'opérateur adjoint, en -calculant un résidu dont les propriétés théoriques sont connues. +Cet algorithme permet de vérifier la qualité de l'adjoint d'un opérateur +:math:`F`, en calculant un résidu dont les propriétés théoriques sont connues. +Le test est applicable à un opérateur quelconque, d'évolution comme +d'observation. + +Pour toutes les formules, avec :math:`\mathbf{x}` le point courant de +vérification, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et +:math:`\mathbf{dx}=\alpha_0*\mathbf{dx}_0` avec :math:`\alpha_0` un paramètre +utilisateur de mise à l'échelle, par défaut à 1. :math:`F` est l'opérateur ou +le code de calcul (qui est ici acquis par la commande d'opérateur d'observation +"*ObservationOperator*"). On observe le résidu suivant, qui est la différence de deux produits scalaires : .. math:: R(\alpha) = | < TangentF_x(\mathbf{dx}) , \mathbf{y} > - < \mathbf{dx} , AdjointF_x(\mathbf{y}) > | -qui doit rester constamment égal à zéro à la précision du calcul. On prend -:math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et -:math:`\mathbf{dx}=\alpha*\mathbf{dx}_0`. :math:`F` est le code de calcul. -:math:`\mathbf{y}` doit être dans l'image de :math:`F`. S'il n'est pas donné, -on prend :math:`\mathbf{y} = F(\mathbf{x})`. +dans lequel la quantité optionnelle :math:`\mathbf{y}` doit être dans l'image +de :math:`F`. Si elle n'est pas donnée, on prend son évaluation par défaut +:math:`\mathbf{y} = F(\mathbf{x})`. + +Ce résidu doit rester constamment égal à zéro à la précision du calcul. .. ------------------------------------ .. .. include:: snippets/Header2Algo02.rst .. include:: snippets/CheckingPoint.rst +.. include:: snippets/Observation.rst + .. include:: snippets/ObservationOperator.rst .. ------------------------------------ .. diff --git a/doc/fr/ref_algorithm_ExtendedBlue.rst b/doc/fr/ref_algorithm_ExtendedBlue.rst index 6ee234f..1ee1821 100644 --- a/doc/fr/ref_algorithm_ExtendedBlue.rst +++ b/doc/fr/ref_algorithm_ExtendedBlue.rst @@ -34,8 +34,9 @@ Cet algorithme réalise une estimation de type BLUE étendu (Best Linear Unbiase Estimator, étendu) de l'état d'un système. Cet algorithme est une généralisation partiellement non-linéaire d'un -:ref:`section_ref_algorithm_Blue`. Il lui est équivalent pour un opérateur -d'observation linéaire. On peut vérifier la linéarité de l'opérateur +:ref:`section_ref_algorithm_Blue`. Si l'opérateur d'observation est +explicitement linéaire, l'algorithme est équivalent à celui du +:ref:`section_ref_algorithm_Blue`. On peut vérifier la linéarité de l'opérateur d'observation à l'aide d'un :ref:`section_ref_algorithm_LinearityTest`. En non-linéaire, ses résultats se rapprochent d'un diff --git a/doc/fr/ref_algorithm_FunctionTest.rst b/doc/fr/ref_algorithm_FunctionTest.rst index fe7436d..26aad55 100644 --- a/doc/fr/ref_algorithm_FunctionTest.rst +++ b/doc/fr/ref_algorithm_FunctionTest.rst @@ -30,11 +30,11 @@ Algorithme de vérification "*FunctionTest*" .. ------------------------------------ .. .. include:: snippets/Header2Algo01.rst -Cet algorithme permet de vérifier qu'un opérateur, dont en particulier celui -d'observation, fonctionne correctement et que son appel se déroule de manière -compatible avec son usage dans les algorithmes d'ADAO. De manière pratique, il -permet d'appeler une ou plusieurs fois l'opérateur, en activant ou non le mode -"debug" lors de l'exécution. +Cet algorithme permet de vérifier qu'un opérateur :math:`F` quelconque, dont en +particulier celui d'observation, fonctionne correctement et que son appel se +déroule de manière compatible avec son usage dans les algorithmes d'ADAO. De +manière pratique, il permet d'appeler une ou plusieurs fois l'opérateur, en +activant ou non le mode "debug" lors de l'exécution. Une statistique sur les vecteurs en entrée et en sortie de chaque exécution de l'opérateur est indiquée, et une autre globale est fournie de manière @@ -59,6 +59,8 @@ elles-mêmes avec le test prévu :ref:`section_ref_algorithm_InputValuesTest`. .. include:: snippets/SetDebug.rst +.. include:: snippets/ShowElementarySummary.rst + StoreSupplementaryCalculations .. index:: single: StoreSupplementaryCalculations @@ -82,7 +84,7 @@ StoreSupplementaryCalculations .. ------------------------------------ .. .. include:: snippets/Header2Algo04.rst -*Aucune* +.. include:: snippets/NoUnconditionalOutput.rst .. ------------------------------------ .. .. include:: snippets/Header2Algo05.rst diff --git a/doc/fr/ref_algorithm_GradientTest.rst b/doc/fr/ref_algorithm_GradientTest.rst index 41249fb..b01d387 100644 --- a/doc/fr/ref_algorithm_GradientTest.rst +++ b/doc/fr/ref_algorithm_GradientTest.rst @@ -32,12 +32,15 @@ Algorithme de vérification "*GradientTest*" Cet algorithme permet de vérifier la qualité du gradient de l'opérateur, en calculant un résidu dont les propriétés théoriques sont connues. Plusieurs -formules de résidu sont disponibles. +formules de résidu sont disponibles. Le test est applicable à un opérateur +quelconque, d'évolution comme d'observation. -Dans tous les cas, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et +Pour toutes les formules, avec :math:`\mathbf{x}` le point courant de +vérification, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et :math:`\mathbf{dx}=\alpha_0*\mathbf{dx}_0` avec :math:`\alpha_0` un paramètre -utilisateur de mise à l'échelle, par défaut à 1. :math:`F` est le code de -calcul. +utilisateur de mise à l'échelle, par défaut à 1. :math:`F` est l'opérateur ou +le code de calcul (qui est ici acquis par la commande d'opérateur d'observation +"*ObservationOperator*"). Résidu "Taylor" *************** diff --git a/doc/fr/ref_algorithm_LinearityTest.rst b/doc/fr/ref_algorithm_LinearityTest.rst index 715813b..3636081 100644 --- a/doc/fr/ref_algorithm_LinearityTest.rst +++ b/doc/fr/ref_algorithm_LinearityTest.rst @@ -30,12 +30,18 @@ Algorithme de vérification "*LinearityTest*" .. ------------------------------------ .. .. include:: snippets/Header2Algo01.rst -Cet algorithme permet de vérifier la qualité de linéarité de l'opérateur, en +Cet algorithme permet de vérifier la qualité de linéarité d'un opérateur, en calculant un résidu dont les propriétés théoriques sont connues. Plusieurs -formules de résidu sont utilisables. - -Dans tous les cas, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et -:math:`\mathbf{dx}=\alpha*\mathbf{dx}_0`. :math:`F` est le code de calcul. +formules de résidu sont utilisables et sont décrites ci-dessous avec leur +interprétation. Le test est applicable à un opérateur quelconque, d'évolution +comme d'observation. + +Pour toutes les formules, avec :math:`\mathbf{x}` le point courant de +vérification, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et +:math:`\mathbf{dx}=\alpha_0*\mathbf{dx}_0` avec :math:`\alpha_0` un paramètre +utilisateur de mise à l'échelle, par défaut à 1. :math:`F` est l'opérateur ou +le code de calcul (qui est ici acquis par la commande d'opérateur d'observation +"*ObservationOperator*"). Résidu "CenteredDL" ******************* diff --git a/doc/fr/ref_algorithm_ParallelFunctionTest.rst b/doc/fr/ref_algorithm_ParallelFunctionTest.rst index b2932df..c363a77 100644 --- a/doc/fr/ref_algorithm_ParallelFunctionTest.rst +++ b/doc/fr/ref_algorithm_ParallelFunctionTest.rst @@ -30,12 +30,12 @@ Algorithme de vérification "*ParallelFunctionTest*" .. ------------------------------------ .. .. include:: snippets/Header2Algo01.rst -Cet algorithme permet de vérifier qu'un opérateur, dont en particulier celui -d'observation, fonctionne -correctement en parallèle et que son appel se déroule de manière compatible -avec son usage dans les algorithmes d'ADAO. De manière pratique, il permet -d'appeler une ou plusieurs fois l'opérateur en parallèle, en activant ou non le -mode "debug" lors de l'exécution. +Cet algorithme permet de vérifier qu'un opérateur :math:`F` quelconque, dont en +particulier celui d'observation, fonctionne correctement en parallèle et que +son appel se déroule de manière compatible avec son usage dans les algorithmes +d'ADAO. De manière pratique, il permet d'appeler une ou plusieurs fois +l'opérateur en parallèle, en activant ou non le mode "debug" lors de +l'exécution. Une statistique sur les vecteurs en entrée et en sortie de chaque exécution de l'opérateur est indiquée, et une autre globale est fournie de manière @@ -60,6 +60,8 @@ elles-mêmes avec le test prévu :ref:`section_ref_algorithm_InputValuesTest`. .. include:: snippets/SetDebug.rst +.. include:: snippets/ShowElementarySummary.rst + StoreSupplementaryCalculations .. index:: single: StoreSupplementaryCalculations @@ -83,7 +85,7 @@ StoreSupplementaryCalculations .. ------------------------------------ .. .. include:: snippets/Header2Algo04.rst -*Aucune* +.. include:: snippets/NoUnconditionalOutput.rst .. ------------------------------------ .. .. include:: snippets/Header2Algo05.rst diff --git a/doc/fr/ref_algorithm_TangentTest.rst b/doc/fr/ref_algorithm_TangentTest.rst index 485ee55..31d2ac0 100644 --- a/doc/fr/ref_algorithm_TangentTest.rst +++ b/doc/fr/ref_algorithm_TangentTest.rst @@ -31,7 +31,15 @@ Algorithme de vérification "*TangentTest*" .. include:: snippets/Header2Algo01.rst Cet algorithme permet de vérifier la qualité de l'opérateur tangent, en -calculant un résidu dont les propriétés théoriques sont connues. +calculant un résidu dont les propriétés théoriques sont connues. Le test est +applicable à un opérateur quelconque, d'évolution comme d'observation. + +Pour toutes les formules, avec :math:`\mathbf{x}` le point courant de +vérification, on prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et +:math:`\mathbf{dx}=\alpha_0*\mathbf{dx}_0` avec :math:`\alpha_0` un paramètre +utilisateur de mise à l'échelle, par défaut à 1. :math:`F` est l'opérateur ou +le code de calcul (qui est ici acquis par la commande d'opérateur d'observation +"*ObservationOperator*"). On observe le résidu suivant, provenant du rapport d'incréments utilisant l'opérateur linéaire tangent : @@ -50,9 +58,6 @@ vraisemblablement linéaire ou quasi-linéaire (ce que l'on peut vérifier par l':ref:`section_ref_algorithm_LinearityTest`), et le tangent est valide jusqu'à ce que l'on atteigne la précision du calcul. -On prend :math:`\mathbf{dx}_0=Normal(0,\mathbf{x})` et -:math:`\mathbf{dx}=\alpha*\mathbf{dx}_0`. :math:`F` est le code de calcul. - .. ------------------------------------ .. .. include:: snippets/Header2Algo02.rst diff --git a/doc/fr/ref_observers_requirements.rst b/doc/fr/ref_observers_requirements.rst index e7aaed1..a9ff7d0 100644 --- a/doc/fr/ref_observers_requirements.rst +++ b/doc/fr/ref_observers_requirements.rst @@ -33,7 +33,7 @@ Conditions requises pour les fonctions décrivant un "*observer*" Certaines variables spéciales, internes à l'optimisation et utilisées au cours des calculs, peuvent être surveillées durant un calcul ADAO. Ces variables peuvent être affichées, tracées, enregistrées, etc. par l'utilisateur. C'est -réalisable en utilisant des "*observer*", parfois aussi appelés des "callback" +réalisable en utilisant des "*observer*", parfois aussi appelés des "callback", sur une variable. Ce sont des fonctions Python spéciales, qui sont chacune associées à une variable donnée, comme décrit conceptuellement dans la figure suivante : @@ -65,10 +65,10 @@ montré dans la figure qui suit : Une fonction "*observer*" peut être fourni sous la forme d'un script explicite (entrée de type "*String*"), d'un script contenu dans un fichier externe (entrée de type "*Script*"), ou en utilisant un modèle (entrée de type -"*Template*"). Les modèles sont fournis par défaut dans ADAO lors de l'usage de -l'éditeur graphique EFICAS d'ADAO ou de l'interface TUI, et sont détaillés dans -la partie :ref:`section_ref_observers_templates` qui suit. Ces derniers sont -des scripts simples qui peuvent être adaptés par l'utilisateur, soit dans +"*Template*"). Les modèles sont fournis par défaut dans ADAO, lors de l'usage +de l'éditeur graphique EFICAS d'ADAO ou de l'interface TUI, et sont détaillés +dans la partie :ref:`section_ref_observers_templates` qui suit. Ces derniers +sont des scripts simples qui peuvent être adaptés par l'utilisateur, soit dans l'étape d'édition intégrée du cas avec EFICAS d'ADAO, soit dans l'étape d'édition du schéma avant l'exécution, pour améliorer la performance du calcul ADAO dans le superviseur d'exécution de SALOME. @@ -96,9 +96,9 @@ Pour pouvoir utiliser directement cette capacité "*observer*", l'utilisateur doit utiliser ou construire un script utilisant en entrée standard (i.e. disponible dans l'espace de nommage) les variables ``var`` et ``info``. La variable ``var`` est à utiliser comme un objet de type liste/tuple, contenant -l'historique de la variable d'intérêt, indicé par les pas d'itérations. Seul le -corps de la fonction "*observer*" doit être spécifié par l'utilisateur, pas -l'appel de fonction lui-même. +l'historique de la variable d'intérêt, indicé par les pas d'itérations et/ou de +temps. Seul le corps de la fonction "*observer*" doit être spécifié par +l'utilisateur, pas l'appel Python ``def`` de fonction lui-même. A titre d'exemple, voici un script très simple (similaire au modèle "*ValuePrinter*"), utilisable pour afficher la valeur d'une variable @@ -107,16 +107,16 @@ surveillée : print(" --->",info," Value =",var[-1]) -Stockées comme un fichier Python ou une chaîne de caractères explicite, ces -lignes de script peuvent être associées à chaque variable présente dans le -mot-clé "*SELECTION*" de la commande "*Observers*" du cas ADAO : "*Analysis*", -"*CurrentState*", "*CostFunction*"... La valeur courante de la variable sera -par exemple affichée à chaque étape de l'algorithme d'optimisation ou -d'assimilation. Les "*observer*" peuvent inclure des capacités d'affichage -graphique, de stockage, de traitement complexe, d'analyse statistique, etc. Si -une variable, à laquelle est lié un "*observer*", n'est pas requise dans le -calcul et par l'utilisateur, l'exécution de cet "*observer*" n'est tout -simplement jamais activée. +Stockées comme un fichier Python ou une chaîne de caractères explicite, cette +ou ces lignes de script peuvent être associées à chaque variable présente dans +le mot-clé "*SELECTION*" de la commande "*Observers*" du cas ADAO : +"*Analysis*", "*CurrentState*", "*CostFunction*"... La valeur courante de la +variable sera par exemple affichée à chaque étape de l'algorithme +d'optimisation ou d'assimilation. Les "*observer*" peuvent inclure des +capacités d'affichage graphique, de stockage, de traitement complexe, d'analyse +statistique, etc. Si une variable, à laquelle est lié un "*observer*", n'est +pas requise dans le calcul et par l'utilisateur, l'exécution de cet +"*observer*" n'est tout simplement jamais activée. .. warning:: @@ -126,6 +126,16 @@ simplement jamais activée. comme une fonction "*observer*". Le débogage peut sinon être vraiment difficile ! +Certains "*observer*" permettent de créer des fichiers ou des figures +successives, qui sont numérotées de manière unique et, le cas échéant, +enregistrées par défaut dans le répertoire standard ``/tmp``. Dans le cas où +ces informations sont à modifier (comme par exemple lorsque le répertoire +``/tmp`` est un dossier virtuel ou local non pérenne, ou lorsque l'on désire +une numérotation en fonction de l'itération), l'utilisateur est invité à +s'inspirer d'un modèle lui convenant pour le modifier en spécifiant +différemment ces informations communes. Ensuite, la fonction modifiée peut être +utilisée dans une entrée de type "*String*" ou de type "*Script*". + On donne ci-après l'identifiant et le contenu de tous les modèles "*observer*" disponibles. diff --git a/doc/fr/scripts/simple_FunctionTest.rst b/doc/fr/scripts/simple_FunctionTest.rst index 74a376b..3152f78 100644 --- a/doc/fr/scripts/simple_FunctionTest.rst +++ b/doc/fr/scripts/simple_FunctionTest.rst @@ -1,16 +1,16 @@ .. index:: single: FunctionTest (exemple) -Cet exemple décrit le test du bon fonctionnement d'un opérateur et que son -appel se déroule de manière compatible avec son usage dans les algorithmes -d'ADAO. Les information nécessaires sont minimales, à savoir ici un opérateur -de type observation :math:`H` et un état :math:`\mathbf{x}^b` sur lequel le -tester (nommé "*CheckingPoint*" pour le test). +Cet exemple décrit le test du bon fonctionnement d'un opérateur quelconque, et +que son appel se déroule de manière compatible avec son usage courant dans les +algorithmes d'ADAO. Les information nécessaires sont minimales, à savoir ici un +opérateur :math:`F` (décrit pour le test par la commande d'observation +"*ObservationOperator*"), et un état :math:`\mathbf{x}^b` sur lequel le tester +(décrit pour le test par la commande "*CheckingPoint*"). Le test est répété un nombre paramétrable de fois, et une statistique finale permet de vérifier rapidement le bon comportement de l'opérateur. Le diagnostic -le plus simple consiste à vérifier, à la fin, l'ordre de grandeur des valeurs -indiquées comme la moyenne des différences entre les sorties répétées et leur -moyenne ("*mean of the differences between the outputs Y and their mean Ym*"). -Pour un opérateur normal, ces valeurs doivent être proches du zéro numérique. - - +le plus simple consiste à vérifier, à la toute fin de l'affichage, l'ordre de +grandeur des valeurs indiquées comme la moyenne des différences entre les +sorties répétées et leur moyenne, sous la partie titrée "*Characteristics of +the mean of the differences between the outputs Y and their mean Ym*". Pour un +opérateur satisfaisant, ces valeurs doivent être proches du zéro numérique. diff --git a/doc/fr/scripts/simple_ParallelFunctionTest.rst b/doc/fr/scripts/simple_ParallelFunctionTest.rst index 74a376b..8548be4 100644 --- a/doc/fr/scripts/simple_ParallelFunctionTest.rst +++ b/doc/fr/scripts/simple_ParallelFunctionTest.rst @@ -1,16 +1,30 @@ .. index:: single: FunctionTest (exemple) -Cet exemple décrit le test du bon fonctionnement d'un opérateur et que son -appel se déroule de manière compatible avec son usage dans les algorithmes -d'ADAO. Les information nécessaires sont minimales, à savoir ici un opérateur -de type observation :math:`H` et un état :math:`\mathbf{x}^b` sur lequel le -tester (nommé "*CheckingPoint*" pour le test). +Cet exemple décrit le test du bon fonctionnement d'un opérateur quelconque, et +que son appel se déroule de manière compatible avec son usage courant en +parallèle dans les algorithmes d'ADAO. Les information nécessaires sont +minimales, à savoir ici un opérateur :math:`F` (décrit pour le test par la +commande d'observation "*ObservationOperator*"), et un état +:math:`\mathbf{x}^b` sur lequel le tester (décrit pour le test par la commande +"*CheckingPoint*"). Le test est répété un nombre paramétrable de fois, et une statistique finale permet de vérifier rapidement le bon comportement de l'opérateur. Le diagnostic -le plus simple consiste à vérifier, à la fin, l'ordre de grandeur des valeurs -indiquées comme la moyenne des différences entre les sorties répétées et leur -moyenne ("*mean of the differences between the outputs Y and their mean Ym*"). -Pour un opérateur normal, ces valeurs doivent être proches du zéro numérique. +le plus simple consiste à vérifier, à la toute fin de l'affichage, l'ordre de +grandeur des valeurs indiquées comme la moyenne des différences entre les +sorties répétées et leur moyenne, sous la partie titrée "*Characteristics of +the mean of the differences between the outputs Y and their mean Ym*". Pour un +opérateur satisfaisant, ces valeurs doivent être proches du zéro numérique. +.. note:: + .. index:: single: EnableMultiProcessingInEvaluation + + Il peut être utile de s'assurer que l'évaluation de l'opérateur est + réalisée réellement en parallèle, et par exemple qu'il n'y a pas + d'utilisation forcée d'une accélération du parallélisme, qui éviterait + ainsi un véritable test parallèle. Pour cela, il est recommandé d'utiliser + systématiquement le paramètre booléen spécial + "*EnableMultiProcessingInEvaluation*", exclusivement réservé à cet usage, + de la commande de déclaration de l'opérateur. L'usage de ce paramètre est + illustré dans l'exemple présent. Il n'est à utiliser dans aucun autre cas. diff --git a/doc/fr/snippets/NoUnconditionalOutput.rst b/doc/fr/snippets/NoUnconditionalOutput.rst new file mode 100644 index 0000000..7edf811 --- /dev/null +++ b/doc/fr/snippets/NoUnconditionalOutput.rst @@ -0,0 +1 @@ + *Aucune* diff --git a/doc/fr/snippets/ShowElementarySummary.rst b/doc/fr/snippets/ShowElementarySummary.rst new file mode 100644 index 0000000..f26e4c0 --- /dev/null +++ b/doc/fr/snippets/ShowElementarySummary.rst @@ -0,0 +1,9 @@ +.. index:: single: ShowElementarySummary + +ShowElementarySummary + *Valeur booléenne*. La variable conduit à l'activation, ou pas, du calcul et + de l'affichage d'un résumé à chaque évaluation élémentaire du test. La valeur + par défaut est "True", les choix sont "True" ou "False". + + Exemple : + ``{"ShowElementarySummary":False}`` diff --git a/src/daComposant/daAlgorithms/Atoms/ecwnlls.py b/src/daComposant/daAlgorithms/Atoms/ecwnlls.py index c73f5fc..4a96ab7 100644 --- a/src/daComposant/daAlgorithms/Atoms/ecwnlls.py +++ b/src/daComposant/daAlgorithms/Atoms/ecwnlls.py @@ -138,6 +138,8 @@ def ecwnlls(selfA, Xb, Y, U, HO, CM, R, B, __storeState = False): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daAlgorithms/Atoms/incr3dvar.py b/src/daComposant/daAlgorithms/Atoms/incr3dvar.py index 42fa7e9..984886a 100644 --- a/src/daComposant/daAlgorithms/Atoms/incr3dvar.py +++ b/src/daComposant/daAlgorithms/Atoms/incr3dvar.py @@ -140,6 +140,8 @@ def incr3dvar(selfA, Xb, Y, U, HO, CM, R, B, __storeState = False): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daAlgorithms/Atoms/lbfgsb19hlt.py b/src/daComposant/daAlgorithms/Atoms/lbfgsb19hlt.py new file mode 100644 index 0000000..d8d7b6c --- /dev/null +++ b/src/daComposant/daAlgorithms/Atoms/lbfgsb19hlt.py @@ -0,0 +1,501 @@ +# Modification de la version 1.9.1 +""" +Functions +--------- +.. autosummary:: + :toctree: generated/ + + fmin_l_bfgs_b + +""" + +## License for the Python wrapper +## ============================== + +## Copyright (c) 2004 David M. Cooke + +## Permission is hereby granted, free of charge, to any person obtaining a +## copy of this software and associated documentation files (the "Software"), +## to deal in the Software without restriction, including without limitation +## the rights to use, copy, modify, merge, publish, distribute, sublicense, +## and/or sell copies of the Software, and to permit persons to whom the +## Software is furnished to do so, subject to the following conditions: + +## The above copyright notice and this permission notice shall be included in +## all copies or substantial portions of the Software. + +## THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +## IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +## FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +## AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +## LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +## FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +## DEALINGS IN THE SOFTWARE. + +## Modifications by Travis Oliphant and Enthought, Inc. for inclusion in SciPy + +import numpy as np +from numpy import array, asarray, float64, zeros +from scipy.optimize import _lbfgsb +from scipy.optimize._optimize import (MemoizeJac, OptimizeResult, + _check_unknown_options, _prepare_scalar_function) +from scipy.optimize._constraints import old_bound_to_new + +from scipy.sparse.linalg import LinearOperator + +__all__ = ['fmin_l_bfgs_b', 'LbfgsInvHessProduct'] + + +def fmin_l_bfgs_b(func, x0, fprime=None, args=(), + approx_grad=0, + bounds=None, m=10, factr=1e7, pgtol=1e-5, + epsilon=1e-8, + iprint=-1, maxfun=15000, maxiter=15000, disp=None, + callback=None, maxls=20): + """ + Minimize a function func using the L-BFGS-B algorithm. + + Parameters + ---------- + func : callable f(x,*args) + Function to minimize. + x0 : ndarray + Initial guess. + fprime : callable fprime(x,*args), optional + The gradient of `func`. If None, then `func` returns the function + value and the gradient (``f, g = func(x, *args)``), unless + `approx_grad` is True in which case `func` returns only ``f``. + args : sequence, optional + Arguments to pass to `func` and `fprime`. + approx_grad : bool, optional + Whether to approximate the gradient numerically (in which case + `func` returns only the function value). + bounds : list, optional + ``(min, max)`` pairs for each element in ``x``, defining + the bounds on that parameter. Use None or +-inf for one of ``min`` or + ``max`` when there is no bound in that direction. + m : int, optional + The maximum number of variable metric corrections + used to define the limited memory matrix. (The limited memory BFGS + method does not store the full hessian but uses this many terms in an + approximation to it.) + factr : float, optional + The iteration stops when + ``(f^k - f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= factr * eps``, + where ``eps`` is the machine precision, which is automatically + generated by the code. Typical values for `factr` are: 1e12 for + low accuracy; 1e7 for moderate accuracy; 10.0 for extremely + high accuracy. See Notes for relationship to `ftol`, which is exposed + (instead of `factr`) by the `scipy.optimize.minimize` interface to + L-BFGS-B. + pgtol : float, optional + The iteration will stop when + ``max{|proj g_i | i = 1, ..., n} <= pgtol`` + where ``pg_i`` is the i-th component of the projected gradient. + epsilon : float, optional + Step size used when `approx_grad` is True, for numerically + calculating the gradient + iprint : int, optional + Controls the frequency of output. ``iprint < 0`` means no output; + ``iprint = 0`` print only one line at the last iteration; + ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations; + ``iprint = 99`` print details of every iteration except n-vectors; + ``iprint = 100`` print also the changes of active set and final x; + ``iprint > 100`` print details of every iteration including x and g. + disp : int, optional + If zero, then no output. If a positive number, then this over-rides + `iprint` (i.e., `iprint` gets the value of `disp`). + maxfun : int, optional + Maximum number of function evaluations. Note that this function + may violate the limit because of evaluating gradients by numerical + differentiation. + maxiter : int, optional + Maximum number of iterations. + callback : callable, optional + Called after each iteration, as ``callback(xk)``, where ``xk`` is the + current parameter vector. + maxls : int, optional + Maximum number of line search steps (per iteration). Default is 20. + + Returns + ------- + x : array_like + Estimated position of the minimum. + f : float + Value of `func` at the minimum. + d : dict + Information dictionary. + + * d['warnflag'] is + + - 0 if converged, + - 1 if too many function evaluations or too many iterations, + - 2 if stopped for another reason, given in d['task'] + + * d['grad'] is the gradient at the minimum (should be 0 ish) + * d['funcalls'] is the number of function calls made. + * d['nit'] is the number of iterations. + + See also + -------- + minimize: Interface to minimization algorithms for multivariate + functions. See the 'L-BFGS-B' `method` in particular. Note that the + `ftol` option is made available via that interface, while `factr` is + provided via this interface, where `factr` is the factor multiplying + the default machine floating-point precision to arrive at `ftol`: + ``ftol = factr * numpy.finfo(float).eps``. + + Notes + ----- + License of L-BFGS-B (FORTRAN code): + + The version included here (in fortran code) is 3.0 + (released April 25, 2011). It was written by Ciyou Zhu, Richard Byrd, + and Jorge Nocedal . It carries the following + condition for use: + + This software is freely available, but we expect that all publications + describing work using this software, or all commercial products using it, + quote at least one of the references given below. This software is released + under the BSD License. + + References + ---------- + * R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound + Constrained Optimization, (1995), SIAM Journal on Scientific and + Statistical Computing, 16, 5, pp. 1190-1208. + * C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B, + FORTRAN routines for large scale bound constrained optimization (1997), + ACM Transactions on Mathematical Software, 23, 4, pp. 550 - 560. + * J.L. Morales and J. Nocedal. L-BFGS-B: Remark on Algorithm 778: L-BFGS-B, + FORTRAN routines for large scale bound constrained optimization (2011), + ACM Transactions on Mathematical Software, 38, 1. + + """ + # handle fprime/approx_grad + if approx_grad: + fun = func + jac = None + elif fprime is None: + fun = MemoizeJac(func) + jac = fun.derivative + else: + fun = func + jac = fprime + + # build options + if disp is None: + disp = iprint + opts = {'disp': disp, + 'iprint': iprint, + 'maxcor': m, + 'ftol': factr * np.finfo(float).eps, + 'gtol': pgtol, + 'eps': epsilon, + 'maxfun': maxfun, + 'maxiter': maxiter, + 'callback': callback, + 'maxls': maxls} + + res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds, + **opts) + d = {'grad': res['jac'], + 'task': res['message'], + 'funcalls': res['nfev'], + 'nit': res['nit'], + 'warnflag': res['status']} + f = res['fun'] + x = res['x'] + + return x, f, d + + +def _minimize_lbfgsb(fun, x0, args=(), jac=None, bounds=None, + disp=None, maxcor=10, ftol=2.2204460492503131e-09, + gtol=1e-5, eps=1e-8, maxfun=15000, maxiter=15000, + iprint=-1, callback=None, maxls=20, + finite_diff_rel_step=None, **unknown_options): + """ + Minimize a scalar function of one or more variables using the L-BFGS-B + algorithm. + + Options + ------- + disp : None or int + If `disp is None` (the default), then the supplied version of `iprint` + is used. If `disp is not None`, then it overrides the supplied version + of `iprint` with the behaviour you outlined. + maxcor : int + The maximum number of variable metric corrections used to + define the limited memory matrix. (The limited memory BFGS + method does not store the full hessian but uses this many terms + in an approximation to it.) + ftol : float + The iteration stops when ``(f^k - + f^{k+1})/max{|f^k|,|f^{k+1}|,1} <= ftol``. + gtol : float + The iteration will stop when ``max{|proj g_i | i = 1, ..., n} + <= gtol`` where ``pg_i`` is the i-th component of the + projected gradient. + eps : float or ndarray + If `jac is None` the absolute step size used for numerical + approximation of the jacobian via forward differences. + maxfun : int + Maximum number of function evaluations. + maxiter : int + Maximum number of iterations. + iprint : int, optional + Controls the frequency of output. ``iprint < 0`` means no output; + ``iprint = 0`` print only one line at the last iteration; + ``0 < iprint < 99`` print also f and ``|proj g|`` every iprint iterations; + ``iprint = 99`` print details of every iteration except n-vectors; + ``iprint = 100`` print also the changes of active set and final x; + ``iprint > 100`` print details of every iteration including x and g. + callback : callable, optional + Called after each iteration, as ``callback(xk)``, where ``xk`` is the + current parameter vector. + maxls : int, optional + Maximum number of line search steps (per iteration). Default is 20. + finite_diff_rel_step : None or array_like, optional + If `jac in ['2-point', '3-point', 'cs']` the relative step size to + use for numerical approximation of the jacobian. The absolute step + size is computed as ``h = rel_step * sign(x) * max(1, abs(x))``, + possibly adjusted to fit into the bounds. For ``method='3-point'`` + the sign of `h` is ignored. If None (default) then step is selected + automatically. + + Notes + ----- + The option `ftol` is exposed via the `scipy.optimize.minimize` interface, + but calling `scipy.optimize.fmin_l_bfgs_b` directly exposes `factr`. The + relationship between the two is ``ftol = factr * numpy.finfo(float).eps``. + I.e., `factr` multiplies the default machine floating-point precision to + arrive at `ftol`. + + """ + _check_unknown_options(unknown_options) + m = maxcor + pgtol = gtol + factr = ftol / np.finfo(float).eps + + x0 = asarray(x0).ravel() + n, = x0.shape + + if bounds is None: + bounds = [(None, None)] * n + if len(bounds) != n: + raise ValueError('length of x0 != length of bounds') + + # unbounded variables must use None, not +-inf, for optimizer to work properly + bounds = [(None if l == -np.inf else l, None if u == np.inf else u) for l, u in bounds] + # LBFGSB is sent 'old-style' bounds, 'new-style' bounds are required by + # approx_derivative and ScalarFunction + new_bounds = old_bound_to_new(bounds) + + # check bounds + if (new_bounds[0] > new_bounds[1]).any(): + raise ValueError("LBFGSB - one of the lower bounds is greater than an upper bound.") + + # initial vector must lie within the bounds. Otherwise ScalarFunction and + # approx_derivative will cause problems + x0 = np.clip(x0, new_bounds[0], new_bounds[1]) + + if disp is not None: + if disp == 0: + iprint = -1 + else: + iprint = disp + + sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps, + bounds=new_bounds, + finite_diff_rel_step=finite_diff_rel_step) + + func_and_grad = sf.fun_and_grad + + fortran_int = _lbfgsb.types.intvar.dtype + + nbd = zeros(n, fortran_int) + low_bnd = zeros(n, float64) + upper_bnd = zeros(n, float64) + bounds_map = {(None, None): 0, + (1, None): 1, + (1, 1): 2, + (None, 1): 3} + for i in range(0, n): + l, u = bounds[i] + if l is not None: + low_bnd[i] = l + l = 1 + if u is not None: + upper_bnd[i] = u + u = 1 + nbd[i] = bounds_map[l, u] + + if not maxls > 0: + raise ValueError('maxls must be positive.') + + x = array(x0, float64) + f = array(0.0, float64) + g = zeros((n,), float64) + wa = zeros(2*m*n + 5*n + 11*m*m + 8*m, float64) + iwa = zeros(3*n, fortran_int) + task = zeros(1, 'S60') + csave = zeros(1, 'S60') + lsave = zeros(4, fortran_int) + isave = zeros(44, fortran_int) + dsave = zeros(29, float64) + + task[:] = 'START' + + n_iterations = 0 + + while 1: + # x, f, g, wa, iwa, task, csave, lsave, isave, dsave = \ + _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr, + pgtol, wa, iwa, task, iprint, csave, lsave, + isave, dsave, maxls) + task_str = task.tobytes() + if task_str.startswith(b'FG'): + # The minimization routine wants f and g at the current x. + # Note that interruptions due to maxfun are postponed + # until the completion of the current minimization iteration. + # Overwrite f and g: + f, g = func_and_grad(x) + if sf.nfev > maxfun: + task[:] = ('STOP: TOTAL NO. of f AND g EVALUATIONS ' + 'EXCEEDS LIMIT') + elif task_str.startswith(b'NEW_X'): + # new iteration + n_iterations += 1 + if callback is not None: + callback(np.copy(x)) + + if n_iterations >= maxiter: + task[:] = 'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT' + elif sf.nfev > maxfun: + task[:] = ('STOP: TOTAL NO. of f AND g EVALUATIONS ' + 'EXCEEDS LIMIT') + else: + break + + task_str = task.tobytes().strip(b'\x00').strip() + if task_str.startswith(b'CONV'): + warnflag = 0 + elif sf.nfev > maxfun or n_iterations >= maxiter: + warnflag = 1 + else: + warnflag = 2 + + # These two portions of the workspace are described in the mainlb + # subroutine in lbfgsb.f. See line 363. + s = wa[0: m*n].reshape(m, n) + y = wa[m*n: 2*m*n].reshape(m, n) + + # See lbfgsb.f line 160 for this portion of the workspace. + # isave(31) = the total number of BFGS updates prior the current iteration; + n_bfgs_updates = isave[30] + + n_corrs = min(n_bfgs_updates, maxcor) + hess_inv = LbfgsInvHessProduct(s[:n_corrs], y[:n_corrs]) + + task_str = task_str.decode() + return OptimizeResult(fun=f, jac=g, nfev=sf.nfev, + njev=sf.ngev, + nit=n_iterations, status=warnflag, message=task_str, + x=x, success=(warnflag == 0), hess_inv=hess_inv) + + +class LbfgsInvHessProduct(LinearOperator): + """Linear operator for the L-BFGS approximate inverse Hessian. + + This operator computes the product of a vector with the approximate inverse + of the Hessian of the objective function, using the L-BFGS limited + memory approximation to the inverse Hessian, accumulated during the + optimization. + + Objects of this class implement the ``scipy.sparse.linalg.LinearOperator`` + interface. + + Parameters + ---------- + sk : array_like, shape=(n_corr, n) + Array of `n_corr` most recent updates to the solution vector. + (See [1]). + yk : array_like, shape=(n_corr, n) + Array of `n_corr` most recent updates to the gradient. (See [1]). + + References + ---------- + .. [1] Nocedal, Jorge. "Updating quasi-Newton matrices with limited + storage." Mathematics of computation 35.151 (1980): 773-782. + + """ + + def __init__(self, sk, yk): + """Construct the operator.""" + if sk.shape != yk.shape or sk.ndim != 2: + raise ValueError('sk and yk must have matching shape, (n_corrs, n)') + n_corrs, n = sk.shape + + super().__init__(dtype=np.float64, shape=(n, n)) + + self.sk = sk + self.yk = yk + self.n_corrs = n_corrs + self.rho = 1 / np.einsum('ij,ij->i', sk, yk) + + def _matvec(self, x): + """Efficient matrix-vector multiply with the BFGS matrices. + + This calculation is described in Section (4) of [1]. + + Parameters + ---------- + x : ndarray + An array with shape (n,) or (n,1). + + Returns + ------- + y : ndarray + The matrix-vector product + + """ + s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho + q = np.array(x, dtype=self.dtype, copy=True) + if q.ndim == 2 and q.shape[1] == 1: + q = q.reshape(-1) + + alpha = np.empty(n_corrs) + + for i in range(n_corrs-1, -1, -1): + alpha[i] = rho[i] * np.dot(s[i], q) + q = q - alpha[i]*y[i] + + r = q + for i in range(n_corrs): + beta = rho[i] * np.dot(y[i], r) + r = r + s[i] * (alpha[i] - beta) + + return r + + def todense(self): + """Return a dense array representation of this operator. + + Returns + ------- + arr : ndarray, shape=(n, n) + An array with the same shape and containing + the same data represented by this `LinearOperator`. + + """ + s, y, n_corrs, rho = self.sk, self.yk, self.n_corrs, self.rho + I = np.eye(*self.shape, dtype=self.dtype) + Hk = I + + for i in range(n_corrs): + A1 = I - s[i][:, np.newaxis] * y[i][np.newaxis, :] * rho[i] + A2 = I - y[i][:, np.newaxis] * s[i][np.newaxis, :] * rho[i] + + Hk = np.dot(A1, np.dot(Hk, A2)) + (rho[i] * s[i][:, np.newaxis] * + s[i][np.newaxis, :]) + return Hk diff --git a/src/daComposant/daAlgorithms/Atoms/mmqr.py b/src/daComposant/daAlgorithms/Atoms/mmqr.py index 48bf735..24228fe 100644 --- a/src/daComposant/daAlgorithms/Atoms/mmqr.py +++ b/src/daComposant/daAlgorithms/Atoms/mmqr.py @@ -47,7 +47,7 @@ def mmqr( Journal of Computational and Graphical Statistics, 9, 1, pp.60-77, 2000. """ # - # Recuperation des donnees et informations initiales + # Récupération des données et informations initiales # -------------------------------------------------- variables = numpy.ravel( x0 ) mesures = numpy.ravel( y ) @@ -56,7 +56,7 @@ def mmqr( n = mesures.size quantile = float(quantile) # - # Calcul des parametres du MM + # Calcul des paramètres du MM # --------------------------- tn = float(toler) / n e0 = -tn / math.log(tn) @@ -70,7 +70,7 @@ def mmqr( lastsurrogate = -numpy.sum(residus*veps) - (1.-2.*quantile)*numpy.sum(residus) iteration = 0 # - # Recherche iterative + # Recherche itérative # ------------------- while (increment > toler) and (iteration < maxfun) : iteration += 1 diff --git a/src/daComposant/daAlgorithms/Atoms/psas3dvar.py b/src/daComposant/daAlgorithms/Atoms/psas3dvar.py index 37b4eec..f69266a 100644 --- a/src/daComposant/daAlgorithms/Atoms/psas3dvar.py +++ b/src/daComposant/daAlgorithms/Atoms/psas3dvar.py @@ -120,6 +120,8 @@ def psas3dvar(selfA, Xb, Y, U, HO, CM, R, B, __storeState = False): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daAlgorithms/Atoms/std3dvar.py b/src/daComposant/daAlgorithms/Atoms/std3dvar.py index 52c7e38..8bb2288 100644 --- a/src/daComposant/daAlgorithms/Atoms/std3dvar.py +++ b/src/daComposant/daAlgorithms/Atoms/std3dvar.py @@ -123,6 +123,8 @@ def std3dvar(selfA, Xb, Y, U, HO, CM, R, B, __storeState = False): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daAlgorithms/Atoms/std4dvar.py b/src/daComposant/daAlgorithms/Atoms/std4dvar.py index c6e6504..6459918 100644 --- a/src/daComposant/daAlgorithms/Atoms/std4dvar.py +++ b/src/daComposant/daAlgorithms/Atoms/std4dvar.py @@ -184,6 +184,8 @@ def std4dvar(selfA, Xb, Y, U, HO, EM, CM, R, B, Q): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daAlgorithms/Atoms/van3dvar.py b/src/daComposant/daAlgorithms/Atoms/van3dvar.py index 587ae94..e3fc307 100644 --- a/src/daComposant/daAlgorithms/Atoms/van3dvar.py +++ b/src/daComposant/daAlgorithms/Atoms/van3dvar.py @@ -132,6 +132,8 @@ def van3dvar(selfA, Xb, Y, U, HO, CM, R, B, __storeState = False): import daAlgorithms.Atoms.lbfgsb17hlt as optimiseur elif "1.8.0" <= scipy.version.version <= "1.8.99": import daAlgorithms.Atoms.lbfgsb18hlt as optimiseur + elif "1.9.0" <= scipy.version.version <= "1.9.99": + import daAlgorithms.Atoms.lbfgsb19hlt as optimiseur else: import scipy.optimize as optimiseur Minimum, J_optimal, Informations = optimiseur.fmin_l_bfgs_b( diff --git a/src/daComposant/daCore/Aidsm.py b/src/daComposant/daCore/Aidsm.py index 99e01ef..31e18e3 100644 --- a/src/daComposant/daCore/Aidsm.py +++ b/src/daComposant/daCore/Aidsm.py @@ -864,7 +864,6 @@ class Aidsm(object): else: return self.__StoredInputs - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/AssimilationStudy.py b/src/daComposant/daCore/AssimilationStudy.py index 4af1e0a..2d28e14 100644 --- a/src/daComposant/daCore/AssimilationStudy.py +++ b/src/daComposant/daCore/AssimilationStudy.py @@ -36,7 +36,6 @@ class AssimilationStudy(_Aidsm): def __init__(self, name = ""): _Aidsm.__init__(self, name) - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/BasicObjects.py b/src/daComposant/daCore/BasicObjects.py index d5b6952..9d12eac 100644 --- a/src/daComposant/daCore/BasicObjects.py +++ b/src/daComposant/daCore/BasicObjects.py @@ -2479,7 +2479,6 @@ def MultiFonction( # logging.debug("MULTF Internal multifonction calculations end") return __multiHX - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/ExtendedLogging.py b/src/daComposant/daCore/ExtendedLogging.py index c9e31e7..2da4d7e 100644 --- a/src/daComposant/daCore/ExtendedLogging.py +++ b/src/daComposant/daCore/ExtendedLogging.py @@ -158,7 +158,6 @@ def logtimer(f): return result return wrapper - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/Interfaces.py b/src/daComposant/daCore/Interfaces.py index 9a7d54e..0876603 100644 --- a/src/daComposant/daCore/Interfaces.py +++ b/src/daComposant/daCore/Interfaces.py @@ -1256,7 +1256,6 @@ class EficasGUI(object): else: logging.debug("Can not launch standalone EFICAS/ADAO interface for path errors.") - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/NumericObjects.py b/src/daComposant/daCore/NumericObjects.py index dfeaab6..ef821f4 100644 --- a/src/daComposant/daCore/NumericObjects.py +++ b/src/daComposant/daCore/NumericObjects.py @@ -1014,7 +1014,6 @@ def multiXOsteps(selfA, Xb, Y, U, HO, EM, CM, R, B, Q, oneCycle, # return 0 - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/Persistence.py b/src/daComposant/daCore/Persistence.py index b6c1173..815b9b1 100644 --- a/src/daComposant/daCore/Persistence.py +++ b/src/daComposant/daCore/Persistence.py @@ -1017,7 +1017,6 @@ class CompositePersistence(object): # return filename - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/src/daComposant/daCore/PlatformInfo.py b/src/daComposant/daCore/PlatformInfo.py index d3f9ea5..4c53c98 100644 --- a/src/daComposant/daCore/PlatformInfo.py +++ b/src/daComposant/daCore/PlatformInfo.py @@ -219,7 +219,6 @@ class PlatformInfo(object): import daCore.version as dav return "%s %s (%s)"%(dav.name,dav.version,dav.date) - # ============================================================================== try: import scipy @@ -494,7 +493,6 @@ class SystemUsage(object): "Renvoie la mémoire totale maximale mesurée" return self._VmB('VmPeak:', unit) - # ============================================================================== if __name__ == "__main__": print('\n AUTODIAGNOSTIC\n') diff --git a/test/test1002/Performances.py b/test/test1002/Performances.py index 4f40074..e4e538e 100644 --- a/test/test1002/Performances.py +++ b/test/test1002/Performances.py @@ -41,10 +41,15 @@ class Test_Adao(unittest.TestCase): import scipy ; print(" - Scipy.............: %s"%scipy.version.version) except: print(" - Scipy.............: %s"%("absent",)) - try: - import numpy.distutils.system_info as sysinfo ; la = sysinfo.get_info('lapack') ; print(" - Lapack............: %s/lib%s.so"%(la['library_dirs'][0],la['libraries'][0])) - except: - print(" - Lapack............: %s"%("absent",)) + if tuple(map(int,numpy.version.version.split("."))) < (1,23,0): + import numpy.distutils.system_info as sysinfo + la = sysinfo.get_info('lapack') + if 'libraries' in la: + print(' - Lapack............: %s/lib%s.so'%(la['library_dirs'][0],la['libraries'][0])) + else: + print(' - Lapack............: absent') + else: + print(' - Lapack............: numpy n\'indique plus où le trouver') print("") return True @@ -102,7 +107,7 @@ class Test_Adao(unittest.TestCase): # t_init = time.time() for i in range(repetitions): - b = scipy.dot(A,x) + b = numpy.dot(A,x) print(" La duree elapsed pour %3i produits est de.......: %4.1f s"%(repetitions, time.time()-t_init)) r = [__d*(__d+1.)*(2.*__d+1.)/6.,]*__d if max(abs(b-r)) > precision: -- 2.39.2