src/daComposant/daAlgorithms/Atoms/ecweim.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 EIM & lcEIM
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import numpy, logging
  29 import daCore.Persistence
  30 from daCore.NumericObjects import FindIndexesFromNames
  31
  32 # ==============================================================================
  33 def EIM_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["ErrorNorm"] == "L2":
  51         MaxNormByColumn = MaxL2NormByColumn
  52     else:
  53         MaxNormByColumn = MaxLinfNormByColumn
  54     #
  55     if selfA._parameters["Variant"] in ["EIM", "PositioningByEIM"]:
  56         __LcCsts = False
  57     else:
  58         __LcCsts = True
  59     if __LcCsts and "ExcludeLocations" in selfA._parameters:
  60         __ExcludedMagicPoints = selfA._parameters["ExcludeLocations"]
  61     else:
  62         __ExcludedMagicPoints = ()
  63     if __LcCsts and "NameOfLocations" in selfA._parameters:
  64         if isinstance(selfA._parameters["NameOfLocations"], (list, numpy.ndarray, tuple)) and len(selfA._parameters["NameOfLocations"]) == __dimS:
  65             __NameOfLocations = selfA._parameters["NameOfLocations"]
  66         else:
  67             __NameOfLocations = ()
  68     else:
  69         __NameOfLocations = ()
  70     if __LcCsts and len(__ExcludedMagicPoints) > 0:
  71         __ExcludedMagicPoints = FindIndexesFromNames( __NameOfLocations, __ExcludedMagicPoints )
  72         __ExcludedMagicPoints = numpy.ravel(numpy.asarray(__ExcludedMagicPoints, dtype=int))
  73         __IncludedMagicPoints = numpy.setdiff1d(
  74             numpy.arange(__EOS.shape[0]),
  75             __ExcludedMagicPoints,
  76             assume_unique = True,
  77             )
  78     else:
  79         __IncludedMagicPoints = []
  80     #
  81     if "MaximumNumberOfLocations" in selfA._parameters and "MaximumRBSize" in selfA._parameters:
  82         selfA._parameters["MaximumRBSize"] = min(selfA._parameters["MaximumNumberOfLocations"],selfA._parameters["MaximumRBSize"])
  83     elif "MaximumNumberOfLocations" in selfA._parameters:
  84         selfA._parameters["MaximumRBSize"] = selfA._parameters["MaximumNumberOfLocations"]
  85     elif "MaximumRBSize" in selfA._parameters:
  86         pass
  87     else:
  88         selfA._parameters["MaximumRBSize"] = __nbmS
  89     __maxM   = min(selfA._parameters["MaximumRBSize"], __dimS, __nbmS)
  90     if "ErrorNormTolerance" in selfA._parameters:
  91         selfA._parameters["EpsilonEIM"] = selfA._parameters["ErrorNormTolerance"]
  92     else:
  93         selfA._parameters["EpsilonEIM"] = 1.e-2
  94     #
  95     __mu     = []
  96     __I      = []
  97     __Q      = numpy.empty(__dimS).reshape((-1,1))
  98     __errors = []
  99     #
 100     __M      = 0
 101     __iM     = -1
 102     __rhoM   = numpy.empty(__dimS)
 103     #
 104     __eM, __muM = MaxNormByColumn(__EOS, __LcCsts, __IncludedMagicPoints)
 105     __residuM = __EOS[:,__muM]
 106     __errors.append(__eM)
 107     #
 108     # Boucle
 109     # ------
 110     while __M < __maxM and __eM > selfA._parameters["EpsilonEIM"]:
 111         __M = __M + 1
 112         #
 113         __mu.append(__muM)
 114         #
 115         # Détermination du point et de la fonction magiques
 116         __abs_residuM = numpy.abs(__residuM)
 117         __iM   = numpy.argmax(__abs_residuM)
 118         __rhoM = __residuM / __residuM[__iM]
 119         #
 120         if __LcCsts and __iM in __ExcludedMagicPoints:
 121             __sIndices = numpy.argsort(__abs_residuM)
 122             __rang = -1
 123             assert __iM == __sIndices[__rang]
 124             while __iM in __ExcludedMagicPoints and __rang >= -len(__abs_residuM):
 125                 __rang = __rang - 1
 126                 __iM   = __sIndices[__rang]
 127         #
 128         if __M > 1:
 129             __Q = numpy.column_stack((__Q, __rhoM))
 130         else:
 131             __Q = __rhoM.reshape((-1,1))
 132         __I.append(__iM)
 133         #
 134         __restrictedQi = __Q[__I,:]
 135         if __M > 1:
 136             __Qi_inv = numpy.linalg.inv(__restrictedQi)
 137         else:
 138             __Qi_inv = 1. / __restrictedQi
 139         #
 140         __restrictedEOSi = __EOS[__I,:]
 141         #
 142         if __M > 1:
 143             __interpolator = numpy.dot(__Q,numpy.dot(__Qi_inv,__restrictedEOSi))
 144         else:
 145             __interpolator = numpy.outer(__Q,numpy.outer(__Qi_inv,__restrictedEOSi))
 146         #
 147         __dataForNextIter = __EOS - __interpolator
 148         __eM, __muM = MaxNormByColumn(__dataForNextIter, __LcCsts, __IncludedMagicPoints)
 149         __errors.append(__eM)
 150         #
 151         __residuM = __dataForNextIter[:,__muM]
 152     #
 153     #--------------------------
 154     if __errors[-1] < selfA._parameters["EpsilonEIM"]:
 155         logging.debug("%s %s (%.1e)"%(selfA._name,"The convergence is obtained when reaching the required EIM tolerance",selfA._parameters["EpsilonEIM"]))
 156     if __M >= __maxM:
 157         logging.debug("%s %s (%i)"%(selfA._name,"The convergence is obtained when reaching the maximum number of RB dimension",__maxM))
 158     logging.debug("%s The RB of size %i has been correctly build"%(selfA._name,__Q.shape[1]))
 159     logging.debug("%s There are %i points that have been excluded from the potential optimal points"%(selfA._name,len(__ExcludedMagicPoints)))
 160     if hasattr(selfA, "StoredVariables"):
 161         selfA.StoredVariables["OptimalPoints"].store( __I )
 162         if selfA._toStore("ReducedBasis"):
 163             selfA.StoredVariables["ReducedBasis"].store( __Q )
 164         if selfA._toStore("Residus"):
 165             selfA.StoredVariables["Residus"].store( __errors )
 166         if selfA._toStore("ExcludedPoints"):
 167             selfA.StoredVariables["ExcludedPoints"].store( __ExcludedMagicPoints )
 168     #
 169     return __mu, __I, __Q, __errors
 170
 171 # ==============================================================================
 172 def EIM_online(selfA, QEIM, gJmu = None, mPoints = None, mu = None, PseudoInverse = True, rbDimension = None, Verbose = False):
 173     """
 174     Reconstruction du champ complet
 175     """
 176     if gJmu is None and mu is None:
 177         raise ValueError("Either measurements or parameters has to be given as a list, both can not be None simultaneously.")
 178     if mPoints is None:
 179         raise ValueError("List of optimal locations for measurements has to be given.")
 180     if gJmu is not None:
 181         if len(gJmu) > len(mPoints):
 182             raise ValueError("The number of measurements (%i) has to be less or equal to the number of optimal locations (%i)."%(len(gJmu),len(mPoints)))
 183         if len(gJmu) > QEIM.shape[1]:
 184             raise ValueError("The number of measurements (%i) in optimal locations has to be less or equal to the dimension of the RB (%i)."%(len(gJmu),QEIM.shape[1]))
 185         __gJmu = numpy.ravel(gJmu)
 186     if mu is not None:
 187         # __gJmu = H(mu)
 188         raise NotImplementedError()
 189     if rbDimension is not None:
 190         rbDimension = min(QEIM.shape[1], rbDimension)
 191     else:
 192         rbDimension = QEIM.shape[1]
 193     __rbDim = min(QEIM.shape[1],len(mPoints),len(gJmu),rbDimension) # Modulation
 194     #--------------------------
 195     #
 196     # Restriction aux mesures
 197     if PseudoInverse:
 198         __QJinv = numpy.linalg.pinv( QEIM[mPoints,0:__rbDim] )
 199         __gammaMu = numpy.dot( __QJinv, __gJmu[0:__rbDim])
 200     else:
 201         __gammaMu = numpy.linalg.solve( QEIM[mPoints,0:__rbDim], __gJmu[0:__rbDim] )
 202     #
 203     # Interpolation du champ complet
 204     __gMmu = numpy.dot( QEIM[:,0:__rbDim], __gammaMu )
 205     #
 206     #--------------------------
 207     logging.debug("%s The full field of size %i has been correctly build"%(selfA._name,__gMmu.size))
 208     if hasattr(selfA, "StoredVariables"):
 209         selfA.StoredVariables["Analysis"].store( __gMmu )
 210         if selfA._toStore("ReducedCoordinates"):
 211             selfA.StoredVariables["ReducedCoordinates"].store( __gammaMu )
 212     #
 213     return __gMmu
 214
 215 # ==============================================================================
 216 def MaxL2NormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
 217     nmax, imax = -1, -1
 218     if LcCsts and len(IncludedPoints) > 0:
 219         for indice in range(Ensemble.shape[1]):
 220             norme = numpy.linalg.norm(
 221                 numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
 222                 )
 223             if norme > nmax:
 224                 nmax, imax, = norme, indice
 225     else:
 226         for indice in range(Ensemble.shape[1]):
 227             norme = numpy.linalg.norm(
 228                 Ensemble[:,indice],
 229                 )
 230             if norme > nmax:
 231                 nmax, imax, = norme, indice
 232     return nmax, imax
 233
 234 def MaxLinfNormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
 235     nmax, imax = -1, -1
 236     if LcCsts and len(IncludedPoints) > 0:
 237         for indice in range(Ensemble.shape[1]):
 238             norme = numpy.linalg.norm(
 239                 numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
 240                 ord=numpy.inf,
 241                 )
 242             if norme > nmax:
 243                 nmax, imax, = norme, indice
 244     else:
 245         for indice in range(Ensemble.shape[1]):
 246             norme = numpy.linalg.norm(
 247                 Ensemble[:,indice],
 248                 ord=numpy.inf,
 249                 )
 250             if norme > nmax:
 251                 nmax, imax, = norme, indice
 252     return nmax, imax
 253
 254 # ==============================================================================
 255 if __name__ == "__main__":
 256     print('\n AUTODIAGNOSTIC\n')