src/daComposant/daAlgorithms/Atoms/ecwdeim.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2023 EDF R&D
   4 #
   5 # This library is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation; either
   8 # version 2.1 of the License.
   9 #
  10 # This library is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 # Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with this library; if not, write to the Free Software
  17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  18 #
  19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  20 #
  21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
  22
  23 __doc__ = """
  24     Empirical Interpolation Method DEIM & lcDEIM
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import numpy, scipy, logging
  29 import daCore.Persistence
  30 from daCore.NumericObjects import FindIndexesFromNames
  31
  32 # ==============================================================================
  33 def DEIM_offline(selfA, EOS = None, Verbose = False):
  34     """
  35     Établissement de la base
  36     """
  37     #
  38     # Initialisations
  39     # ---------------
  40     if isinstance(EOS, (numpy.ndarray, numpy.matrix)):
  41         __EOS = numpy.asarray(EOS)
  42     elif isinstance(EOS, (list, tuple, daCore.Persistence.Persistence)):
  43         __EOS = numpy.stack([numpy.ravel(_sn) for _sn in EOS], axis=1)
  44         # __EOS = numpy.asarray(EOS).T
  45     else:
  46         raise ValueError("EnsembleOfSnapshots has to be an array/matrix (each column being a vector) or a list/tuple (each element being a vector).")
  47     __dimS, __nbmS = __EOS.shape
  48     logging.debug("%s Building a RB using a collection of %i snapshots of individual size of %i"%(selfA._name,__nbmS,__dimS))
  49     #
  50     if selfA._parameters["Variant"] in ["DEIM", "PositioningByDEIM"]:
  51         __LcCsts = False
  52     else:
  53         __LcCsts = True
  54     if __LcCsts and "ExcludeLocations" in selfA._parameters:
  55         __ExcludedMagicPoints = selfA._parameters["ExcludeLocations"]
  56     else:
  57         __ExcludedMagicPoints = ()
  58     if __LcCsts and "NameOfLocations" in selfA._parameters:
  59         if isinstance(selfA._parameters["NameOfLocations"], (list, numpy.ndarray, tuple)) and len(selfA._parameters["NameOfLocations"]) == __dimS:
  60             __NameOfLocations = selfA._parameters["NameOfLocations"]
  61         else:
  62             __NameOfLocations = ()
  63     else:
  64         __NameOfLocations = ()
  65     if __LcCsts and len(__ExcludedMagicPoints) > 0:
  66         __ExcludedMagicPoints = FindIndexesFromNames( __NameOfLocations, __ExcludedMagicPoints )
  67         __ExcludedMagicPoints = numpy.ravel(numpy.asarray(__ExcludedMagicPoints, dtype=int))
  68         __IncludedMagicPoints = numpy.setdiff1d(
  69             numpy.arange(__EOS.shape[0]),
  70             __ExcludedMagicPoints,
  71             assume_unique = True,
  72             )
  73     else:
  74         __IncludedMagicPoints = []
  75     #
  76     if "MaximumNumberOfLocations" in selfA._parameters and "MaximumRBSize" in selfA._parameters:
  77         selfA._parameters["MaximumRBSize"] = min(selfA._parameters["MaximumNumberOfLocations"],selfA._parameters["MaximumRBSize"])
  78     elif "MaximumNumberOfLocations" in selfA._parameters:
  79         selfA._parameters["MaximumRBSize"] = selfA._parameters["MaximumNumberOfLocations"]
  80     elif "MaximumRBSize" in selfA._parameters:
  81         pass
  82     else:
  83         selfA._parameters["MaximumRBSize"] = __nbmS
  84     __maxM   = min(selfA._parameters["MaximumRBSize"], __dimS, __nbmS)
  85     if "ErrorNormTolerance" in selfA._parameters:
  86         selfA._parameters["EpsilonEIM"] = selfA._parameters["ErrorNormTolerance"]
  87     else:
  88         selfA._parameters["EpsilonEIM"] = 1.e-2
  89     #
  90     __U, __vs, _ = scipy.linalg.svd( __EOS )
  91     __rhoM = numpy.compress(__vs > selfA._parameters["EpsilonEIM"], __U, axis=1)
  92     __lVs, __svdM = __rhoM.shape
  93     assert __lVs == __dimS, "Différence entre lVs et dim(EOS)"
  94     __qivs = (1. - __vs[:__svdM].cumsum()/__vs.sum())
  95     __maxM   = min(__maxM,__svdM)
  96     #
  97     if __LcCsts and len(__IncludedMagicPoints) > 0:
  98         __im = numpy.argmax( numpy.abs(
  99             numpy.take(__rhoM[:,0], __IncludedMagicPoints, mode='clip')
 100             ))
 101     else:
 102         __im = numpy.argmax( numpy.abs(
 103             __rhoM[:,0]
 104             ))
 105     #
 106     __mu     = [None,] # Convention
 107     __I      = [__im,]
 108     __Q      = __rhoM[:,0].reshape((-1,1))
 109     __errors = []
 110     #
 111     __M      = 1 # Car le premier est déjà construit
 112     __errors.append(__qivs[0])
 113     #
 114     # Boucle
 115     # ------
 116     while __M < __maxM:
 117         #
 118         __restrictedQi = __Q[__I,:]
 119         if __M > 1:
 120             __Qi_inv = numpy.linalg.inv(__restrictedQi)
 121         else:
 122             __Qi_inv = 1. / __restrictedQi
 123         #
 124         __restrictedrhoMi = __rhoM[__I,__M].reshape((-1,1))
 125         #
 126         if __M > 1:
 127             __interpolator = numpy.dot(__Q,numpy.dot(__Qi_inv,__restrictedrhoMi))
 128         else:
 129             __interpolator = numpy.outer(__Q,numpy.outer(__Qi_inv,__restrictedrhoMi))
 130         #
 131         __residuM = __rhoM[:,__M].reshape((-1,1)) - __interpolator
 132         #
 133         if __LcCsts and len(__IncludedMagicPoints) > 0:
 134             __im = numpy.argmax( numpy.abs(
 135                 numpy.take(__residuM, __IncludedMagicPoints, mode='clip')
 136                 ))
 137         else:
 138             __im = numpy.argmax( numpy.abs(
 139                 __residuM
 140                 ))
 141         __Q = numpy.column_stack((__Q, __rhoM[:,__M]))
 142         __I.append(__im)
 143         #
 144         __errors.append(__qivs[__M])
 145         __mu.append(None) # Convention
 146         #
 147         __M = __M + 1
 148     #
 149     #--------------------------
 150     if __errors[-1] < selfA._parameters["EpsilonEIM"]:
 151         logging.debug("%s %s (%.1e)"%(selfA._name,"The convergence is obtained when reaching the required EIM tolerance",selfA._parameters["EpsilonEIM"]))
 152     if __M >= __maxM:
 153         logging.debug("%s %s (%i)"%(selfA._name,"The convergence is obtained when reaching the maximum number of RB dimension",__maxM))
 154     logging.debug("%s The RB of size %i has been correctly build"%(selfA._name,__Q.shape[1]))
 155     logging.debug("%s There are %i points that have been excluded from the potential optimal points"%(selfA._name,len(__ExcludedMagicPoints)))
 156     if hasattr(selfA, "StoredVariables"):
 157         selfA.StoredVariables["OptimalPoints"].store( __I )
 158         if selfA._toStore("ReducedBasis"):
 159             selfA.StoredVariables["ReducedBasis"].store( __Q )
 160         if selfA._toStore("Residus"):
 161             selfA.StoredVariables["Residus"].store( __errors )
 162         if selfA._toStore("ExcludedPoints"):
 163             selfA.StoredVariables["ExcludedPoints"].store( __ExcludedMagicPoints )
 164         if selfA._toStore("SingularValues"):
 165             selfA.StoredVariables["SingularValues"].store( __vs )
 166     #
 167     return __mu, __I, __Q, __errors
 168
 169 # ==============================================================================
 170 # DEIM_online == EIM_online
 171 # ==============================================================================
 172 if __name__ == "__main__":
 173     print('\n AUTODIAGNOSTIC\n')