src/daComposant/daAlgorithms/ReducedModelingTest.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 # Copyright (C) 2008-2024 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 import math, numpy, logging
  24 import daCore
  25 from daCore import BasicObjects, PlatformInfo
  26 from daCore.NumericObjects import FindIndexesFromNames, SingularValuesEstimation
  27 from daAlgorithms.Atoms import eosg
  28 mpr = PlatformInfo.PlatformInfo().MachinePrecision()
  29 mfp = PlatformInfo.PlatformInfo().MaximumPrecision()
  30
  31 # ==============================================================================
  32 class ElementaryAlgorithm(BasicObjects.Algorithm):
  33     def __init__(self):
  34         BasicObjects.Algorithm.__init__(self, "REDUCEDMODELINGTEST")
  35         self.defineRequiredParameter(
  36             name     = "EnsembleOfSnapshots",
  37             default  = [],
  38             typecast = numpy.array,
  39             message  = "Ensemble de vecteurs d'état physique (snapshots), 1 état par colonne (Training Set)",
  40         )
  41         self.defineRequiredParameter(
  42             name     = "MaximumNumberOfLocations",
  43             default  = 1,
  44             typecast = int,
  45             message  = "Nombre maximal de positions",
  46             minval   = 0,
  47         )
  48         self.defineRequiredParameter(
  49             name     = "ExcludeLocations",
  50             default  = [],
  51             typecast = tuple,
  52             message  = "Liste des indices ou noms de positions exclues selon l'ordre interne d'un snapshot",
  53         )
  54         self.defineRequiredParameter(
  55             name     = "NameOfLocations",
  56             default  = [],
  57             typecast = tuple,
  58             message  = "Liste des noms de positions selon l'ordre interne d'un snapshot",
  59         )
  60         self.defineRequiredParameter(
  61             name     = "SampleAsnUplet",
  62             default  = [],
  63             typecast = tuple,
  64             message  = "Points de calcul définis par une liste de n-uplet",
  65         )
  66         self.defineRequiredParameter(
  67             name     = "SampleAsExplicitHyperCube",
  68             default  = [],
  69             typecast = tuple,
  70             message  = "Points de calcul définis par un hyper-cube dont on donne la liste des échantillonnages explicites de chaque variable comme une liste",  # noqa: E501
  71         )
  72         self.defineRequiredParameter(
  73             name     = "SampleAsMinMaxStepHyperCube",
  74             default  = [],
  75             typecast = tuple,
  76             message  = "Points de calcul définis par un hyper-cube dont on donne la liste des échantillonnages implicites de chaque variable par un triplet [min,max,step]",  # noqa: E501
  77         )
  78         self.defineRequiredParameter(
  79             name     = "SampleAsMinMaxLatinHyperCube",
  80             default  = [],
  81             typecast = tuple,
  82             message  = "Points de calcul définis par un hyper-cube Latin dont on donne les bornes de chaque variable par une paire [min,max], suivi du nombre de points demandés",  # noqa: E501
  83         )
  84         self.defineRequiredParameter(
  85             name     = "SampleAsMinMaxSobolSequence",
  86             default  = [],
  87             typecast = tuple,
  88             message  = "Points de calcul définis par une séquence de Sobol dont on donne les bornes de chaque variable par une paire [min,max], suivi de la paire [dimension, nombre minimal de points demandés]",  # noqa: E501
  89         )
  90         self.defineRequiredParameter(
  91             name     = "SampleAsIndependantRandomVariables",
  92             default  = [],
  93             typecast = tuple,
  94             message  = "Points de calcul définis par un hyper-cube dont les points sur chaque axe proviennent de l'échantillonnage indépendant de la variable selon la spécification ['distribution',[parametres],nombre]",  # noqa: E501
  95         )
  96         self.defineRequiredParameter(
  97             name     = "SetDebug",
  98             default  = False,
  99             typecast = bool,
 100             message  = "Activation du mode debug lors de l'exécution",
 101         )
 102         self.defineRequiredParameter(
 103             name     = "StoreSupplementaryCalculations",
 104             default  = [],
 105             typecast = tuple,
 106             message  = "Liste de calculs supplémentaires à stocker et/ou effectuer",
 107             listval  = [
 108                 "EnsembleOfSimulations",
 109                 "EnsembleOfStates",
 110                 "Residus",
 111                 "SingularValues",
 112             ]
 113         )
 114         self.defineRequiredParameter(
 115             name     = "SetSeed",
 116             typecast = numpy.random.seed,
 117             message  = "Graine fixée pour le générateur aléatoire",
 118         )
 119         self.defineRequiredParameter(
 120             name     = "MaximumNumberOfModes",
 121             default  = 1000000,
 122             typecast = int,
 123             message  = "Nombre maximal de modes pour l'analyse",
 124             minval   = 0,
 125         )
 126         self.defineRequiredParameter(
 127             name     = "ShowElementarySummary",
 128             default  = True,
 129             typecast = bool,
 130             message  = "Calcule et affiche un résumé à chaque évaluation élémentaire",
 131         )
 132         self.defineRequiredParameter(
 133             name     = "NumberOfPrintedDigits",
 134             default  = 5,
 135             typecast = int,
 136             message  = "Nombre de chiffres affichés pour les impressions de réels",
 137             minval   = 0,
 138         )
 139         self.defineRequiredParameter(
 140             name     = "ResultTitle",
 141             default  = "",
 142             typecast = str,
 143             message  = "Titre du tableau et de la figure",
 144         )
 145         self.defineRequiredParameter(
 146             name     = "ResultFile",
 147             default  = self._name + "_result_file.pdf",
 148             typecast = str,
 149             message  = "Nom de base (y.c. extension) des fichiers de sauvegarde des résultats",
 150         )
 151         self.defineRequiredParameter(
 152             name     = "PlotAndSave",
 153             default  = False,
 154             typecast = bool,
 155             message  = "Trace et sauve les résultats",
 156         )
 157         self.requireInputArguments(
 158             mandatory= (),
 159             optional = ("Xb", "HO"),
 160         )
 161         self.setAttributes(
 162             tags=(
 163                 "Reduction",
 164                 "Checking",
 165             ),
 166             features=(
 167                 "DerivativeFree",
 168             ),
 169         )
 170
 171     def run(self, Xb=None, Y=None, U=None, HO=None, EM=None, CM=None, R=None, B=None, Q=None, Parameters=None):
 172         self._pre_run(Parameters, Xb, Y, U, HO, EM, CM, R, B, Q)
 173         #
 174         if len(self._parameters["EnsembleOfSnapshots"]) > 0:
 175             if self._toStore("EnsembleOfSimulations"):
 176                 self.StoredVariables["EnsembleOfSimulations"].store( self._parameters["EnsembleOfSnapshots"] )
 177             EOS = self._parameters["EnsembleOfSnapshots"]
 178         elif isinstance(HO, dict):
 179             EOS = eosg.eosg(self, Xb, HO)
 180         else:
 181             raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
 182         #
 183         if isinstance(EOS, (numpy.ndarray, numpy.matrix)):
 184             __EOS = numpy.asarray(EOS)
 185         elif isinstance(EOS, (list, tuple, daCore.Persistence.Persistence)):
 186             __EOS = numpy.stack([numpy.ravel(_sn) for _sn in EOS], axis=1)
 187         else:
 188             raise ValueError("EnsembleOfSnapshots has to be an array/matrix (each column being a vector) or a list/tuple (each element being a vector).")  # noqa: E501
 189         __dimS, __nbmS = __EOS.shape
 190         logging.debug("%s Using a collection of %i snapshots of individual size of %i"%(self._name, __nbmS, __dimS))
 191         #
 192         __fdim, __nsn = __EOS.shape
 193         #
 194         # --------------------------
 195         __s = self._parameters["ShowElementarySummary"]
 196         __p = self._parameters["NumberOfPrintedDigits"]
 197         #
 198         __marge = 5 * u" "
 199         __flech = 3 * "=" + "> "
 200         __ordre = int(math.log10(__nsn)) + 1
 201         msgs  = ("\n")  # 1
 202         if len(self._parameters["ResultTitle"]) > 0:
 203             __rt = str(self._parameters["ResultTitle"])
 204             msgs += (__marge + "====" + "=" * len(__rt) + "====\n")
 205             msgs += (__marge + "    " + __rt + "\n")
 206             msgs += (__marge + "====" + "=" * len(__rt) + "====\n")
 207         else:
 208             msgs += (__marge + "%s\n"%self._name)
 209             msgs += (__marge + "%s\n"%("=" * len(self._name),))
 210         #
 211         msgs += ("\n")
 212         msgs += (__marge + "This test allows to analyze the characteristics of the collection of\n")
 213         msgs += (__marge + "states from a reduction point of view. Using an SVD, it measures how\n")
 214         msgs += (__marge + "the information decreases with the number of singular values, either\n")
 215         msgs += (__marge + "as values or, with a statistical point of view, as remaining variance.\n")
 216         msgs += ("\n")
 217         msgs += (__flech + "Information before launching:\n")
 218         msgs += (__marge + "-----------------------------\n")
 219         msgs += ("\n")
 220         msgs += (__marge + "Characteristics of input data:\n")
 221         msgs += (__marge + "  State dimension................: %i\n")%__fdim
 222         msgs += (__marge + "  Number of snapshots to test....: %i\n")%__nsn
 223         #
 224         if "ExcludeLocations" in self._parameters:
 225             __ExcludedPoints = self._parameters["ExcludeLocations"]
 226         else:
 227             __ExcludedPoints = ()
 228         if "NameOfLocations" in self._parameters:
 229             if isinstance(self._parameters["NameOfLocations"], (list, numpy.ndarray, tuple)) \
 230                     and len(self._parameters["NameOfLocations"]) == __dimS:
 231                 __NameOfLocations = self._parameters["NameOfLocations"]
 232             else:
 233                 __NameOfLocations = ()
 234         else:
 235             __NameOfLocations = ()
 236         if len(__ExcludedPoints) > 0:
 237             __ExcludedPoints = FindIndexesFromNames( __NameOfLocations, __ExcludedPoints )
 238             __ExcludedPoints = numpy.ravel(numpy.asarray(__ExcludedPoints, dtype=int))
 239             __IncludedPoints = numpy.setdiff1d(
 240                 numpy.arange(__EOS.shape[0]),
 241                 __ExcludedPoints,
 242                 assume_unique = True,
 243             )
 244         else:
 245             __IncludedPoints = []
 246         if len(__IncludedPoints) > 0:
 247             __Ensemble = numpy.take(__EOS, __IncludedPoints, axis=0, mode='clip')
 248         else:
 249             __Ensemble = __EOS
 250         #
 251         __sv, __svsq, __tisv, __qisv = SingularValuesEstimation( __Ensemble )
 252         if self._parameters["MaximumNumberOfModes"] < len(__sv):
 253             __sv   = __sv[:self._parameters["MaximumNumberOfModes"]]
 254             __tisv = __tisv[:self._parameters["MaximumNumberOfModes"]]
 255             __qisv = __qisv[:self._parameters["MaximumNumberOfModes"]]
 256         #
 257         if self._toStore("SingularValues"):
 258             self.StoredVariables["SingularValues"].store( __sv )
 259         if self._toStore("Residus"):
 260             self.StoredVariables["Residus"].store( __qisv )
 261         #
 262         nbsv = min(5, self._parameters["MaximumNumberOfModes"])
 263         msgs += ("\n")
 264         msgs += (__flech + "Summary of the %i first singular values:\n"%nbsv)
 265         msgs += (__marge + "---------------------------------------\n")
 266         msgs += ("\n")
 267         msgs += (__marge + "Singular values σ:\n")
 268         for i in range(nbsv):
 269             msgs += __marge + ("  σ[%i] = %." + str(__p) + "e\n")%(i + 1, __sv[i])
 270         msgs += ("\n")
 271         msgs += (__marge + "Singular values σ divided by the first one σ[1]:\n")
 272         for i in range(nbsv):
 273             msgs += __marge + ("  σ[%i] / σ[1] = %." + str(__p) + "e\n")%(i + 1, __sv[i] / __sv[0])
 274         #
 275         if __s:
 276             msgs += ("\n")
 277             msgs += (__flech + "Ordered singular values and remaining variance:\n")
 278             msgs += (__marge + "-----------------------------------------------\n")
 279             __entete = ("  %" + str(__ordre) + "s  | %16s | %16s | Variance: part, remaining")%("i", "Singular value σ", "σ[i]/σ[1]")  # noqa: E501
 280             #
 281             __nbtirets = len(__entete) + 2
 282             msgs += "\n" + __marge + "-" * __nbtirets
 283             msgs += "\n" + __marge + __entete
 284             msgs += "\n" + __marge + "-" * __nbtirets
 285             msgs += ("\n")
 286         #
 287         cut1pd, cut1pc, cut1pm, cut1pi = 1, 1, 1, 1
 288         for ns in range(len(__sv)):
 289             svalue = __sv[ns]
 290             rvalue = __sv[ns] / __sv[0]
 291             vsinfo = 100 * __tisv[ns]
 292             rsinfo = max(100 * __qisv[ns], 0.)
 293             if __s:
 294                 msgs += (__marge + "  %0" + str(__ordre) + "i  | %16." + str(__p) + "e | %16." + str(__p) + "e |           %2i%s ,    %4.1f%s\n")%(ns, svalue, rvalue, vsinfo, "%", rsinfo, "%")  # noqa: E501
 295             if rsinfo > 10:
 296                 cut1pd = ns + 2  # 10%
 297             if rsinfo > 1:
 298                 cut1pc = ns + 2  # 1%
 299             if rsinfo > 0.1:
 300                 cut1pm = ns + 2  # 1‰
 301             if rsinfo > 0.01:
 302                 cut1pi = ns + 2  # 0.1‰
 303         #
 304         if __s:
 305             msgs += __marge + "-" * __nbtirets + "\n"
 306         msgs += ("\n")
 307         msgs += (__flech + "Summary of variance cut-off:\n")
 308         msgs += (__marge + "----------------------------\n")
 309         if cut1pd > 0:
 310             msgs += __marge + "Representing more than 90%s    of variance requires at least %i mode(s).\n"%("%", cut1pd)
 311         if cut1pc > 0:
 312             msgs += __marge + "Representing more than 99%s    of variance requires at least %i mode(s).\n"%("%", cut1pc)
 313         if cut1pm > 0:
 314             msgs += __marge + "Representing more than 99.9%s  of variance requires at least %i mode(s).\n"%("%", cut1pm)
 315         if cut1pi > 0:
 316             msgs += __marge + "Representing more than 99.99%s of variance requires at least %i mode(s).\n"%("%", cut1pi)
 317         #
 318         if PlatformInfo.has_matplotlib and self._parameters["PlotAndSave"]:
 319             # Evite les message debug de matplotlib
 320             dL = logging.getLogger().getEffectiveLevel()
 321             logging.getLogger().setLevel(logging.WARNING)
 322             try:
 323                 msgs += ("\n")
 324                 msgs += (__marge + "Plot and save results in a file named \"%s\"\n"%str(self._parameters["ResultFile"]))
 325                 #
 326                 import matplotlib.pyplot as plt
 327                 fig = plt.figure(figsize=(10, 15))
 328                 plt.tight_layout()
 329                 if len(self._parameters["ResultTitle"]) > 0:
 330                     fig.suptitle(self._parameters["ResultTitle"])
 331                 else:
 332                     fig.suptitle("Singular values analysis on an ensemble of %i snapshots\n"%__nsn)
 333                 # ----
 334                 ax = fig.add_subplot(3, 1, 1)
 335                 ax.set_xlabel("Singular values index, numbered from 1 (first %i ones)"%len(__qisv))
 336                 ax.set_ylabel("Remaining variance to be explained (%, linear scale)", color="tab:blue")
 337                 ax.grid(True, which='both', color="tab:blue")
 338                 ax.set_xlim(1, 1 + len(__qisv))
 339                 ax.set_ylim(0, 100)
 340                 ax.plot(range(1, 1 + len(__qisv)), 100 * __qisv, linewidth=2, color="b", label="On linear scale")
 341                 ax.tick_params(axis='y', labelcolor="tab:blue")
 342                 ax.yaxis.set_major_formatter('{x:.0f}%')
 343                 #
 344                 rg = ax.twinx()
 345                 rg.set_ylabel("Remaining variance to be explained (%, log scale)", color="tab:red")
 346                 rg.grid(True, which='both', color="tab:red")
 347                 rg.set_xlim(1, 1 + len(__qisv))
 348                 rg.set_yscale("log")
 349                 rg.plot(range(1, 1 + len(__qisv)), 100 * __qisv, linewidth=2, color="r", label="On log10 scale")
 350                 rg.set_ylim(rg.get_ylim()[0], 101)
 351                 rg.tick_params(axis='y', labelcolor="tab:red")
 352                 # ----
 353                 ax = fig.add_subplot(3, 1, 2)
 354                 ax.set_ylabel("Singular values")
 355                 ax.set_xlim(1, 1 + len(__sv))
 356                 ax.plot(range(1, 1 + len(__sv)), __sv, linewidth=2)
 357                 ax.grid(True)
 358                 # ----
 359                 ax = fig.add_subplot(3, 1, 3)
 360                 ax.set_ylabel("Singular values (log scale)")
 361                 ax.grid(True, which='both')
 362                 ax.set_xlim(1, 1 + len(__sv))
 363                 ax.set_xscale("log")
 364                 ax.set_yscale("log")
 365                 ax.plot(range(1, 1 + len(__sv)), __sv, linewidth=2)
 366                 # ----
 367                 plt.savefig(str(self._parameters["ResultFile"]))
 368                 plt.close(fig)
 369             except Exception:
 370                 msgs += ("\n")
 371                 msgs += (__marge + "Saving figure fail, please update your Matplolib version.\n")
 372                 msgs += ("\n")
 373             logging.getLogger().setLevel(dL)
 374             #
 375         msgs += ("\n")
 376         msgs += (__marge + "%s\n"%("-" * 75,))
 377         msgs += ("\n")
 378         msgs += (__marge + "End of the \"%s\" verification\n\n"%self._name)
 379         msgs += (__marge + "%s\n"%("-" * 75,))
 380         print(msgs)  # 3
 381         #
 382         self._post_run(HO, EM)
 383         return 0
 384
 385 # ==============================================================================
 386 if __name__ == "__main__":
 387     print("\n AUTODIAGNOSTIC\n")