This algorithm provides optimal positioning of measurement points by an EIM
(Empirical Interpolation Method) analysis, in a iterative greedy way from a set
-of state vectors (usually called "snapshots" in reduced basis methodology).
-
-Each of these state vectors is usually (but not necessarily) the result
-:math:`\mathbf{y}` of a simulation :math:`H` for a given set of parameters
-:math:`\mathbf{x}=\mu`. In its simplest use, if the set of state vectors is
-pre-existing, it is only necessary to provide it through the algorithm options.
-
-It is also possible to exclude a priori potential locations for optimal
-measurement points, using the "*lcEIM*" analysis for a constrained positioning
+of given state vectors (usually called "snapshots" in reduced basis
+methodology) or obtained by a direct simulation during the algorithm. Each of
+these state vectors are usaully (but not necessarily) the result
+:math:`\mathbf{y}` of a simulation :math:`H` for a given set of paramters
+:math:`\mathbf{x}=\mu`.
+
+There are two ways to use this algorithm:
+
+#. In its simplest use, if the set of state vectors is pre-existing, it is only
+ necessary to provide it by the option "*SetOfSnapshots*" of algorithm.
+#. If the set of state vectors is to be obtained by simulations during the
+ course of the algorithm, then one must provide the :math:`H` simulation
+ operator and the parametric :math:`\mathbf{x}` state space design of
+ experiments.
+
+The sample of states :math:`\mathbf{x}` can be provided explicitly or in the
+form of hyper-cubes, explicit or sampled according to standard laws. Beware of
+the size of the hyper-cube (and thus the number of computations) that can be
+reached, it can quickly become important.
+
+It is possible to exclude a priori potential positions for optimal measurement
+points, using the analysis variant "*lcEIM*" for a constrained positioning
search.
.. ------------------------------------ ..
.. include:: snippets/MaximumNumberOfLocations.rst
+.. include:: snippets/SampleAsExplicitHyperCube.rst
+
+.. include:: snippets/SampleAsIndependantRandomVariables.rst
+
+.. include:: snippets/SampleAsMinMaxStepHyperCube.rst
+
+.. include:: snippets/SampleAsnUplet.rst
+
+.. include:: snippets/SetDebug.rst
+
+.. include:: snippets/SetSeed.rst
+
StoreSupplementaryCalculations
.. index:: single: StoreSupplementaryCalculations
the following part of this specific algorithmic documentation, in the
sub-section "*Information and variables available at the end of the
algorithm*"): [
- "EnsembleOfSnapshots",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"OptimalPoints",
"ReducedBasis",
"Residus",
.. ------------------------------------ ..
.. include:: snippets/Header2Algo05.rst
-.. include:: snippets/EnsembleOfSnapshots.rst
+.. include:: snippets/EnsembleOfSimulations.rst
+
+.. include:: snippets/EnsembleOfStates.rst
.. include:: snippets/OptimalPoints.rst
"CostFunctionJb",
"CostFunctionJo",
"CurrentState",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"InnovationAtCurrentState",
"SimulatedObservationAtCurrentState",
].
.. include:: snippets/CurrentState.rst
+.. include:: snippets/EnsembleOfSimulations.rst
+
+.. include:: snippets/EnsembleOfStates.rst
+
.. include:: snippets/InnovationAtCurrentState.rst
.. include:: snippets/SimulatedObservationAtCurrentState.rst
.. index:: single: setEvolutionModel
.. index:: single: setControlModel
-The operators for observation and evolution are required to implement the data
-assimilation or optimization procedures. They include the physical simulation
-by numerical calculations, but also the filtering and restriction to compare
-the simulation to observation. The evolution operator is considered here in its
-incremental form, representing the transition between two successive states,
-and is then similar to the observation operator.
-
-Schematically, an operator :math:`O` has to give a output solution for
-specified input parameters. Part of the input parameters can be modified during
-the optimization procedure. So the mathematical representation of such a
-process is a function. It was briefly described in the section
-:ref:`section_theory` and is generalized here by the relation:
+The availability of the operators of observation, and sometimes of evolution,
+are required to implement the data assimilation or optimization procedures. As
+the evolution operator is considered in its incremental form, which represents
+the transition between two successive states, it is then formally similar to
+the observation operator and the way to describe them is unique.
+
+These operators include the **physical simulation by numerical calculations**.
+But they also include **filtering, projection or restriction** of simulated
+quantities, which are necessary to compare the simulation to the observation.
+
+Schematically, an operator :math:`O` has to give a output simulation or
+solution for specified input parameters. Part of the input parameters can be
+modified during the optimization procedure. So the mathematical representation
+of such a process is a function. It was briefly described in the section
+:ref:`section_theory`. It is generalized here by the relation:
.. math:: \mathbf{y} = O( \mathbf{x} )
between the pseudo-observations outputs :math:`\mathbf{y}` and the input
-parameters :math:`\mathbf{x}` using the observation or evolution operator
-:math:`O`. The same functional representation can be used for the linear
-tangent model :math:`\mathbf{O}` of :math:`O` and its adjoint
-:math:`\mathbf{O}^*`, also required by some data assimilation or optimization
-algorithms.
+parameters :math:`\mathbf{x}` using the observation or evolution :math:`O`
+operator. The same functional representation can be used for the linear tangent
+model :math:`\mathbf{O}` of :math:`O` and its adjoint :math:`\mathbf{O}^*`,
+also required by some data assimilation or optimization algorithms.
On input and output of these operators, the :math:`\mathbf{x}` and
:math:`\mathbf{y}` variables, or their increments, are mathematically vectors,
Then, **to fully describe an operator, the user has only to provide a function
that completely and only realize the functional operation**.
-This function is usually given as a Python function or script, that can be in
-particular executed as an independent Python function or in a YACS node. These
-function or script can, with no differences, launch external codes or use
+This function is usually given as a **Python function or script**, that can be
+in particular executed as an independent Python function or in a YACS node.
+These function or script can, with no differences, launch external codes or use
internal Python or SALOME calls and methods. If the algorithm requires the 3
aspects of the operator (direct form, tangent form and adjoint form), the user
has to give the 3 functions or to approximate them using ADAO.
...
...
...
- return Y=O(X)
+ # Result: Y = O(X)
+ return "a vector similar to Y"
In this case, the user has also provide a value for the differential increment
(or keep the default value), using through the graphical interface (GUI) or
a default value of 1%. This coefficient will be used in the finite differences
approximation to build the tangent and adjoint operators. The finite
differences approximation order can also be chosen through the GUI, using the
-keyword "*CenteredFiniteDifference*", with 0 for an uncentered schema of first
-order (which is the default value), and with 1 for a centered schema of second
-order (and of twice the first order computational cost). If necessary and if
-possible, :ref:`subsection_ref_parallel_df` can be used. In all cases, an
-internal cache mechanism is used to restrict the number of operator evaluations
-at the minimum possible in a sequential or parallel execution scheme for
-numerical approximations of the tangent and adjoint operators, to avoid
-redundant calculations. One can refer to the section dealing with
-:ref:`subsection_iterative_convergence_control` to discover the interaction
-with the convergence parameters.
+keyword "*CenteredFiniteDifference*", with ``False`` or 0 for an uncentered
+schema of first order (which is the default value), and with ``True`` or 1 for
+a centered schema of second order (and of twice the first order computational
+cost). If necessary and if possible, :ref:`subsection_ref_parallel_df` can be
+used. In all cases, an internal cache mechanism is used to restrict the number
+of operator evaluations at the minimum possible in a sequential or parallel
+execution scheme for numerical approximations of the tangent and adjoint
+operators, to avoid redundant calculations. One can refer to the section
+dealing with :ref:`subsection_iterative_convergence_control` to discover the
+interaction with the convergence parameters.
This first operator definition form allows easily to test the functional form
before its use in an ADAO case, greatly reducing the complexity of operator
implementation. One can then use the "*FunctionTest*" ADAO checking algorithm
-(see the section on the :ref:`section_ref_algorithm_FunctionTest`) for this
-test.
+(see the section on the :ref:`section_ref_algorithm_FunctionTest`) specifically
+designed for this test.
+
+**Important:** the name "*DirectOperator*" is mandatory when using an
+independant Python script. The type of the input ``X`` argument can be either a
+list of float values, a Numpy array or a Numpy matrix, and the user function
+has to accept and treat all these cases. The type of the output argument ``Y``
+must also be equivalent to a list of real values.
-**Important warning:** the name "*DirectOperator*" is mandatory, and the type
-of the ``X`` argument can be either a list of float values, a Numpy array or a
-Numpy matrix. The user function has to accept and treat all these cases.
+Various forms of operators are available in several scripts included in the
+:ref:`section_docu_examples`.
.. _section_ref_operator_funcs:
else:
return "a vector similar to X"
-**Important warning:** the names "*DirectOperator*", "*TangentOperator*" and
-"*AdjointOperator*" are mandatory, and the type of the ``X``, Y``, ``dX``
-arguments can be either a list of float values, a Numpy array or a Numpy
-matrix. The user function has to treat these cases in his script.
+**Important:** the names "*DirectOperator*", "*TangentOperator*" and
+"*AdjointOperator*" are mandatory when using an independant Python script. The
+type of the ``X``, Y``, ``dX`` input or output arguments can be either a list
+of float values, a Numpy array or a Numpy matrix. The user function has to
+treat these cases in his script.
.. _section_ref_operator_switch:
--- /dev/null
+.. index:: single: EnsembleOfSimulations
+
+EnsembleOfSimulations
+ *List of vectors or matrix*. This key contains a set of physical or simulated
+ state vectors :math:`\mathbf{y}` (called "*snapshots*" in reduced basis
+ terminology), with 1 state per column if it is a matrix, or 1 state per
+ element if it is a list. Caution: the numbering of the points, to which a
+ state value is given in each vector, is implicitly that of the natural order
+ of numbering of the state vector, from 0 to the "size minus 1" of this
+ vector.
+
+ Example :
+ ``{"EnsembleOfSimulations":[y1, y2, y3...]}``
--- /dev/null
+.. index:: single: EnsembleOfStates
+
+EnsembleOfStates
+ *List of vectors or matrix*. Each element is a set of physical or parameter
+ state vectors :math:`\mathbf{x}`, with 1 state per column if it is a matrix,
+ or 1 state per element if it is a list. Caution: the numbering of the points,
+ to which a state value is given in each vector, is implicitly that of the
+ natural order of numbering of the state vector, from 0 to the "size minus 1"
+ of this vector.
+
+ Example :
+ ``{"EnsembleOfStates":[x1, x2, x3...]}``
Cet algorithme permet d'établir la position de points de mesures optimaux par
une analyse EIM (Empirical Interpolation Method), de manière itérative à partir
-d'un ensemble de vecteurs d'état (usuellement appelés "*snapshots*" en
-méthodologie de bases réduites).
-
-Chacun de ces vecteurs d'état est habituellement (mais pas obligatoirement) le
-résultat :math:`\mathbf{y}` d'une simulation :math:`H` pour un jeu de
-paramètres donné :math:`\mathbf{x}=\mu`. Dans son usage le plus simple, si
-l'ensemble des vecteurs d'état est pré-existant, il suffit de le fournir par
-les options d'algorithme.
-
-Il est aussi possible d'exclure a priori des positions potentielles pour les
-points de mesures optimaux, en utilisant l'analyse "*lcEIM*" pour une recherche
-de positionnement contraint.
+d'un ensemble de vecteurs d'état établis (usuellement appelés "*snapshots*" en
+méthodologie de bases réduites) ou obtenus par une simulation directe au cours
+de l'algorithme. Chacun de ces vecteurs d'état est habituellement (mais pas
+obligatoirement) le résultat :math:`\mathbf{y}` d'une simulation :math:`H` pour
+un jeu de paramètres donné :math:`\mathbf{x}=\mu`.
+
+Il y a deux manières d'utiliser cet algorithme:
+
+#. Dans son usage le plus simple, si l'ensemble des vecteurs d'état est
+ pré-existant, il suffit de le fournir par l'option "*EnsembleOfSnapshots*"
+ d'algorithme.
+#. Si l'ensemble des vecteurs d'état doit être obtenu par des simulations au
+ cours de l'algorithme, alors on doit fournir l'opérateur de simulation
+ :math:`H` et le plan d'expérience de l'espace des états :math:`\mathbf{x}`
+ paramétriques.
+
+L'échantillon des états :math:`\mathbf{x}` peut être fourni explicitement ou
+sous la forme d'hyper-cubes, explicites ou échantillonnés selon des lois
+courantes. Attention à la taille de l'hyper-cube (et donc au nombre de calculs)
+qu'il est possible d'atteindre, elle peut rapidement devenir importante.
+
+Il est possible d'exclure a priori des positions potentielles pour les points
+de mesures optimaux, en utilisant le variant "*lcEIM*" d'analyse pour une
+recherche de positionnement contraint.
.. ------------------------------------ ..
.. include:: snippets/Header2Algo02.rst
.. include:: snippets/MaximumNumberOfLocations.rst
+.. include:: snippets/SampleAsExplicitHyperCube.rst
+
+.. include:: snippets/SampleAsIndependantRandomVariables.rst
+
+.. include:: snippets/SampleAsMinMaxStepHyperCube.rst
+
+.. include:: snippets/SampleAsnUplet.rst
+
+.. include:: snippets/SetDebug.rst
+
+.. include:: snippets/SetSeed.rst
+
StoreSupplementaryCalculations
.. index:: single: StoreSupplementaryCalculations
(la description détaillée de chaque variable nommée est donnée dans la suite
de cette documentation par algorithme spécifique, dans la sous-partie
"*Informations et variables disponibles à la fin de l'algorithme*") : [
- "EnsembleOfSnapshots",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"OptimalPoints",
"ReducedBasis",
"Residus",
.. ------------------------------------ ..
.. include:: snippets/Header2Algo05.rst
-.. include:: snippets/EnsembleOfSnapshots.rst
+.. include:: snippets/EnsembleOfSimulations.rst
+
+.. include:: snippets/EnsembleOfStates.rst
.. include:: snippets/OptimalPoints.rst
"CostFunctionJb",
"CostFunctionJo",
"CurrentState",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"InnovationAtCurrentState",
"SimulatedObservationAtCurrentState",
].
.. include:: snippets/CurrentState.rst
+.. include:: snippets/EnsembleOfSimulations.rst
+
+.. include:: snippets/EnsembleOfStates.rst
+
.. include:: snippets/InnovationAtCurrentState.rst
.. include:: snippets/SimulatedObservationAtCurrentState.rst
.. index:: single: setEvolutionModel
.. index:: single: setControlModel
-Les opérateurs d'observation et d'évolution sont nécessaires pour mettre en
-oeuvre les procédures d'assimilation de données ou d'optimisation. Ils
-comprennent la simulation physique par des calculs numériques, mais aussi le
-filtrage et de restriction pour comparer la simulation à l'observation.
-L'opérateur d'évolution est ici considéré dans sa forme incrémentale, qui
-représente la transition entre deux états successifs, et il est alors similaire
-à l'opérateur d'observation.
-
-Schématiquement, un opérateur :math:`O` a pour objet de restituer une solution
-pour des paramètres d'entrée spécifiés. Une partie des paramètres d'entrée peut
-être modifiée au cours de la procédure d'optimisation. Ainsi, la représentation
-mathématique d'un tel processus est une fonction. Il a été brièvement décrit
-dans la section :ref:`section_theory` et il est généralisé ici par la relation:
+La disponibilité des opérateurs d'observation et parfois d'évolution sont
+nécessaires pour mettre en oeuvre les procédures d'assimilation de données ou
+d'optimisation. Comme l'opérateur d'évolution est considéré dans sa forme
+incrémentale, qui représente la transition entre deux états successifs, il est
+alors formellement similaire à l'opérateur d'observation et la manière de les
+décrire est unique.
+
+Ces opérateurs comprennent la **simulation physique par des calculs
+numériques**. Mais ils comprennent aussi **le filtrage, la projection ou la
+restriction** des grandeurs simulées, qui sont nécessaires pour comparer la
+simulation à l'observation.
+
+Schématiquement, un opérateur :math:`O` a pour objet de restituer une
+simulation ou une solution pour des paramètres d'entrée spécifiés. Une partie
+des paramètres d'entrée peut être modifiée au cours de la procédure
+d'optimisation. Ainsi, la représentation mathématique d'un tel processus est
+une fonction. Il a été brièvement décrit dans la section :ref:`section_theory`.
+Il est généralisé ici par la relation:
.. math:: \mathbf{y} = O( \mathbf{x} )
entre les pseudo-observations en sortie :math:`\mathbf{y}` et les paramètres
-d'entrée :math:`\mathbf{x}` en utilisant l'opérateur d'observation ou
-d'évolution :math:`O`. La même représentation fonctionnelle peut être utilisée
+d'entrée :math:`\mathbf{x}` en utilisant l'opérateur :math:`O` d'observation ou
+d'évolution. La même représentation fonctionnelle peut être utilisée
pour le modèle linéaire tangent :math:`\mathbf{O}` de :math:`O` et son adjoint
:math:`\mathbf{O}^*` qui sont aussi requis par certains algorithmes
d'assimilation de données ou d'optimisation.
fournir une fonction qui réalise complètement et uniquement l'opération
fonctionnelle**.
-Cette fonction est généralement donnée comme une fonction ou un script Python,
-qui peuvent en particulier être exécuté comme une fonction Python indépendante
-ou dans un noeud YACS. Cette fonction ou ce script peuvent, sans différences,
-lancer des codes externes ou utiliser des appels et des méthodes internes
-Python ou SALOME. Si l'algorithme nécessite les 3 aspects de l'opérateur (forme
-directe, forme tangente et forme adjointe), l'utilisateur doit donner les 3
-fonctions ou les approximer grâce à ADAO.
+Cette fonction est généralement donnée comme une **fonction ou un script
+Python**, qui peuvent en particulier être exécuté comme une fonction Python
+indépendante ou dans un noeud YACS. Cette fonction ou ce script peuvent, sans
+différences, lancer des codes externes ou utiliser des appels et des méthodes
+internes Python ou SALOME. Si l'algorithme nécessite les 3 aspects de
+l'opérateur (forme directe, forme tangente et forme adjointe), l'utilisateur
+doit donner les 3 fonctions ou les approximer grâce à ADAO.
Il existe pour l'utilisateur 3 méthodes effectives de fournir une représentation
fonctionnelle de l'opérateur, qui diffèrent selon le type d'argument choisi:
...
...
...
- return Y=O(X)
+ # Résultat : Y = O(X)
+ return "un vecteur similaire à Y"
Dans ce cas, l'utilisateur doit aussi fournir une valeur pour l'incrément
-différentiel (ou conserver la valeur par défaut), en utilisant dans l'interface
-graphique (GUI) ou textuelle (TUI) le mot-clé "*DifferentialIncrement*" comme
-paramètre, qui a une valeur par défaut de 1%. Ce coefficient est utilisé dans
-l'approximation différences finies pour construire les opérateurs tangent et
-adjoint. L'ordre de l'approximation différences finies peut aussi être choisi à
-travers l'interface, en utilisant le mot-clé "*CenteredFiniteDifference*", avec
-0 pour un schéma non centré du premier ordre (qui est la valeur par défaut), et
-avec 1 pour un schéma centré du second ordre (et qui coûte numériquement deux
-fois plus cher que le premier ordre). Si nécessaire et si possible, on peut
+différentiel ou conserver la valeur par défaut. Cela se réalise en utilisant
+dans l'interface graphique (GUI) ou textuelle (TUI) le mot-clé
+"*DifferentialIncrement*" comme paramètre, qui a une valeur par défaut de 1%.
+Ce coefficient est utilisé dans l'approximation différences finies pour
+construire les opérateurs tangent et adjoint. L'ordre de l'approximation
+différences finies peut aussi être choisi à travers l'interface, en utilisant
+le mot-clé "*CenteredFiniteDifference*", avec ``False`` ou 0 pour un schéma non
+centré du premier ordre (qui est la valeur par défaut), et avec ``True`` ou 1
+pour un schéma centré du second ordre (et qui coûte numériquement deux fois
+plus cher que le premier ordre). Si nécessaire et si possible, on peut
:ref:`subsection_ref_parallel_df`. Dans tous les cas, un mécanisme de cache
interne permet de limiter le nombre d'évaluations de l'opérateur au minimum
possible du point de vue de l'exécution séquentielle ou parallèle des
forme fonctionnelle avant son usage dans un cas ADAO, réduisant notablement la
complexité de l'implémentation de l'opérateur. On peut ainsi utiliser
l'algorithme ADAO de vérification "*FunctionTest*" (voir la section sur
-l':ref:`section_ref_algorithm_FunctionTest`) pour ce test.
+l':ref:`section_ref_algorithm_FunctionTest`) spécifiquement prévu pour ce test.
+
+**Important :** le nom "*DirectOperator*" est obligatoire lorsque l'on utilise
+un script Python indépendant. Le type de l'argument ``X`` en entrée peut être
+une liste de valeurs réelles, un vecteur Numpy ou une matrice Numpy, et la
+fonction utilisateur doit accepter et traiter tous ces cas. Le type de
+l'argument ``Y`` en sortie doit aussi être équivalent à une liste de valeurs
+réelles.
-**Avertissement important :** le nom "*DirectOperator*" est obligatoire, et le
-type de l'argument ``X`` peut être une liste de valeur réelles, un vecteur
-Numpy ou une matrice Numpy. La fonction utilisateur doit accepter et traiter
-tous ces cas.
+Des formes variées d'opérateurs sont disponibles dans les divers scripts inclus
+dans les :ref:`section_docu_examples`.
.. _section_ref_operator_funcs:
else:
return "un vecteur similaire à X"
-**Avertissement important :** les noms "*DirectOperator*", "*TangentOperator*"
-et "*AdjointOperator*" sont obligatoires, et le type des arguments ``X``,
-``Y``, ``dX`` peut être une liste de valeur réelles, un vecteur Numpy ou une
-matrice Numpy. La fonction utilisateur doit accepter et traiter tous ces cas.
+**Important :** les noms "*DirectOperator*", "*TangentOperator*" et
+"*AdjointOperator*" sont obligatoires lorsque l'on utilise un script Python
+indépendant. Le type des arguments en entrée ou en sortie ``X``, ``Y``, ``dX``
+peut être une liste de valeur réelles, un vecteur Numpy ou une matrice Numpy.
+La fonction utilisateur doit accepter et traiter tous ces cas.
.. _section_ref_operator_switch:
--- /dev/null
+.. index:: single: EnsembleOfSimulations
+
+EnsembleOfSimulations
+ *Liste de vecteurs ou matrice*. Chaque élément est un ensemble de vecteurs
+ d'état physique ou d'état simulé :math:`\mathbf{y}` (nommés "*snapshots*" en
+ terminologie de bases réduites), avec 1 état par colonne si c'est une
+ matrice, ou 1 état par élément si c'est une liste. Important : la
+ numérotation des points, auxquels sont fournis une valeur d'état dans chaque
+ vecteur, est implicitement celle de l'ordre naturel de numérotation du
+ vecteur d'état, de 0 à la "taille moins 1" de ce vecteur.
+
+ Exemple :
+ ``{"EnsembleOfSimulations":[y1, y2, y3...]}``
--- /dev/null
+.. index:: single: EnsembleOfStates
+
+EnsembleOfStates
+ *Liste de vecteurs ou matrice*. Chaque élément est un ensemble de vecteurs
+ d'état physique ou d'état paramétrique :math:`\mathbf{x}`, avec 1 état par
+ colonne si c'est une matrice, ou 1 état par élément si c'est une liste.
+ Important : la numérotation des points, auxquels sont fournis une valeur
+ d'état dans chaque vecteur, est implicitement celle de l'ordre naturel de
+ numérotation du vecteur d'état, de 0 à la "taille moins 1" de ce vecteur.
+
+ Exemple :
+ ``{"EnsembleOfStates":[x1, x2, x3...]}``
__author__ = "Jean-Philippe ARGAUD"
import numpy
+import daCore.Persistence
# ==============================================================================
-def EIM_offline(selfA, Verbose = False):
+def EIM_offline(selfA, EOS = None, Verbose = False):
"""
Établissement de base par Empirical Interpolation Method (EIM)
"""
#
# Initialisations
# ---------------
- if isinstance(selfA._parameters["EnsembleOfSnapshots"], (numpy.ndarray,numpy.matrix)):
- __EOS = numpy.asarray(selfA._parameters["EnsembleOfSnapshots"])
- elif isinstance(selfA._parameters["EnsembleOfSnapshots"], (list,tuple)):
- __EOS = numpy.asarray(selfA._parameters["EnsembleOfSnapshots"]).T
+ if isinstance(EOS, (numpy.ndarray, numpy.matrix)):
+ __EOS = numpy.asarray(EOS)
+ elif isinstance(EOS, (list, tuple, daCore.Persistence.Persistence)):
+ __EOS = numpy.stack([numpy.ravel(_sn) for _sn in EOS], axis=1)
+ # __EOS = numpy.asarray(EOS).T
else:
raise ValueError("EnsembleOfSnapshots has to be an array/matrix (each column being a vector) or a list/tuple (each element being a vector).")
#
--- /dev/null
+# -*- coding: utf-8 -*-
+#
+# Copyright (C) 2008-2023 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
+
+__doc__ = """
+ Ensemble Of Simulations Generation
+"""
+__author__ = "Jean-Philippe ARGAUD"
+
+import numpy, logging
+import daCore.NumericObjects
+
+# ==============================================================================
+def eosg(selfA, Xb, HO, outputEOX = False):
+ """
+ Ensemble Of Simulations Generation
+ """
+ #
+ __seed = numpy.random.get_state()
+ sampleList = daCore.NumericObjects.BuildComplexSampleList(
+ selfA._parameters["SampleAsnUplet"],
+ selfA._parameters["SampleAsExplicitHyperCube"],
+ selfA._parameters["SampleAsMinMaxStepHyperCube"],
+ selfA._parameters["SampleAsIndependantRandomVariables"],
+ Xb,
+ )
+ #
+ # ----------
+ if selfA._parameters["SetDebug"]:
+ CUR_LEVEL = logging.getLogger().getEffectiveLevel()
+ logging.getLogger().setLevel(logging.DEBUG)
+ print("===> Beginning of evaluation, activating debug\n")
+ print(" %s\n"%("-"*75,))
+ #
+ Hm = HO["Direct"].appliedTo
+ EOS = Hm(
+ sampleList,
+ argsAsSerie = True,
+ returnSerieAsArrayMatrix = True,
+ )
+ #
+ if selfA._parameters["SetDebug"]:
+ print("\n %s\n"%("-"*75,))
+ print("===> End evaluation, deactivating debug if necessary\n")
+ logging.getLogger().setLevel(CUR_LEVEL)
+ # ----------
+ #
+ if outputEOX or selfA._toStore("EnsembleOfStates"):
+ # Attention la liste s'épuise donc il faut la recréer
+ numpy.random.set_state(__seed)
+ sampleList = daCore.NumericObjects.BuildComplexSampleList(
+ selfA._parameters["SampleAsnUplet"],
+ selfA._parameters["SampleAsExplicitHyperCube"],
+ selfA._parameters["SampleAsMinMaxStepHyperCube"],
+ selfA._parameters["SampleAsIndependantRandomVariables"],
+ Xb,
+ )
+ # Il faut passer la liste en tuple/list pour stack
+ EOX = numpy.stack(tuple(sampleList), axis=1)
+ assert EOX.shape[1] == EOS.shape[1], " Error of number of states in Ensemble Of Simulations Generation"
+ if selfA._toStore("EnsembleOfStates"):
+ selfA.StoredVariables["EnsembleOfStates"].store( EOX )
+ if selfA._toStore("EnsembleOfSimulations"):
+ selfA.StoredVariables["EnsembleOfSimulations"].store( EOS )
+ #
+ if outputEOX:
+ return EOX, EOS
+ else:
+ return EOS
+
+# ==============================================================================
+if __name__ == "__main__":
+ print('\n AUTODIAGNOSTIC\n')
import numpy
from daCore import BasicObjects
from daAlgorithms.Atoms import ecweim
+from daAlgorithms.Atoms import eosg
# ==============================================================================
class ElementaryAlgorithm(BasicObjects.Algorithm):
message = "Valeur limite inférieure du critère d'optimalité forçant l'arrêt",
minval = 0.,
)
+ self.defineRequiredParameter(
+ name = "SampleAsnUplet",
+ default = [],
+ typecast = tuple,
+ message = "Points de calcul définis par une liste de n-uplet",
+ )
+ self.defineRequiredParameter(
+ name = "SampleAsExplicitHyperCube",
+ default = [],
+ typecast = tuple,
+ message = "Points de calcul définis par un hyper-cube dont on donne la liste des échantillonnages de chaque variable comme une liste",
+ )
+ self.defineRequiredParameter(
+ name = "SampleAsMinMaxStepHyperCube",
+ default = [],
+ typecast = tuple,
+ message = "Points de calcul définis par un hyper-cube dont on donne la liste des échantillonnages de chaque variable par un triplet [min,max,step]",
+ )
+ self.defineRequiredParameter(
+ name = "SampleAsIndependantRandomVariables",
+ default = [],
+ typecast = tuple,
+ message = "Points de calcul définis par un hyper-cube dont les points sur chaque axe proviennent de l'échantillonnage indépendant de la variable selon la spécification ['distribution',[parametres],nombre]",
+ )
+ self.defineRequiredParameter(
+ name = "SetDebug",
+ default = False,
+ typecast = bool,
+ message = "Activation du mode debug lors de l'exécution",
+ )
self.defineRequiredParameter(
name = "StoreSupplementaryCalculations",
default = [],
typecast = tuple,
message = "Liste de calculs supplémentaires à stocker et/ou effectuer",
listval = [
- "EnsembleOfSnapshots",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"OptimalPoints",
"ReducedBasis",
"Residus",
]
)
+ self.defineRequiredParameter(
+ name = "SetSeed",
+ typecast = numpy.random.seed,
+ message = "Graine fixée pour le générateur aléatoire",
+ )
self.requireInputArguments(
mandatory= (),
optional = ("Xb", "HO"),
#--------------------------
if self._parameters["Variant"] == "PositioningBylcEIM":
if len(self._parameters["EnsembleOfSnapshots"]) > 0:
- if self._toStore("EnsembleOfSnapshots"):
- self.StoredVariables["EnsembleOfSnapshots"].store( self._parameters["EnsembleOfSnapshots"] )
- ecweim.EIM_offline(self)
+ if self._toStore("EnsembleOfSimulations"):
+ self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
+ ecweim.EIM_offline(self, self._parameters["EnsembleOfSnapshots"])
+ elif isinstance(HO, dict):
+ ecweim.EIM_offline(self, eosg.eosg(self, Xb, HO))
else:
- raise ValueError("Snapshots have to be given in order to launch the positionning analysis")
+ raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
#
elif self._parameters["Variant"] == "PositioningByEIM":
if len(self._parameters["EnsembleOfSnapshots"]) > 0:
- if self._toStore("EnsembleOfSnapshots"):
- self.StoredVariables["EnsembleOfSnapshots"].store( self._parameters["EnsembleOfSnapshots"] )
- ecweim.EIM_offline(self)
+ if self._toStore("EnsembleOfSimulations"):
+ self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
+ ecweim.EIM_offline(self, self._parameters["EnsembleOfSnapshots"])
+ elif isinstance(HO, dict):
+ ecweim.EIM_offline(self, eosg.eosg(self, Xb, HO))
else:
- raise ValueError("Snapshots have to be given in order to launch the positionning analysis")
+ raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
#
#--------------------------
else:
#
# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
-import numpy, logging, itertools
+import numpy, logging
from daCore import BasicObjects, NumericObjects
+from daAlgorithms.Atoms import eosg
from daCore.PlatformInfo import PlatformInfo
mfp = PlatformInfo().MaximumPrecision()
"CostFunctionJb",
"CostFunctionJo",
"CurrentState",
+ "EnsembleOfSimulations",
+ "EnsembleOfStates",
"InnovationAtCurrentState",
"SimulatedObservationAtCurrentState",
]
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(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
#
- Hm = HO["Direct"].appliedTo
- #
- X0 = numpy.ravel( Xb )
- Y0 = numpy.ravel( Y )
- #
- # ---------------------------
- sampleList = NumericObjects.BuildComplexSampleList(
- self._parameters["SampleAsnUplet"],
- self._parameters["SampleAsExplicitHyperCube"],
- self._parameters["SampleAsMinMaxStepHyperCube"],
- self._parameters["SampleAsIndependantRandomVariables"],
- X0,
- )
- # ----------
+ if hasattr(Y,"store"):
+ Yb = numpy.asarray( Y[-1] ).reshape((-1,1)) # Y en Vector ou VectorSerie
+ else:
+ Yb = numpy.asarray( Y ).reshape((-1,1)) # Y en Vector ou VectorSerie
BI = B.getI()
RI = R.getI()
def CostFunction(x, HmX, QualityMeasure="AugmentedWeightedLeastSquares"):
_HX = numpy.nan
Jb, Jo, J = numpy.nan, numpy.nan, numpy.nan
else:
- _X = numpy.ravel( x )
- _HX = numpy.ravel( HmX )
+ _X = numpy.asarray( x ).reshape((-1,1))
+ _HX = numpy.asarray( HmX ).reshape((-1,1))
+ _Innovation = Yb - _HX
+ assert Yb.size == _HX.size
+ assert Yb.size == _Innovation.size
if QualityMeasure in ["AugmentedWeightedLeastSquares","AWLS","AugmentedPonderatedLeastSquares","APLS","DA"]:
if BI is None or RI is None:
raise ValueError("Background and Observation error covariance matrix has to be properly defined!")
- Jb = float( 0.5 * (_X - X0).T * (BI * (_X - X0)) )
- Jo = float( 0.5 * (Y0 - _HX).T * (RI * (Y0 - _HX)) )
+ Jb = float( 0.5 * (_X - Xb).T * (BI * (_X - Xb)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
elif QualityMeasure in ["WeightedLeastSquares","WLS","PonderatedLeastSquares","PLS"]:
if RI is None:
raise ValueError("Observation error covariance matrix has to be properly defined!")
Jb = 0.
- Jo = float( 0.5 * (Y0 - _HX).T * (RI * (Y0 - _HX)) )
+ Jo = float( 0.5 * _Innovation.T * (RI * _Innovation) )
elif QualityMeasure in ["LeastSquares","LS","L2"]:
Jb = 0.
- Jo = float( 0.5 * (Y0 - _HX).T @ (Y0 - _HX) )
+ Jo = float( 0.5 * _Innovation.T @ _Innovation )
elif QualityMeasure in ["AbsoluteValue","L1"]:
Jb = 0.
- Jo = float( numpy.sum( numpy.abs(Y0 - _HX), dtype=mfp ) )
+ Jo = float( numpy.sum( numpy.abs(_Innovation), dtype=mfp ) )
elif QualityMeasure in ["MaximumError","ME", "Linf"]:
Jb = 0.
- Jo = numpy.max( numpy.abs(Y0 - _HX) )
+ Jo = numpy.max( numpy.abs(_Innovation) )
#
J = Jb + Jo
if self._toStore("CurrentState"):
self.StoredVariables["CurrentState"].store( _X )
if self._toStore("InnovationAtCurrentState"):
- self.StoredVariables["InnovationAtCurrentState"].store( Y0 - _HX )
+ self.StoredVariables["InnovationAtCurrentState"].store( _Innovation )
if self._toStore("SimulatedObservationAtCurrentState"):
self.StoredVariables["SimulatedObservationAtCurrentState"].store( _HX )
self.StoredVariables["CostFunctionJb"].store( Jb )
self.StoredVariables["CostFunctionJo"].store( Jo )
self.StoredVariables["CostFunctionJ" ].store( J )
return J, Jb, Jo
- # ----------
- if self._parameters["SetDebug"]:
- CUR_LEVEL = logging.getLogger().getEffectiveLevel()
- logging.getLogger().setLevel(logging.DEBUG)
- print("===> Beginning of evaluation, activating debug\n")
- print(" %s\n"%("-"*75,))
#
# ----------
- for i,Xx in enumerate(sampleList):
- if self._parameters["SetDebug"]:
- print("===> Launching evaluation for state %i"%i)
- try:
- Yy = Hm( numpy.ravel( Xx ) )
- except:
- Yy = numpy.nan
- #
- J, Jb, Jo = CostFunction( Xx, Yy, self._parameters["QualityCriterion"])
- # ----------
+ EOX, EOS = eosg.eosg(self, Xb, HO, True)
#
- if self._parameters["SetDebug"]:
- print("\n %s\n"%("-"*75,))
- print("===> End evaluation, deactivating debug if necessary\n")
- logging.getLogger().setLevel(CUR_LEVEL)
+ for i in range(EOS.shape[1]):
+ J, Jb, Jo = CostFunction( EOX[:,i], EOS[:,i], self._parameters["QualityCriterion"])
+ # ----------
#
self._post_run(HO)
return 0
- CurrentOptimum : état optimal courant lors d'itérations
- CurrentState : état courant lors d'itérations
- CurrentStepNumber : pas courant d'avancement dans les algorithmes en évolution, à partir de 0
+ - EnsembleOfSimulations : ensemble d'états (sorties, simulations) rangés par colonne dans une matrice
- EnsembleOfSnapshots : ensemble d'états rangés par colonne dans une matrice
+ - EnsembleOfStates : ensemble d'états (entrées, paramètres) rangés par colonne dans une matrice
- ForecastCovariance : covariance de l'état prédit courant lors d'itérations
- ForecastState : état prédit courant lors d'itérations
- GradientOfCostFunctionJ : gradient de la fonction-coût globale
self.StoredVariables["CurrentOptimum"] = Persistence.OneVector(name = "CurrentOptimum")
self.StoredVariables["CurrentState"] = Persistence.OneVector(name = "CurrentState")
self.StoredVariables["CurrentStepNumber"] = Persistence.OneIndex(name = "CurrentStepNumber")
+ self.StoredVariables["EnsembleOfSimulations"] = Persistence.OneMatrix(name = "EnsembleOfSimulations")
self.StoredVariables["EnsembleOfSnapshots"] = Persistence.OneMatrix(name = "EnsembleOfSnapshots")
+ self.StoredVariables["EnsembleOfStates"] = Persistence.OneMatrix(name = "EnsembleOfStates")
self.StoredVariables["ForecastCovariance"] = Persistence.OneMatrix(name = "ForecastCovariance")
self.StoredVariables["ForecastState"] = Persistence.OneVector(name = "ForecastState")
self.StoredVariables["GradientOfCostFunctionJ"] = Persistence.OneVector(name = "GradientOfCostFunctionJ")
"""
try:
return numpy.asarray(self.__values).transpose()
+ # Eqvlt: return numpy.stack([numpy.ravel(sv) for sv in self.__values], axis=1)
except Exception:
raise TypeError("Base type is incompatible with numpy")
