1 # -*- coding: utf-8 -*-
3 # Copyright (C) 2008-2024 EDF R&D
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.
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.
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
19 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
23 import math, numpy, logging
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()
31 # ==============================================================================
32 class ElementaryAlgorithm(BasicObjects.Algorithm):
34 BasicObjects.Algorithm.__init__(self, "REDUCEDMODELINGTEST")
35 self.defineRequiredParameter(
36 name = "EnsembleOfSnapshots",
38 typecast = numpy.array,
39 message = "Ensemble de vecteurs d'état physique (snapshots), 1 état par colonne (Training Set)",
41 self.defineRequiredParameter(
42 name = "MaximumNumberOfLocations",
45 message = "Nombre maximal de positions",
48 self.defineRequiredParameter(
49 name = "ExcludeLocations",
52 message = "Liste des indices ou noms de positions exclues selon l'ordre interne d'un snapshot",
54 self.defineRequiredParameter(
55 name = "NameOfLocations",
58 message = "Liste des noms de positions selon l'ordre interne d'un snapshot",
60 self.defineRequiredParameter(
61 name = "SampleAsnUplet",
64 message = "Points de calcul définis par une liste de n-uplet",
66 self.defineRequiredParameter(
67 name = "SampleAsExplicitHyperCube",
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
72 self.defineRequiredParameter(
73 name = "SampleAsMinMaxStepHyperCube",
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
78 self.defineRequiredParameter(
79 name = "SampleAsMinMaxLatinHyperCube",
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
84 self.defineRequiredParameter(
85 name = "SampleAsMinMaxSobolSequence",
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
90 self.defineRequiredParameter(
91 name = "SampleAsIndependantRandomVariables",
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
96 self.defineRequiredParameter(
100 message = "Activation du mode debug lors de l'exécution",
102 self.defineRequiredParameter(
103 name = "StoreSupplementaryCalculations",
106 message = "Liste de calculs supplémentaires à stocker et/ou effectuer",
108 "EnsembleOfSimulations",
114 self.defineRequiredParameter(
116 typecast = numpy.random.seed,
117 message = "Graine fixée pour le générateur aléatoire",
119 self.defineRequiredParameter(
120 name = "MaximumNumberOfModes",
123 message = "Nombre maximal de modes pour l'analyse",
126 self.defineRequiredParameter(
127 name = "ShowElementarySummary",
130 message = "Calcule et affiche un résumé à chaque évaluation élémentaire",
132 self.defineRequiredParameter(
133 name = "NumberOfPrintedDigits",
136 message = "Nombre de chiffres affichés pour les impressions de réels",
139 self.defineRequiredParameter(
140 name = "ResultTitle",
143 message = "Titre du tableau et de la figure",
145 self.defineRequiredParameter(
147 default = self._name + "_result_file.pdf",
149 message = "Nom de base (y.c. extension) des fichiers de sauvegarde des résultats",
151 self.defineRequiredParameter(
152 name = "PlotAndSave",
155 message = "Trace et sauve les résultats",
157 self.requireInputArguments(
159 optional = ("Xb", "HO"),
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)
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)
181 raise ValueError("Snapshots or Operator have to be given in order to launch the analysis")
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)
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))
192 __fdim, __nsn = __EOS.shape
194 # --------------------------
195 __s = self._parameters["ShowElementarySummary"]
196 __p = self._parameters["NumberOfPrintedDigits"]
199 __flech = 3 * "=" + "> "
200 __ordre = int(math.log10(__nsn)) + 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")
208 msgs += (__marge + "%s\n"%self._name)
209 msgs += (__marge + "%s\n"%("=" * len(self._name),))
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")
217 msgs += (__flech + "Information before launching:\n")
218 msgs += (__marge + "-----------------------------\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
224 if "ExcludeLocations" in self._parameters:
225 __ExcludedPoints = self._parameters["ExcludeLocations"]
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"]
233 __NameOfLocations = ()
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]),
242 assume_unique = True,
245 __IncludedPoints = []
246 if len(__IncludedPoints) > 0:
247 __Ensemble = numpy.take(__EOS, __IncludedPoints, axis=0, mode='clip')
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"]]
257 if self._toStore("SingularValues"):
258 self.StoredVariables["SingularValues"].store( __sv )
259 if self._toStore("Residus"):
260 self.StoredVariables["Residus"].store( __qisv )
262 nbsv = min(5, self._parameters["MaximumNumberOfModes"])
264 msgs += (__flech + "Summary of the %i first singular values:\n"%nbsv)
265 msgs += (__marge + "---------------------------------------\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])
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])
277 msgs += (__flech + "Ordered singular values and remaining variance:\n")
278 msgs += (__marge + "-----------------------------------------------\n")
279 __entete = (" %" + str(__ordre) + "s | %22s | %22s | Variance: part, remaining")%("i", "Singular value σ", "σ[i]/σ[1]") # noqa: E501
281 __nbtirets = len(__entete) + 2
282 msgs += "\n" + __marge + "-" * __nbtirets
283 msgs += "\n" + __marge + __entete
284 msgs += "\n" + __marge + "-" * __nbtirets
287 cut1pd, cut1pc, cut1pm, cut1pi = 1, 1, 1, 1
288 for ns in range(len(__sv)):
290 rvalue = __sv[ns] / __sv[0]
291 vsinfo = 100 * __tisv[ns]
292 rsinfo = max(100 * __qisv[ns], 0.)
294 msgs += (__marge + " %0" + str(__ordre) + "i | %22." + str(__p) + "e | %22." + str(__p) + "e | %2i%s , %4.1f%s\n")%(ns, svalue, rvalue, vsinfo, "%", rsinfo, "%") # noqa: E501
296 cut1pd = ns + 2 # 10%
302 cut1pi = ns + 2 # 0.1‰
305 msgs += __marge + "-" * __nbtirets + "\n"
307 msgs += (__flech + "Summary of variance cut-off:\n")
308 msgs += (__marge + "----------------------------\n")
310 msgs += __marge + "Representing more than 90%s of variance requires at least %i mode(s).\n"%("%", cut1pd)
312 msgs += __marge + "Representing more than 99%s of variance requires at least %i mode(s).\n"%("%", cut1pc)
314 msgs += __marge + "Representing more than 99.9%s of variance requires at least %i mode(s).\n"%("%", cut1pm)
316 msgs += __marge + "Representing more than 99.99%s of variance requires at least %i mode(s).\n"%("%", cut1pi)
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)
324 msgs += (__marge + "Plot and save results in a file named \"%s\"\n"%str(self._parameters["ResultFile"]))
326 import matplotlib.pyplot as plt
327 fig = plt.figure(figsize=(10, 15))
329 if len(self._parameters["ResultTitle"]) > 0:
330 fig.suptitle(self._parameters["ResultTitle"])
332 fig.suptitle("Singular values analysis on an ensemble of %i snapshots\n"%__nsn)
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))
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}%')
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))
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")
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)
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))
365 ax.plot(range(1, 1 + len(__sv)), __sv, linewidth=2)
367 plt.savefig(str(self._parameters["ResultFile"]))
371 msgs += (__marge + "Saving figure fail, please update your Matplolib version.\n")
373 logging.getLogger().setLevel(dL)
376 msgs += (__marge + "%s\n"%("-" * 75,))
378 msgs += (__marge + "End of the \"%s\" verification\n\n"%self._name)
379 msgs += (__marge + "%s\n"%("-" * 75,))
382 self._post_run(HO, EM)
385 # ==============================================================================
386 if __name__ == "__main__":
387 print("\n AUTODIAGNOSTIC\n")