doc/en/scripts/simple_MeasurementsOptimalPositioningTask3.py

   1 # -*- coding: utf-8 -*-
   2 #
   3 import numpy as np
   4 import matplotlib.pyplot as plt
   5 np.random.seed(123456789)
   6 #
   7 dimension = 20
   8 nbmeasures = 15
   9 #
  10 print("Defining a set of artificial physical fields")
  11 from Models.TwoDimensionalInverseDistanceCS2010 \
  12      import TwoDimensionalInverseDistanceCS2010 as Equation
  13 Eq = Equation(dimension, dimension)
  14 print()
  15 #
  16 print("Search for optimal measurement positions")
  17 from adao import adaoBuilder
  18 case = adaoBuilder.New()
  19 case.setAlgorithmParameters(
  20     Algorithm = 'MeasurementsOptimalPositioningTask',
  21     Parameters = {
  22         "Variant":"DEIM",
  23         "SampleAsnUplet":Eq.get_sample_of_mu(15, 15),
  24         "MaximumNumberOfLocations":nbmeasures,
  25         "ErrorNorm":"Linf",
  26         "ErrorNormTolerance":0.,
  27         "StoreSupplementaryCalculations":[
  28             "SingularValues",
  29             "OptimalPoints",
  30             "EnsembleOfSimulations",
  31             ],
  32         }
  33     )
  34 case.setBackground(Vector = [1,1] )
  35 case.setObservationOperator(OneFunction = Eq.OneRealisation)
  36 case.execute()
  37 #
  38 #-------------------------------------------------------------------------------
  39 print()
  40 print("Graphical display of the results")
  41 #
  42 sp = case.get("EnsembleOfSimulations")[-1]
  43 sv = case.get("SingularValues")[-1]
  44 op = case.get("OptimalPoints")[-1]
  45 #
  46 x1, x2 = Eq.get_x()
  47 x1, x2 = np.meshgrid(x1, x2)
  48 posx1 = [x1.reshape((-1,))[ip] for ip in op]
  49 posx2 = [x2.reshape((-1,))[ip] for ip in op]
  50 Omega = Eq.get_bounds_on_space()
  51 #
  52 fig = plt.figure(figsize=(18, 6))
  53 name = "Measurement optimal points search by "
  54 name += '"ADAO/MeasurementsOptimalPositioningTask/DEIM"'
  55 fig.suptitle(name, fontsize=20, fontstyle='italic')
  56 #
  57 ax = fig.add_subplot(1, 3, 1)
  58 name = "Singular values of the set of %i G simulations"%sp.shape[1]
  59 print("  -", name)
  60 ax.set_title(name, fontstyle='italic', color='red')
  61 ax.set_xlabel("Index of singular values, numbered from 1")
  62 ax.set_ylabel("Magnitude of singular values")
  63 ax.set_xlim(1, len(sv))
  64 ax.set_yscale("log")
  65 ax.grid(True)
  66 ax.plot(range(1, 1 + len(sv)), sv)
  67 #
  68 ax = fig.add_subplot(1, 3, 2, projection = "3d")
  69 name = "A few simulations of G out of a total of %i"%sp.shape[1]
  70 print("  -", name)
  71 ax.set_title(name, fontstyle='italic', color='red')
  72 ax.set_xlabel("Position x1")
  73 ax.set_ylabel("Position x2")
  74 ax.set_zlabel("Magnitude of G")
  75 for i in range(sp.shape[1]):
  76     if i % 44 != 0: continue
  77     Gfield = sp[:,i].reshape((dimension, dimension))
  78     ax.plot_surface(x1, x2, Gfield, cmap='coolwarm')
  79 #
  80 ax = fig.add_subplot(1, 3, 3)
  81 name = "Set of %i first optimal points of measurement"%nbmeasures
  82 print("  -", name)
  83 ax.set_title(name, fontstyle='italic', color='red')
  84 ax.set_xlabel("Position x1")
  85 ax.set_ylabel("Position x2")
  86 ax.set_xlim(Omega[0][0] - 0.05, Omega[0][1] + 0.05)
  87 ax.set_ylim(Omega[1][0] - 0.05, Omega[1][1] + 0.05)
  88 ax.grid(True, which='both', linestyle=(0, (1, 5)), linewidth=0.5)
  89 ax.plot(posx1, posx2, markersize=6, marker="o", linestyle='')
  90 for i in range(len(posx1)):
  91     ax.text(posx1[i] + 0.005, posx2[i] + 0.01, str(i + 1), fontweight='bold')
  92 #
  93 plt.tight_layout()
  94 fig.savefig("simple_MeasurementsOptimalPositioningTask3.png")
  95 plt.close()