.. include:: snippets/Header2Algo11.rst
-.. _simple_MeasurementsOptimalPositioningTask31:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask31.png
+.. _simple_MeasurementsOptimalPositioningTask3:
+.. image:: scripts/simple_MeasurementsOptimalPositioningTask3.png
:align: center
- :width: 60%
-
-.. _simple_MeasurementsOptimalPositioningTask32:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask32.png
- :align: center
- :width: 60%
-
-.. _simple_MeasurementsOptimalPositioningTask33:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask33.png
- :align: center
- :width: 60%
+ :width: 100%
.. ------------------------------------ ..
.. include:: snippets/Header2Algo06.rst
np.random.seed(123456789)
#
dimension = 20
+nbmeasures = 15
#
print("Defining a set of artificial physical fields")
from Models.TwoDimensionalInverseDistanceCS2010 \
Parameters = {
"Variant":"DEIM",
"SampleAsnUplet":Eq.get_sample_of_mu(15, 15),
- "MaximumNumberOfLocations":50,
+ "MaximumNumberOfLocations":nbmeasures,
"ErrorNorm":"Linf",
"ErrorNormTolerance":0.,
"StoreSupplementaryCalculations":[
case.setObservationOperator(OneFunction = Eq.OneRealisation)
case.execute()
#
+#-------------------------------------------------------------------------------
print()
print("Graphical display of the results")
#
sp = case.get("EnsembleOfSimulations")[-1]
+sv = case.get("SingularValues")[-1]
+op = case.get("OptimalPoints")[-1]
+#
x1, x2 = Eq.get_x()
x1, x2 = np.meshgrid(x1, x2)
-name = "Representation of a few snapshots of G out of a total of %i"%sp.shape[1]
-print(" -", name)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(6, 6))
-fig.suptitle(name)
-ax.set_xlabel("Position x1", fontweight='bold', color='red')
-ax.set_ylabel("Position x2", fontweight='bold', color='red')
-ax.set_zlabel("Magnitude of G", fontweight='bold', color='red')
-for i in range(sp.shape[1]):
- if i % 44 != 0: continue
- ax.plot_surface(x1, x2, sp[:,i].reshape((dimension, dimension)), cmap='coolwarm')
-fig.savefig("simple_MeasurementsOptimalPositioningTask31.png")
-plt.close()
+posx1 = [x1.reshape((-1,))[ip] for ip in op]
+posx2 = [x2.reshape((-1,))[ip] for ip in op]
+Omega = Eq.get_bounds_on_space()
#
-sv = case.get("SingularValues")[-1]
-name = "Singular values of the set of G simulations"
+fig = plt.figure(figsize=(18, 6))
+name = "Measurement optimal points search by "
+name += '"ADAO/MeasurementsOptimalPositioningTask/DEIM"'
+fig.suptitle(name, fontsize=20, fontstyle='italic')
+#
+ax = fig.add_subplot(1, 3, 1)
+name = "Singular values of the set of %i G simulations"%sp.shape[1]
print(" -", name)
-fig, ax = plt.subplots(figsize=(6, 6))
-fig.suptitle(name)
+ax.set_title(name, fontstyle='italic', color='red')
ax.set_xlabel("Index of singular values, numbered from 1")
ax.set_ylabel("Magnitude of singular values")
ax.set_xlim(1, len(sv))
ax.set_yscale("log")
ax.grid(True)
ax.plot(range(1, 1 + len(sv)), sv)
-fig.savefig("simple_MeasurementsOptimalPositioningTask32.png")
-plt.tight_layout()
-plt.close()
#
-nbmax = 15
-op = case.get("OptimalPoints")[-1]
-posx1 = [x1.reshape((-1,))[ip] for ip in op[:nbmax]]
-posx2 = [x2.reshape((-1,))[ip] for ip in op[:nbmax]]
-name = "Set of %i first optimal points of measurement"%nbmax
-Omega = Eq.get_bounds_on_space()
+ax = fig.add_subplot(1, 3, 2, projection = "3d")
+name = "A few simulations of G out of a total of %i"%sp.shape[1]
print(" -", name)
-fig, ax = plt.subplots(figsize=(6, 6))
-fig.suptitle(name)
-ax.set_xlabel("Position x1", fontweight='bold', color='red')
-ax.set_ylabel("Position x2", fontweight='bold', color='red')
+ax.set_title(name, fontstyle='italic', color='red')
+ax.set_xlabel("Position x1")
+ax.set_ylabel("Position x2")
+ax.set_zlabel("Magnitude of G")
+for i in range(sp.shape[1]):
+ if i % 44 != 0: continue
+ Gfield = sp[:,i].reshape((dimension, dimension))
+ ax.plot_surface(x1, x2, Gfield, cmap='coolwarm')
+#
+ax = fig.add_subplot(1, 3, 3)
+name = "Set of %i first optimal points of measurement"%nbmeasures
+print(" -", name)
+ax.set_title(name, fontstyle='italic', color='red')
+ax.set_xlabel("Position x1")
+ax.set_ylabel("Position x2")
ax.set_xlim(Omega[0][0] - 0.05, Omega[0][1] + 0.05)
ax.set_ylim(Omega[1][0] - 0.05, Omega[1][1] + 0.05)
ax.grid(True, which='both', linestyle=(0, (1, 5)), linewidth=0.5)
ax.plot(posx1, posx2, markersize=6, marker="o", linestyle='')
for i in range(len(posx1)):
ax.text(posx1[i] + 0.005, posx2[i] + 0.01, str(i + 1), fontweight='bold')
-fig.savefig("simple_MeasurementsOptimalPositioningTask33.png")
+#
plt.tight_layout()
+fig.savefig("simple_MeasurementsOptimalPositioningTask3.png")
plt.close()
Search for optimal measurement positions
Graphical display of the results
- - Representation of a few snapshots of G out of a total of 225
- - Singular values of the set of G simulations
+ - Singular values of the set of 225 G simulations
+ - A few simulations of G out of a total of 225
- Set of 15 first optimal points of measurement
The function is represented on a regular :math:`\Omega_G` spatial grid of size
20x20 points. It is available in ADAO built-in test models under the name
-`TwoDimensionalInverseDistanceCS2010`. So here we first build a set of
-simulations of :math:`G`, then apply the DEIM-type decomposition algorithm to
-it, and derive some simple illustrations.
+``TwoDimensionalInverseDistanceCS2010``, together with the spatial and
+parametric domain default definition. So here we first build a set of
+simulations of :math:`G`, then we look for the best locations for measurements
+to obtain an DEIM interpolation representation of the fields, by applying the
+DEIM-type decomposition algorithm to it, and then derive some simple
+illustrations. We choose to look for an arbitrary number ``nbmeasures`` of 15
+measurement locations.
It can be seen that the singular values decrease steadily down to numerical
noise, indicating that around a hundred basis elements are needed to fully
Python, 3.6.5, 3.12.3
Numpy, 1.14.3, 1.26.4
- Scipy, 0.19.1, 1.13.0
+ Scipy, 0.19.1, 1.13.1
MatplotLib, 2.2.2, 3.8.4
GnuplotPy, 1.8, 1.8
NLopt, 2.4.2, 2.7.1
.. include:: snippets/Header2Algo11.rst
-.. _simple_MeasurementsOptimalPositioningTask31:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask31.png
+.. _simple_MeasurementsOptimalPositioningTask3:
+.. image:: scripts/simple_MeasurementsOptimalPositioningTask3.png
:align: center
- :width: 60%
-
-.. _simple_MeasurementsOptimalPositioningTask32:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask32.png
- :align: center
- :width: 60%
-
-.. _simple_MeasurementsOptimalPositioningTask33:
-.. image:: scripts/simple_MeasurementsOptimalPositioningTask33.png
- :align: center
- :width: 60%
+ :width: 100%
.. ------------------------------------ ..
.. include:: snippets/Header2Algo06.rst
np.random.seed(123456789)
#
dimension = 20
+nbmeasures = 15
#
print("Définition d'un ensemble artificiel de champs physiques")
from Models.TwoDimensionalInverseDistanceCS2010 \
Parameters = {
"Variant":"DEIM",
"SampleAsnUplet":Eq.get_sample_of_mu(15, 15),
- "MaximumNumberOfLocations":50,
+ "MaximumNumberOfLocations":nbmeasures,
"ErrorNorm":"Linf",
"ErrorNormTolerance":0.,
"StoreSupplementaryCalculations":[
case.setObservationOperator(OneFunction = Eq.OneRealisation)
case.execute()
#
+#-------------------------------------------------------------------------------
print()
print("Affichage graphique des résultats")
#
sp = case.get("EnsembleOfSimulations")[-1]
+sv = case.get("SingularValues")[-1]
+op = case.get("OptimalPoints")[-1]
+#
x1, x2 = Eq.get_x()
x1, x2 = np.meshgrid(x1, x2)
-name = "Représentation de quelques snapshots de G sur un total de %i"%sp.shape[1]
-print(" -", name)
-fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(6, 6))
-fig.suptitle(name)
-ax.set_xlabel("Position x1", fontweight='bold', color='red')
-ax.set_ylabel("Position x2", fontweight='bold', color='red')
-ax.set_zlabel("Amplitude de G", fontweight='bold', color='red')
-for i in range(sp.shape[1]):
- if i % 44 != 0: continue
- ax.plot_surface(x1, x2, sp[:,i].reshape((dimension, dimension)), cmap='coolwarm')
-fig.savefig("simple_MeasurementsOptimalPositioningTask31.png")
-plt.close()
+posx1 = [x1.reshape((-1,))[ip] for ip in op]
+posx2 = [x2.reshape((-1,))[ip] for ip in op]
+Omega = Eq.get_bounds_on_space()
#
-sv = case.get("SingularValues")[-1]
-name = "Valeurs singulières de l'ensemble des simulations de G"
+fig = plt.figure(figsize=(18, 6))
+name = "Recherche de points optimaux de mesure par "
+name += '"ADAO/MeasurementsOptimalPositioningTask/DEIM"'
+fig.suptitle(name, fontsize=20, fontstyle='italic')
+#
+ax = fig.add_subplot(1, 3, 1)
+name = "Valeurs singulières de l'ensemble des %i simulations de G"%sp.shape[1]
print(" -", name)
-fig, ax = plt.subplots(figsize=(6, 6))
-fig.suptitle(name)
+ax.set_title(name, fontstyle='italic', color='red')
ax.set_xlabel("Index des valeurs singulières, numérotées à partir de 1")
ax.set_ylabel("Amplitude des valeurs singulières")
ax.set_xlim(1, len(sv))
ax.set_yscale("log")
ax.grid(True)
ax.plot(range(1, 1 + len(sv)), sv)
-fig.savefig("simple_MeasurementsOptimalPositioningTask32.png")
-plt.tight_layout()
-plt.close()
#
-nbmax = 15
-op = case.get("OptimalPoints")[-1]
-posx1 = [x1.reshape((-1,))[ip] for ip in op[:nbmax]]
-posx2 = [x2.reshape((-1,))[ip] for ip in op[:nbmax]]
-name = "Ensemble des %i premiers points optimaux de mesure"%nbmax
-Omega = Eq.get_bounds_on_space()
+ax = fig.add_subplot(1, 3, 2, projection = "3d")
+name = "Quelques simulations de G sur un total de %i"%sp.shape[1]
print(" -", name)
-fig, ax = plt.subplots(figsize=(6, 6))
-fig.suptitle(name)
-ax.set_xlabel("Position x1", fontweight='bold', color='red')
-ax.set_ylabel("Position x2", fontweight='bold', color='red')
+ax.set_title(name, fontstyle='italic', color='red')
+ax.set_xlabel("Position x1")
+ax.set_ylabel("Position x2")
+ax.set_zlabel("Amplitude de G")
+for i in range(sp.shape[1]):
+ if i % 44 != 0: continue
+ Gfield = sp[:,i].reshape((dimension, dimension))
+ ax.plot_surface(x1, x2, Gfield, cmap='coolwarm')
+#
+ax = fig.add_subplot(1, 3, 3)
+name = "Ensemble des %i premiers points optimaux de mesure"%nbmeasures
+print(" -", name)
+ax.set_title(name, fontstyle='italic', color='red')
+ax.set_xlabel("Position x1")
+ax.set_ylabel("Position x2")
ax.set_xlim(Omega[0][0] - 0.05, Omega[0][1] + 0.05)
ax.set_ylim(Omega[1][0] - 0.05, Omega[1][1] + 0.05)
ax.grid(True, which='both', linestyle=(0, (1, 5)), linewidth=0.5)
ax.plot(posx1, posx2, markersize=6, marker="o", linestyle='')
for i in range(len(posx1)):
ax.text(posx1[i] + 0.005, posx2[i] + 0.01, str(i + 1), fontweight='bold')
-fig.savefig("simple_MeasurementsOptimalPositioningTask33.png")
+#
plt.tight_layout()
+fig.savefig("simple_MeasurementsOptimalPositioningTask3.png")
plt.close()
Recherche des positions optimales de mesure
Affichage graphique des résultats
- - Représentation de quelques snapshots de G sur un total de 225
- - Valeurs singulières de l'ensemble des simulations de G
+ - Valeurs singulières de l'ensemble des 225 simulations de G
+ - Quelques simulations de G sur un total de 225
- Ensemble des 15 premiers points optimaux de mesure
La fonction est représenté sur une grille spatiale régulière :math:`\Omega_G`
de taille 20x20 points. Elle est disponible dans les modèles de tests intégrés
-pour ADAO sous le nom `TwoDimensionalInverseDistanceCS2010`. On construit donc
-ici tout d'abord un ensemble de simulations de :math:`G`, pour lui appliquer
-ensuite l'algorithme de décomposition de type DEIM, et en tirer des
-illustrations simples.
+pour ADAO sous le nom ``TwoDimensionalInverseDistanceCS2010``. On construit
+donc ici tout d'abord un ensemble de simulations de :math:`G`, puis on cherche
+les meilleurs positions de mesures pour obtenir une représentation par
+interpolation DEIM des champs, en appliquant l'algorithme de décomposition de
+type DEIM, et on en tire ensuite des illustrations simples. On choisit de
+rechercher un nombre arbitraire ``nbmeasures`` de 15 positions de mesures.
On observe ainsi que les valeurs singulières décroissent régulièrement jusqu'au
bruit numérique, indiquant qu'il faut environ une centaine d'éléments de base
Python, 3.6.5, 3.12.3
Numpy, 1.14.3, 1.26.4
- Scipy, 0.19.1, 1.13.0
+ Scipy, 0.19.1, 1.13.1
MatplotLib, 2.2.2, 3.8.4
GnuplotPy, 1.8, 1.8
NLopt, 2.4.2, 2.7.1
SIAM Journal on Scientific Computing, 32(5), pp. 2737-2764 (2010).
"""
def __init__(self, nx: int = 20, ny: int = 20):
+ "Définition du maillage spatial"
self.nx = max(1, nx)
self.ny = max(1, ny)
self.x = numpy.linspace(0.1, 0.9, self.nx, dtype=float)
self.y = numpy.linspace(0.1, 0.9, self.ny, dtype=float)
def FunctionH(self, XX ):
+ "Fonction simulation pour un paramètre donné"
__mu1, __mu2 = numpy.ravel(XX)
#
__x, __y = numpy.meshgrid( self.x, self.y )
return __sxymu
def get_x(self):
+ "Renvoie le maillage spatial"
return self.x, self.y
def get_sample_of_mu(self, ns1: int = 20, ns2: int = 20):
+ "Renvoie l'échantillonnage paramétrique régulier"
smu1 = numpy.linspace(-1, -0.01, ns1, dtype=float)
smu2 = numpy.linspace(-1, -0.01, ns2, dtype=float)
smu = numpy.array([(mu1, mu2) for mu1 in smu1 for mu2 in smu2])
return smu
def get_random_sample_of_mu(self, ns1: int = 1, ns2: int = 1):
+ "Renvoie l'échantillonnage paramétrique aléatoire"
smu = []
for i in range(ns1 * ns2):
smu1 = numpy.random.uniform(-1, -0.01)
return smu
def get_bounds_on_space(self):
+ "Renvoie les bornes sur le maillage spatial"
return [[min(self.x), max(self.x)], [min(self.y), max(self.y)]]
def get_bounds_on_parameter(self):
+ "Renvoie les bornes sur le maillage paramétrique"
return [[-1, -0.01]] * 2
OneRealisation = FunctionH
