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     EIM
  25 """
  26 __author__ = "Jean-Philippe ARGAUD"
  27
  28 import numpy
  29 import daCore.Persistence
  30
  31 # ==============================================================================
  32 def EIM_offline(selfA, EOS = None, Verbose = False):
  33     """
  34     Établissement de base par Empirical Interpolation Method (EIM)
  35     """
  36     #
  37     # Initialisations
  38     # ---------------
  39     if isinstance(EOS, (numpy.ndarray, numpy.matrix)):
  40         __EOS = numpy.asarray(EOS)
  41     elif isinstance(EOS, (list, tuple, daCore.Persistence.Persistence)):
  42         __EOS = numpy.stack([numpy.ravel(_sn) for _sn in EOS], axis=1)
  43         # __EOS = numpy.asarray(EOS).T
  44     else:
  45         raise ValueError("EnsembleOfSnapshots has to be an array/matrix (each column being a vector) or a list/tuple (each element being a vector).")
  46     #
  47     if   selfA._parameters["ErrorNorm"] == "L2":
  48         MaxNormByColumn = MaxL2NormByColumn
  49     else:
  50         MaxNormByColumn = MaxLinfNormByColumn
  51     #
  52     if selfA._parameters["Variant"] == "PositioningByEIM":
  53         __LcCsts = False
  54     else:
  55         __LcCsts = True
  56     if __LcCsts and "ExcludeLocations" in selfA._parameters:
  57         __ExcludedMagicPoints = selfA._parameters["ExcludeLocations"]
  58     else:
  59         __ExcludedMagicPoints = []
  60     if __LcCsts and len(__ExcludedMagicPoints) > 0:
  61         __ExcludedMagicPoints = numpy.ravel(numpy.asarray(__ExcludedMagicPoints, dtype=int))
  62         __IncludedMagicPoints = numpy.setdiff1d(
  63             numpy.arange(__EOS.shape[0]),
  64             __ExcludedMagicPoints,
  65             assume_unique = True,
  66             )
  67     else:
  68         __IncludedMagicPoints = []
  69     #
  70     __dimS, __nbmS = __EOS.shape
  71     if "MaximumNumberOfLocations" in selfA._parameters and "MaximumRBSize" in selfA._parameters:
  72         selfA._parameters["MaximumRBSize"] = min(selfA._parameters["MaximumNumberOfLocations"],selfA._parameters["MaximumRBSize"])
  73     elif "MaximumNumberOfLocations" in selfA._parameters:
  74         selfA._parameters["MaximumRBSize"] = selfA._parameters["MaximumNumberOfLocations"]
  75     elif "MaximumRBSize" in selfA._parameters:
  76         pass
  77     else:
  78         selfA._parameters["MaximumRBSize"] = __nbmS
  79     __maxM   = min(selfA._parameters["MaximumRBSize"], __dimS, __nbmS)
  80     if "ErrorNormTolerance" in selfA._parameters:
  81         selfA._parameters["EpsilonEIM"] = selfA._parameters["ErrorNormTolerance"]
  82     else:
  83         selfA._parameters["EpsilonEIM"] = 1.e-2
  84     #
  85     __mu     = []
  86     __I      = []
  87     __Q      = numpy.empty(__dimS)
  88     __errors = []
  89     #
  90     __M      = 0
  91     __iM     = -1
  92     __rhoM   = numpy.empty(__dimS)
  93     #
  94     __eM, __muM = MaxNormByColumn(__EOS, __LcCsts, __IncludedMagicPoints)
  95     __residuM = __EOS[:,__muM]
  96     __errors.append(__eM)
  97     #
  98     # Boucle
  99     # ------
 100     while __M < __maxM and __eM > selfA._parameters["EpsilonEIM"]:
 101         __M = __M + 1
 102         #
 103         __mu.append(__muM)
 104         #
 105         # Détermination du point et de la fonction magiques
 106         __abs_residuM = numpy.abs(__residuM)
 107         __iM   = numpy.argmax(__abs_residuM)
 108         __rhoM = __residuM / __abs_residuM[__iM]
 109         #
 110         if __LcCsts and __iM in __ExcludedMagicPoints:
 111             __sIndices = numpy.argsort(__abs_residuM)
 112             __rang = -1
 113             assert __iM == __sIndices[__rang]
 114             while __iM in __ExcludedMagicPoints and __rang >= -len(__abs_residuM):
 115                 __rang = __rang - 1
 116                 __iM   = __sIndices[__rang]
 117         #
 118         if __M > 1:
 119             __Q = numpy.column_stack((__Q, __rhoM))
 120         else:
 121             __Q = __rhoM
 122         __I.append(__iM)
 123         #
 124         __restrictedQi = __Q[__I]
 125         if __M > 1:
 126             __Qi_inv = numpy.linalg.inv(__restrictedQi)
 127         else:
 128             __Qi_inv = 1. / __restrictedQi
 129         #
 130         __restrictedEOSi = __EOS[__I]
 131         #
 132         __interpolator = numpy.empty(__EOS.shape)
 133         if __M > 1:
 134             __interpolator = numpy.dot(__Q,numpy.dot(__Qi_inv,__restrictedEOSi))
 135         else:
 136             __interpolator = numpy.outer(__Q,numpy.outer(__Qi_inv,__restrictedEOSi))
 137         #
 138         __dataForNextIter = __EOS - __interpolator
 139         __eM, __muM = MaxNormByColumn(__dataForNextIter, __LcCsts, __IncludedMagicPoints)
 140         __errors.append(__eM)
 141         #
 142         __residuM = __dataForNextIter[:,__muM]
 143     #
 144     #--------------------------
 145     if hasattr(selfA, "StoredVariables"):
 146         selfA.StoredVariables["OptimalPoints"].store( __I )
 147         if selfA._toStore("ReducedBasis"):
 148             selfA.StoredVariables["ReducedBasis"].store( __Q )
 149         if selfA._toStore("Residus"):
 150             selfA.StoredVariables["Residus"].store( __errors )
 151     #
 152     return __mu, __I, __Q, __errors
 153
 154 # ==============================================================================
 155 def EIM_online(selfA, QEIM, mu, iEIM):
 156     raise NotImplementedError()
 157
 158 # ==============================================================================
 159 def MaxL2NormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
 160     nmax, imax = -1, -1
 161     if LcCsts and len(IncludedPoints) > 0:
 162         for indice in range(Ensemble.shape[1]):
 163             norme = numpy.linalg.norm(
 164                 numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
 165                 )
 166             if norme > nmax:
 167                 nmax, imax, = norme, indice
 168     else:
 169         for indice in range(Ensemble.shape[1]):
 170             norme = numpy.linalg.norm(
 171                 Ensemble[:,indice],
 172                 )
 173             if norme > nmax:
 174                 nmax, imax, = norme, indice
 175     return nmax, imax
 176
 177 def MaxLinfNormByColumn(Ensemble, LcCsts = False, IncludedPoints = []):
 178     nmax, imax = -1, -1
 179     if LcCsts and len(IncludedPoints) > 0:
 180         for indice in range(Ensemble.shape[1]):
 181             norme = numpy.linalg.norm(
 182                 numpy.take(Ensemble[:,indice], IncludedPoints, mode='clip'),
 183                 ord=numpy.inf,
 184                 )
 185             if norme > nmax:
 186                 nmax, imax, = norme, indice
 187     else:
 188         for indice in range(Ensemble.shape[1]):
 189             norme = numpy.linalg.norm(
 190                 Ensemble[:,indice],
 191                 ord=numpy.inf,
 192                 )
 193             if norme > nmax:
 194                 nmax, imax, = norme, indice
 195     return nmax, imax
 196
 197 # ==============================================================================
 198 if __name__ == "__main__":
 199     print('\n AUTODIAGNOSTIC\n')