Minor documentation and code review corrections (42)

diff --git a/doc/en/ref_algorithm_InterpolationByReducedModelTask.rst b/doc/en/ref_algorithm_InterpolationByReducedModelTask.rst
index 6ca9a67f54b58962c02adec274b16a0a9236e68e..989af820eedacea4906fa3371925ed2174f33ff1 100644 (file)
--- a/doc/en/ref_algorithm_InterpolationByReducedModelTask.rst
+++ b/doc/en/ref_algorithm_InterpolationByReducedModelTask.rst
@@ -49,7 +49,8 @@ field reconstruction using measurements and a reduced numerical model.
 
 To interpolate these measurements, a method of Empirical Interpolation Method
 (EIM [Barrault04]_) type is used, which uses a reduced model of type Reduced
-Order Model (ROM), with or without measurement positioning constraints.
+Order Model (ROM) issued from EIM or DEIM decomposition, with or without
+measurement positioning constraints.
 
 To use this algorithm, you need the optimally positioned measurements and the
 associated reduced basis for model representation. This can be achieved as

diff --git a/doc/fr/ref_algorithm_InterpolationByReducedModelTask.rst b/doc/fr/ref_algorithm_InterpolationByReducedModelTask.rst
index 0e188b163227d355a8c68cba59a6bd30bf1a4542..e123e8fd8ec939fef1bfd4e709844fa9ad4253bf 100644 (file)
--- a/doc/fr/ref_algorithm_InterpolationByReducedModelTask.rst
+++ b/doc/fr/ref_algorithm_InterpolationByReducedModelTask.rst
@@ -50,8 +50,8 @@ numérique réduit.
 
 Pour interpoler ces mesures, on utilise une méthode de type Empirical
 Interpolation Method (EIM [Barrault04]_), qui utilise un modèle réduit de type
-Reduced Order Model (ROM), avec ou sans contraintes de positionnement de
-mesures.
+Reduced Order Model (ROM) provenant d'une décomposition EIM or DEIM, avec ou
+sans contraintes de positionnement de mesures.
 
 Pour utiliser cet algorithme, il faut disposer des mesures optimalement
 positionnées et de la base réduite associée pour la représentation du modèle.

diff --git a/src/daComposant/daAlgorithms/Atoms/ecwdeim.py b/src/daComposant/daAlgorithms/Atoms/ecwdeim.py
index 44151778d2e24f76eba14ad66426a29c07047e12..e0b9883718f37e6e9b268c3270275aa0a687e05d 100644 (file)
--- a/src/daComposant/daAlgorithms/Atoms/ecwdeim.py
+++ b/src/daComposant/daAlgorithms/Atoms/ecwdeim.py
@@ -46,7 +46,7 @@ def DEIM_offline(selfA, EOS = None, Verbose = False):
     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))
+    logging.debug("%s Using a collection of %i snapshots of individual size of %i"%(selfA._name,__nbmS,__dimS))
     #
     if selfA._parameters["Variant"] in ["DEIM", "PositioningByDEIM"]:
         __LcCsts = False
@@ -97,20 +97,21 @@ def DEIM_offline(selfA, EOS = None, Verbose = False):
         __rhoM = scipy.linalg.orth( __EOS, selfA._parameters["EpsilonEIM"] )
     __lVs, __svdM = __rhoM.shape
     assert __lVs == __dimS, "Différence entre lVs et dim(EOS)"
-    __qivs = (1. - __vs[:__svdM].cumsum()/__vs.sum())
+    __vs2 = __vs**2
+    __qivs = 1. - __vs2[:__svdM].cumsum()/__vs2.sum()
     __maxM   = min(__maxM,__svdM)
     #
     if __LcCsts and len(__IncludedMagicPoints) > 0:
-        __im = numpy.argmax( numpy.abs(
+        __iM = numpy.argmax( numpy.abs(
             numpy.take(__rhoM[:,0], __IncludedMagicPoints, mode='clip')
             ))
     else:
-        __im = numpy.argmax( numpy.abs(
+        __iM = numpy.argmax( numpy.abs(
             __rhoM[:,0]
             ))
     #
     __mu     = [None,] # Convention
-    __I      = [__im,]
+    __I      = [__iM,]
     __Q      = __rhoM[:,0].reshape((-1,1))
     __errors = []
     #
@@ -137,17 +138,17 @@ def DEIM_offline(selfA, EOS = None, Verbose = False):
         __residuM = __rhoM[:,__M].reshape((-1,1)) - __interpolator
         #
         if __LcCsts and len(__IncludedMagicPoints) > 0:
-            __im = numpy.argmax( numpy.abs(
+            __iM = numpy.argmax( numpy.abs(
                 numpy.take(__residuM, __IncludedMagicPoints, mode='clip')
                 ))
         else:
-            __im = numpy.argmax( numpy.abs(
+            __iM = numpy.argmax( numpy.abs(
                 __residuM
                 ))
         __Q = numpy.column_stack((__Q, __rhoM[:,__M]))
-        __I.append(__im)
         #
         __errors.append(__qivs[__M])
+        __I.append(__iM)
         __mu.append(None) # Convention
         #
         __M = __M + 1

diff --git a/src/daComposant/daAlgorithms/Atoms/ecweim.py b/src/daComposant/daAlgorithms/Atoms/ecweim.py
index aab8f9ce6408e3d4a42dfc7830e0bdd5633bd704..8152493725b7ab7eb03e8354d648e69c25e947b9 100644 (file)
--- a/src/daComposant/daAlgorithms/Atoms/ecweim.py
+++ b/src/daComposant/daAlgorithms/Atoms/ecweim.py
@@ -45,7 +45,7 @@ def EIM_offline(selfA, EOS = None, Verbose = False):
     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))
+    logging.debug("%s Using a collection of %i snapshots of individual size of %i"%(selfA._name,__nbmS,__dimS))
     #
     if   selfA._parameters["ErrorNorm"] == "L2":
         MaxNormByColumn = MaxL2NormByColumn
@@ -98,7 +98,6 @@ def EIM_offline(selfA, EOS = None, Verbose = False):
     __errors = []
     #
     __M      = 0
-    __iM     = -1
     __rhoM   = numpy.empty(__dimS)
     #
     __eM, __muM = MaxNormByColumn(__EOS, __LcCsts, __IncludedMagicPoints)
@@ -131,7 +130,7 @@ def EIM_offline(selfA, EOS = None, Verbose = False):
             __Q = __rhoM.reshape((-1,1))
         __I.append(__iM)
         #
-        __restrictedQi = __Q[__I,:]
+        __restrictedQi = numpy.tril( __Q[__I,:] )
         if __M > 1:
             __Qi_inv = numpy.linalg.inv(__restrictedQi)
         else:
@@ -214,41 +213,27 @@ def EIM_online(selfA, QEIM, gJmu = None, mPoints = None, mu = None, PseudoInvers
 
 # ==============================================================================
 def MaxL2NormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
-    nmax, imax = -1, -1
     if LcCsts and len(IncludedPoints) > 0:
-        for indice in range(Ensemble.shape[1]):
-            norme = numpy.linalg.norm(
-                numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
-                )
-            if norme > nmax:
-                nmax, imax, = norme, indice
+        normes = numpy.linalg.norm(
+            numpy.take(Ensemble, IncludedPoints, axis=0, mode='clip'),
+            axis = 0,
+            )
     else:
-        for indice in range(Ensemble.shape[1]):
-            norme = numpy.linalg.norm(
-                Ensemble[:,indice],
-                )
-            if norme > nmax:
-                nmax, imax, = norme, indice
+        normes = numpy.linalg.norm( Ensemble, axis = 0)
+    nmax = numpy.max(normes)
+    imax = numpy.argmax(normes)
     return nmax, imax
 
 def MaxLinfNormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
-    nmax, imax = -1, -1
     if LcCsts and len(IncludedPoints) > 0:
-        for indice in range(Ensemble.shape[1]):
-            norme = numpy.linalg.norm(
-                numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
-                ord=numpy.inf,
-                )
-            if norme > nmax:
-                nmax, imax, = norme, indice
+        normes = numpy.linalg.norm(
+            numpy.take(Ensemble, IncludedPoints, axis=0, mode='clip'),
+            axis = 0, ord=numpy.inf,
+            )
     else:
-        for indice in range(Ensemble.shape[1]):
-            norme = numpy.linalg.norm(
-                Ensemble[:,indice],
-                ord=numpy.inf,
-                )
-            if norme > nmax:
-                nmax, imax, = norme, indice
+        normes = numpy.linalg.norm( Ensemble, axis = 0, ord=numpy.inf)
+    nmax = numpy.max(normes)
+    imax = numpy.argmax(normes)
     return nmax, imax
 
 # ==============================================================================

diff --git a/src/daComposant/daAlgorithms/InterpolationByReducedModelTest.py b/src/daComposant/daAlgorithms/InterpolationByReducedModelTest.py
index b366d0f4eb517ec19b3e7176e0d915cdc915fc58..2cfd86dfffd9105f7778f83fc12f78d23ce01181 100644 (file)
--- a/src/daComposant/daAlgorithms/InterpolationByReducedModelTest.py
+++ b/src/daComposant/daAlgorithms/InterpolationByReducedModelTest.py
@@ -140,21 +140,25 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         msgs += ("\n")
         msgs += (__marge + "%s\n\n"%("-"*75,))
         #
-        msgs += (__flech + "Interpolation error test for all given states:\n")
-        msgs += (__marge + "----------------------------------------------\n")
+        st = "Normalized interpolation error test using \"%s\" norm for all given states:"%self._parameters["ErrorNorm"]
+        msgs += (__flech + "%s\n"%st)
+        msgs += (__marge + "%s\n"%("-"*len(st),))
         msgs += ("\n")
-        Es = []
+        Ns, Es = [], []
         for ns in range(__nsn):
             __rm = __eos[__ip,ns]
             __im = ecweim.EIM_online(self, __rb, __eos[__ip,ns], __ip)
             #
             if   self._parameters["ErrorNorm"] == "L2":
-                __ecart = vfloat(numpy.linalg.norm( __eos[:,ns] - __im ) / numpy.linalg.norm( __eos[:,ns] ))
+                __norms = numpy.linalg.norm( __eos[:,ns] )
+                __ecart = vfloat(numpy.linalg.norm( __eos[:,ns] - __im ) / __norms )
             else:
-                __ecart = vfloat(numpy.linalg.norm( __eos[:,ns] - __im, ord=numpy.inf ) / numpy.linalg.norm( __eos[:,ns], ord=numpy.inf ))
+                __norms = numpy.linalg.norm( __eos[:,ns], ord=numpy.inf )
+                __ecart = vfloat(numpy.linalg.norm( __eos[:,ns] - __im, ord=numpy.inf ) / __norms )
+            Ns.append( __norms )
             Es.append( __ecart )
             if __s:
-                msgs += (__marge + "Normalized interpolation error (%s) for state number %0"+str(__ordre)+"i..: %."+str(__p)+"e\n")%(self._parameters["ErrorNorm"],ns,__ecart)
+                msgs += (__marge + "State %0"+str(__ordre)+"i: error of %."+str(__p)+"e for a state norm of %."+str(__p)+"e (= %3i%s)\n")%(ns,__ecart,__norms,100*__ecart/__norms,"%")
         msgs += ("\n")
         msgs += (__marge + "%s\n"%("-"*75,))
         #

diff --git a/src/daComposant/daAlgorithms/MeasurementsOptimalPositioningTask.py b/src/daComposant/daAlgorithms/MeasurementsOptimalPositioningTask.py
index 71989bc05baaeb03a055bb987a2a6fcc4122ba7d..129047273c2af0b4786653d1a439bfcc80ab6437 100644 (file)
--- a/src/daComposant/daAlgorithms/MeasurementsOptimalPositioningTask.py
+++ b/src/daComposant/daAlgorithms/MeasurementsOptimalPositioningTask.py
@@ -154,17 +154,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
         self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
         #
         #--------------------------
-        if   self._parameters["Variant"] in ["lcEIM", "PositioningBylcEIM"]:
-            if len(self._parameters["EnsembleOfSnapshots"]) > 0:
-                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 or Operator have to be given in order to launch the analysis")
-        #
-        elif self._parameters["Variant"] in ["EIM", "PositioningByEIM"]:
+        if   self._parameters["Variant"] in ["lcEIM", "PositioningBylcEIM", "EIM", "PositioningByEIM"]:
             if len(self._parameters["EnsembleOfSnapshots"]) > 0:
                 if self._toStore("EnsembleOfSimulations"):
                     self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
@@ -175,17 +165,7 @@ class ElementaryAlgorithm(BasicObjects.Algorithm):
                 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"]:
+        elif self._parameters["Variant"] in ["lcDEIM", "PositioningBylcDEIM", "DEIM", "PositioningByDEIM"]:
             if len(self._parameters["EnsembleOfSnapshots"]) > 0:
                 if self._toStore("EnsembleOfSimulations"):
                     self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )