.. [Chakraborty08] Chakraborty U.K., *Advances in differential evolution*, Studies in computational intelligence, Vol.143, Springer, 2008
+.. [Chaturantabut10] Chaturantabut S., Sorensen D.C., *Nonlinear model reduction via discrete empirical interpolation*, SIMA Journal of Scientific Computing, 32(5), pp.2737-2764, 2010
+
.. [Cohn98] Cohn S. E., Da Silva A., Guo J., Sienkiewicz M., Lamich D., *Assessing the effects of data selection with the DAO Physical-space Statistical Analysis System*, Monthly Weather Review, 126, pp.2913–2926, 1998
.. [Courtier94] Courtier P., Thépaut J.-N., Hollingsworth A., *A strategy for operational implementation of 4D-Var, using an incremental approach*, Quarterly Journal of the Royal Meteorological Society, 120(519), pp.1367–1387, 1994
or of an explicit observation of the complete field(s) :math:`\mathbf{y}`.
To determine the optimum positioning of measurements, an Empirical
-Interpolation Method (EIM [Barrault04]_) is used, which establishes a reduced
-model of type Reduced Order Model (ROM), with (variant "*lcEIM*") or without
-(variant "*EIM*") positioning constraints.
+Interpolation Method (EIM [Barrault04]_) or Discrete Empirical Interpolation
+Method (DEIM [Chaturantabut10]_) is used, which establishes a reduced model of
+type Reduced Order Model (ROM), with (variant "*lcEIM*" or "*lcDEIM*") or
+without (variant "*EIM*" or "*DEIM*") positioning constraints. For performance,
+we recommend using the variant "*lcEIM*" or "*EIM*" when the dimension of the
+full fields space is large.
There are two ways to use this algorithm:
**General scheme for using the algorithm**
It is possible to exclude a priori potential positions for measurement
-positioning, using the analysis variant "*lcEIM*" for a constrained positioning
-search.
+positioning, using the analysis variant "*lcEIM*" or "*lcDEIM*" for a
+constrained positioning search.
.. ------------------------------------ ..
.. include:: snippets/Header2Algo02.rst
"OptimalPoints",
"ReducedBasis",
"Residus",
+ "SingularValues",
].
Example :
.. include:: snippets/Residus.rst
+.. include:: snippets/SingularValues.rst
+
.. ------------------------------------ ..
.. _section_ref_algorithm_MeasurementsOptimalPositioningTask_examples:
.. include:: snippets/Header2Algo07.rst
- [Barrault04]_
+- [Chaturantabut10]_
- [Gong18]_
- [Quarteroni16]_
of the particular residue checked during the running of the algorithm.
Example :
- ``rs = ADD.get("Residus")[:]``
+ ``rs = ADD.get("Residus")[-1]``
--- /dev/null
+.. index:: single: SingularValues
+
+SingularValues
+ *List of real value series*. Each element is a series, containing the
+ singular values obtained through a SVD decomposition of a collection of full
+ state vectors.
+
+ Example :
+ ``sv = ADD.get("SingularValues")[-1]``
.. index::
single: Variant
+ pair: Variant ; EIM
+ pair: Variant ; DEIM
+ pair: Variant ; lcEIM
+ pair: Variant ; lcDEIM
pair: Variant ; PositioningByEIM
+ pair: Variant ; PositioningByDEIM
pair: Variant ; PositioningBylcEIM
+ pair: Variant ; PositioningBylcDEIM
Variant
*Predefined name*. This key allows to choose one of the possible variants
for the optimal positioning search. The default variant is the constrained by
- excluded locations "PositioningBylcEIM", and the possible choices are
- "PositioningByEIM" (using the original EIM algorithm),
- "PositioningBylcEIM" (using the constrained by excluded locations EIM, named "Location Constrained EIM").
+ excluded locations "lcEIM" or "PositioningBylcEIM", and the possible choices
+ are "EIM" or "PositioningByEIM" (using the original EIM algorithm), "lcEIM"
+ or "PositioningBylcEIM" (using the constrained by excluded locations EIM,
+ named "Location Constrained EIM"), "DEIM" or "PositioningByDEIM" (using the
+ original DEIM algorithm), "lcDEIM" or "PositioningBylcDEIM" (using the
+ constrained by excluded locations DEIM, named "Location Constrained DEIM").
It is highly recommended to keep the default value.
Example :
- ``{"Variant":"PositioningBylcEIM"}``
+ ``{"Variant":"lcEIM"}``
.. [Chakraborty08] Chakraborty U.K., *Advances in differential evolution*, Studies in computational intelligence, Vol.143, Springer, 2008
+.. [Chaturantabut10] Chaturantabut S., Sorensen D.C., *Nonlinear model reduction via discrete empirical interpolation*, SIMA Journal of Scientific Computing, 32(5), pp.2737-2764, 2010
+
.. [Cohn98] Cohn S. E., Da Silva A., Guo J., Sienkiewicz M., Lamich D., *Assessing the effects of data selection with the DAO Physical-space Statistical Analysis System*, Monthly Weather Review, 126, pp.2913–2926, 1998
.. [Courtier94] Courtier P., Thépaut J.-N., Hollingsworth A., *A strategy for operational implementation of 4D-Var, using an incremental approach*, Quarterly Journal of the Royal Meteorological Society, 120(519), pp.1367–1387, 1994
du (ou des) champ(s) complet(s) :math:`\mathbf{y}`.
Pour établir la position optimale de mesures, on utilise une méthode de type
-Empirical Interpolation Method (EIM [Barrault04]_), qui établit un modèle
-réduit de type Reduced Order Model (ROM), avec contraintes (variant "*lcEIM*")
-ou sans contraintes (variant "*EIM*") de positionnement.
+Empirical Interpolation Method (EIM [Barrault04]_) ou Discrete Empirical
+Interpolation Method (DEIM [Chaturantabut10]_), qui établit un modèle réduit de
+type Reduced Order Model (ROM), avec contraintes (variante "*lcEIM*" ou
+"*lcDEIM*") ou sans contraintes (variante "*EIM*" ou "*DEIM*") de
+positionnement. Pour la performance, il est recommandé d'utiliser la variante
+"*lcEIM*" ou "*EIM*" lorsque la dimension de l'espace des champs complets est
+grande.
Il y a deux manières d'utiliser cet algorithme:
**Schéma général d'utilisation de l'algorithme**
Il est possible d'exclure a priori des positions potentielles pour le
-positionnement des mesures, en utilisant le variant "*lcEIM*" d'analyse pour
-une recherche de positionnement contraint.
+positionnement des mesures, en utilisant le variant "*lcEIM*" ou "*lcDEIM*"
+d'analyse pour une recherche de positionnement contraint.
.. ------------------------------------ ..
.. include:: snippets/Header2Algo02.rst
"OptimalPoints",
"ReducedBasis",
"Residus",
+ "SingularValues",
].
Exemple :
.. include:: snippets/Residus.rst
+.. include:: snippets/SingularValues.rst
+
.. ------------------------------------ ..
.. _section_ref_algorithm_MeasurementsOptimalPositioningTask_examples:
.. include:: snippets/Header2Algo07.rst
- [Barrault04]_
+- [Chaturantabut10]_
- [Gong18]_
- [Quarteroni16]_
l'algorithme.
Exemple :
- ``rs = ADD.get("Residus")[:]``
+ ``rs = ADD.get("Residus")[-1]``
--- /dev/null
+.. index:: single: SingularValues
+
+SingularValues
+ *Liste de série de valeurs réelles*. Chaque élément est une série, contenant
+ les valeurs singulières obtenues par une décomposition SVD d'un ensemble de
+ vecteurs d'états complets.
+
+ Exemple :
+ ``sv = ADD.get("SingularValues")[-1]``
.. index::
single: Variant
+ pair: Variant ; EIM
+ pair: Variant ; DEIM
+ pair: Variant ; lcEIM
+ pair: Variant ; lcDEIM
pair: Variant ; PositioningByEIM
+ pair: Variant ; PositioningByDEIM
pair: Variant ; PositioningBylcEIM
+ pair: Variant ; PositioningBylcDEIM
Variant
*Nom prédéfini*. Cette clé permet de choisir l'une des variantes possibles
pour la recherche du positionnement optimal. La variante par défaut est la
- version contrainte par des positions exclues "PositioningBylcEIM", et les
- choix possibles sont
- "PositioningByEIM" (utilisant l'algorithme EIM original),
- "PositioningBylcEIM" (utilisant l'algorithme EIM contraint par des positions exclues, nommé "Location Constrained EIM").
- Il est fortement recommandé de conserver la valeur par défaut.
+ version contrainte par des positions exclues "lcEIM" ou "PositioningBylcEIM",
+ et les choix possibles sont "EIM" ou "PositioningByEIM" (utilisant
+ l'algorithme EIM original), "lcEIM" ou "PositioningBylcEIM" (utilisant
+ l'algorithme EIM contraint par des positions exclues, nommé "Location
+ Constrained EIM"), "DEIM" ou "PositioningByDEIM" (utilisant l'algorithme DEIM
+ original), "lcDEIM" ou "PositioningBylcDEIM" (utilisant l'algorithme DEIM
+ contraint par des positions exclues, nommé "Location Constrained DEIM"). Il
+ est fortement recommandé de conserver la valeur par défaut.
Exemple :
- ``{"Variant":"PositioningBylcEIM"}``
+ ``{"Variant":"lcEIM"}``
--- /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__ = """
+ Empirical Interpolation Method DEIM & lcDEIM
+"""
+__author__ = "Jean-Philippe ARGAUD"
+
+import numpy, scipy, logging
+import daCore.Persistence
+from daCore.NumericObjects import FindIndexesFromNames
+
+# ==============================================================================
+def DEIM_offline(selfA, EOS = None, Verbose = False):
+ """
+ Établissement de la base
+ """
+ #
+ # Initialisations
+ # ---------------
+ 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).")
+ __dimS, __nbmS = __EOS.shape
+ logging.debug("%s Building a RB using a collection of %i snapshots of individual size of %i"%(selfA._name,__nbmS,__dimS))
+ #
+ if selfA._parameters["Variant"] in ["DEIM", "PositioningByDEIM"]:
+ __LcCsts = False
+ else:
+ __LcCsts = True
+ if __LcCsts and "ExcludeLocations" in selfA._parameters:
+ __ExcludedMagicPoints = selfA._parameters["ExcludeLocations"]
+ else:
+ __ExcludedMagicPoints = ()
+ if __LcCsts and "NameOfLocations" in selfA._parameters:
+ if isinstance(selfA._parameters["NameOfLocations"], (list, numpy.ndarray, tuple)) and len(selfA._parameters["NameOfLocations"]) == __dimS:
+ __NameOfLocations = selfA._parameters["NameOfLocations"]
+ else:
+ __NameOfLocations = ()
+ else:
+ __NameOfLocations = ()
+ if __LcCsts and len(__ExcludedMagicPoints) > 0:
+ __ExcludedMagicPoints = FindIndexesFromNames( __NameOfLocations, __ExcludedMagicPoints )
+ __ExcludedMagicPoints = numpy.ravel(numpy.asarray(__ExcludedMagicPoints, dtype=int))
+ __IncludedMagicPoints = numpy.setdiff1d(
+ numpy.arange(__EOS.shape[0]),
+ __ExcludedMagicPoints,
+ assume_unique = True,
+ )
+ else:
+ __IncludedMagicPoints = []
+ #
+ if "MaximumNumberOfLocations" in selfA._parameters and "MaximumRBSize" in selfA._parameters:
+ selfA._parameters["MaximumRBSize"] = min(selfA._parameters["MaximumNumberOfLocations"],selfA._parameters["MaximumRBSize"])
+ elif "MaximumNumberOfLocations" in selfA._parameters:
+ selfA._parameters["MaximumRBSize"] = selfA._parameters["MaximumNumberOfLocations"]
+ elif "MaximumRBSize" in selfA._parameters:
+ pass
+ else:
+ selfA._parameters["MaximumRBSize"] = __nbmS
+ __maxM = min(selfA._parameters["MaximumRBSize"], __dimS, __nbmS)
+ if "ErrorNormTolerance" in selfA._parameters:
+ selfA._parameters["EpsilonEIM"] = selfA._parameters["ErrorNormTolerance"]
+ else:
+ selfA._parameters["EpsilonEIM"] = 1.e-2
+ #
+ __U, __vs, _ = scipy.linalg.svd( __EOS )
+ __rhoM = numpy.compress(__vs > selfA._parameters["EpsilonEIM"], __U, axis=1)
+ __lVs, __svdM = __rhoM.shape
+ assert __lVs == __dimS, "Différence entre lVs et dim(EOS)"
+ __qivs = (1. - __vs[:__svdM].cumsum()/__vs.sum())
+ __maxM = min(__maxM,__svdM)
+ #
+ if __LcCsts and len(__IncludedMagicPoints) > 0:
+ __im = numpy.argmax( numpy.abs(
+ numpy.take(__rhoM[:,0], __IncludedMagicPoints, mode='clip')
+ ))
+ else:
+ __im = numpy.argmax( numpy.abs(
+ __rhoM[:,0]
+ ))
+ #
+ __mu = [None,] # Convention
+ __I = [__im,]
+ __Q = __rhoM[:,0].reshape((-1,1))
+ __errors = []
+ #
+ __M = 1 # Car le premier est déjà construit
+ __errors.append(__qivs[0])
+ #
+ # Boucle
+ # ------
+ while __M < __maxM:
+ #
+ __restrictedQi = __Q[__I,:]
+ if __M > 1:
+ __Qi_inv = numpy.linalg.inv(__restrictedQi)
+ else:
+ __Qi_inv = 1. / __restrictedQi
+ #
+ __restrictedrhoMi = __rhoM[__I,__M].reshape((-1,1))
+ #
+ if __M > 1:
+ __interpolator = numpy.dot(__Q,numpy.dot(__Qi_inv,__restrictedrhoMi))
+ else:
+ __interpolator = numpy.outer(__Q,numpy.outer(__Qi_inv,__restrictedrhoMi))
+ #
+ __residuM = __rhoM[:,__M].reshape((-1,1)) - __interpolator
+ #
+ if __LcCsts and len(__IncludedMagicPoints) > 0:
+ __im = numpy.argmax( numpy.abs(
+ numpy.take(__residuM, __IncludedMagicPoints, mode='clip')
+ ))
+ else:
+ __im = numpy.argmax( numpy.abs(
+ __residuM
+ ))
+ __Q = numpy.column_stack((__Q, __rhoM[:,__M]))
+ __I.append(__im)
+ #
+ __errors.append(__qivs[__M])
+ __mu.append(None) # Convention
+ #
+ __M = __M + 1
+ #
+ #--------------------------
+ if __errors[-1] < selfA._parameters["EpsilonEIM"]:
+ logging.debug("%s %s (%.1e)"%(selfA._name,"The convergence is obtained when reaching the required EIM tolerance",selfA._parameters["EpsilonEIM"]))
+ if __M >= __maxM:
+ logging.debug("%s %s (%i)"%(selfA._name,"The convergence is obtained when reaching the maximum number of RB dimension",__maxM))
+ logging.debug("%s The RB of size %i has been correctly build"%(selfA._name,__Q.shape[1]))
+ logging.debug("%s There are %i points that have been excluded from the potential optimal points"%(selfA._name,len(__ExcludedMagicPoints)))
+ if hasattr(selfA, "StoredVariables"):
+ selfA.StoredVariables["OptimalPoints"].store( __I )
+ if selfA._toStore("ReducedBasis"):
+ selfA.StoredVariables["ReducedBasis"].store( __Q )
+ if selfA._toStore("Residus"):
+ selfA.StoredVariables["Residus"].store( __errors )
+ if selfA._toStore("ExcludedPoints"):
+ selfA.StoredVariables["ExcludedPoints"].store( __ExcludedMagicPoints )
+ if selfA._toStore("SingularValues"):
+ selfA.StoredVariables["SingularValues"].store( __vs )
+ #
+ return __mu, __I, __Q, __errors
+
+# ==============================================================================
+# DEIM_online == EIM_online
+# ==============================================================================
+if __name__ == "__main__":
+ print('\n AUTODIAGNOSTIC\n')
#
__restrictedEOSi = __EOS[__I,:]
#
- __interpolator = numpy.empty(__EOS.shape)
if __M > 1:
__interpolator = numpy.dot(__Q,numpy.dot(__Qi_inv,__restrictedEOSi))
else:
__residuM = __dataForNextIter[:,__muM]
#
#--------------------------
- if __eM < selfA._parameters["EpsilonEIM"]:
+ if __errors[-1] < selfA._parameters["EpsilonEIM"]:
logging.debug("%s %s (%.1e)"%(selfA._name,"The convergence is obtained when reaching the required EIM tolerance",selfA._parameters["EpsilonEIM"]))
if __M >= __maxM:
logging.debug("%s %s (%i)"%(selfA._name,"The convergence is obtained when reaching the maximum number of RB dimension",__maxM))
import numpy
from daCore import BasicObjects
-from daAlgorithms.Atoms import ecweim, eosg
+from daAlgorithms.Atoms import ecweim
# ==============================================================================
class ElementaryAlgorithm(BasicObjects.Algorithm):
def __init__(self):
- # ModelInterpolationByROM
- # ModelEvaluationByReducedInterpolation
- # MeasuresInterpolationByReducedModel
- #
BasicObjects.Algorithm.__init__(self, "INTERPOLATIONBYREDUCEDMODEL")
self.defineRequiredParameter(
name = "ReducedBasis",
import numpy
from daCore import BasicObjects
-from daAlgorithms.Atoms import ecweim, eosg
+from daAlgorithms.Atoms import ecweim, ecwdeim, eosg
# ==============================================================================
class ElementaryAlgorithm(BasicObjects.Algorithm):
BasicObjects.Algorithm.__init__(self, "MEASUREMENTSOPTIMALPOSITIONING")
self.defineRequiredParameter(
name = "Variant",
- default = "PositioningBylcEIM",
+ default = "lcEIM",
typecast = str,
message = "Variant ou formulation de la méthode",
listval = [
- "EIM", "PositioningByEIM",
- "lcEIM", "PositioningBylcEIM",
+ "EIM", "PositioningByEIM",
+ "lcEIM", "PositioningBylcEIM",
+ "DEIM", "PositioningByDEIM",
+ "lcDEIM", "PositioningBylcDEIM",
],
)
self.defineRequiredParameter(
"OptimalPoints",
"ReducedBasis",
"Residus",
+ "SingularValues",
]
)
self.defineRequiredParameter(
raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
#
#--------------------------
+ elif self._parameters["Variant"] in ["lcDEIM", "PositioningBylcDEIM"]:
+ if len(self._parameters["EnsembleOfSnapshots"]) > 0:
+ if self._toStore("EnsembleOfSimulations"):
+ self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
+ ecwdeim.DEIM_offline(self, self._parameters["EnsembleOfSnapshots"])
+ elif isinstance(HO, dict):
+ ecwdeim.DEIM_offline(self, eosg.eosg(self, Xb, HO))
+ else:
+ raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
+ #
+ elif self._parameters["Variant"] in ["DEIM", "PositioningByDEIM"]:
+ if len(self._parameters["EnsembleOfSnapshots"]) > 0:
+ if self._toStore("EnsembleOfSimulations"):
+ self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
+ ecwdeim.DEIM_offline(self, self._parameters["EnsembleOfSnapshots"])
+ elif isinstance(HO, dict):
+ ecwdeim.DEIM_offline(self, eosg.eosg(self, Xb, HO))
+ else:
+ raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
+ #
+ #--------------------------
else:
raise ValueError("Error in Variant name: %s"%self._parameters["Variant"])
#
- SimulatedObservationAtCurrentState : l'état observé H(X) à l'état courant
- SimulatedObservationAtOptimum : l'état observé H(Xa) à l'optimum
- SimulationQuantiles : états observés H(X) pour les quantiles demandés
+ - SingularValues : valeurs singulières provenant d'une décomposition SVD
On peut rajouter des variables à stocker dans l'initialisation de
l'algorithme élémentaire qui va hériter de cette classe
"""
self.StoredVariables["SimulatedObservationAtCurrentState"] = Persistence.OneVector(name = "SimulatedObservationAtCurrentState")
self.StoredVariables["SimulatedObservationAtOptimum"] = Persistence.OneVector(name = "SimulatedObservationAtOptimum")
self.StoredVariables["SimulationQuantiles"] = Persistence.OneMatrix(name = "SimulationQuantiles")
+ self.StoredVariables["SingularValues"] = Persistence.OneVector(name = "SingularValues")
#
for k in self.StoredVariables:
self.__canonical_stored_name[k.lower()] = k
