convergence to projection for PDE, closes #36

Author: mathematicalmichael

.gitignore

@@ -62,3 +62,4 @@ pde_*D/
 *-data/
 *.pkl
 
+data

setup.cfg

@@ -64,8 +64,8 @@ pub =
 console_scripts =
     mud_examples = mud_examples.runner:run
     generate_poisson_data = mud_examples.poisson:run
-    mud_run_inv = mud_examples.runner:run_inv
-    mud_run_lin = mud_examples.runner:run_lin
+    mud_run_inv = mud_examples.runner:run_monomial
+    mud_run_lin = mud_examples.runner:run_linear
     mud_run_pde = mud_examples.runner:run_pde
     mud_run_ode = mud_examples.runner:run_ode
     mud_run_all = mud_examples.runner:run_all

src/mud_examples/experiments.py

@@ -21,8 +21,11 @@
 _mpl_logger.setLevel(logging.WARNING)
 
 
-def experiment_measurements(fun, num_measurements,
-                            sd, num_trials, seed=21):
+def experiment_measurements(fun,
+                            num_measurements,
+                            sd,
+                            num_trials,
+                            seed=21):
     """
     Fixed sensors, varying how much data is incorporated into the solution.
     """
@@ -44,79 +47,84 @@ def experiment_measurements(fun, num_measurements,
     return experiments, solutions
 
 
-def experiment_equipment(fun, num_measure,
-                         sd_vals, num_trials,
-                         reference_value):
+def experiment_equipment(fun,
+                         num_measure,
+                         sd_vals,
+                         num_trials,
+                         seed=21):
     """
     Fixed number of sensors, varying the quality of equipment.
     """
-    sd_err = []
-    sd_var = []
+    experiments = {}
+    solutions = {}
     for sd in sd_vals:
         _logger.debug(f'Equipment Experiment. Std Dev: {sd}')
-        temp_err = []
+        discretizations = []
+        estimates = []
         for t in range(num_trials):
+            np.random.seed(seed + t)
             _d = fun(sd=sd, num_obs=num_measure)
             estimate = _d.estimate()
-            temp_err.append(np.linalg.norm(estimate - reference_value))
-        sd_err.append(np.mean(temp_err))
-        sd_var.append(np.var(temp_err))
+            discretizations.append(_d)
+            estimates.append(estimate)
+        experiments[sd] = discretizations
+        solutions[sd] = estimates
 
-    return sd_err, sd_var
+    return experiments, solutions
 
 
 def plot_experiment_equipment(tolerances, res, prefix, fsize=32, linewidth=5,
                               title="Variance of MUD Error", save=True):
     print("Plotting experiments involving equipment differences...")
     plt.figure(figsize=(10, 10))
     for _res in res:
-        _prefix, _in, _rm, _re = _res
+        _example, _in, _rm, _re, _fname = _res
         regression_err_mean, slope_err_mean, \
             regression_err_vars, slope_err_vars, \
             sd_means, sd_vars, num_sensors = _re
         plt.plot(tolerances, regression_err_mean,
-                 label=f"{_prefix:10s} slope: {slope_err_mean:1.4f}",
+                 label=f"{_example.upper()} slope: {slope_err_mean:1.4f}",
                  lw=linewidth)
         plt.scatter(tolerances, sd_means, marker='x', lw=20)
 
     plt.yscale('log')
     plt.xscale('log')
     plt.Axes.set_aspect(plt.gca(), 1)
-    plt.ylim(2E-3, 2E-2)
+    # plt.ylim(2E-3, 2E-2)
     # plt.ylabel("Absolute Error", fontsize=fsize)
     plt.xlabel('Tolerance', fontsize=fsize)
     plt.legend()
     plt.title(f"Mean of MUD Error for N={num_sensors}", fontsize=1.25 * fsize)
     if save:
         fdir = ''.join(prefix.split('/')[::-1])
-        check_dir(fdir)
+        check_dir(f'figures/{_fname}/{fdir}')
         _logger.info("Saving equipment experiments: mean convergence.")
-        plt.savefig(f'{prefix}_convergence_mud_std_mean.png',
+        plt.savefig(f'figures/{_fname}/{prefix}_convergence_mud_std_mean.png',
                     bbox_inches='tight')
     else:
         plt.show()
 
     plt.figure(figsize=(10, 10))
     for _res in res:
-        _prefix, _in, _rm, _re = _res
+        _example, _in, _rm, _re, _fname = _res
         regression_err_mean, slope_err_mean, \
             regression_err_vars, slope_err_vars, \
             sd_means, sd_vars, num_sensors = _re
         plt.plot(tolerances, regression_err_vars,
-                 label=f"{_prefix:10s} slope: {slope_err_vars:1.4f}",
+                 label=f"{_example.upper()} slope: {slope_err_vars:1.4f}",
                  lw=linewidth)
         plt.scatter(tolerances, sd_vars, marker='x', lw=20)
     plt.xscale('log')
     plt.yscale('log')
-    plt.ylim(2E-5, 2E-4)
+    # plt.ylim(2E-5, 2E-4)
     plt.Axes.set_aspect(plt.gca(), 1)
     # plt.ylabel("Absolute Error", fontsize=fsize)
     plt.xlabel('Tolerance', fontsize=fsize)
     plt.legend()
     plt.title(title, fontsize=1.25 * fsize)
     if save:
         _logger.info("Saving equipment experiments: variance convergence.")
-        plt.savefig(f'{prefix}_convergence_mud_std_var.png',
+        plt.savefig(f'figures/{_fname}/{prefix}_convergence_mud_std_var.png',
                     bbox_inches='tight')
     else:
         plt.show()
@@ -125,25 +133,25 @@ def plot_experiment_equipment(tolerances, res, prefix, fsize=32, linewidth=5,
 def plot_experiment_measurements(res, prefix,
                                  fsize=32, linewidth=5,
                                  xlabel='Number of Measurements',
-                                 save=True, legend=False):
+                                 save=True, legend=True):
     print("Plotting experiments involving increasing # of measurements.")
     plt.figure(figsize=(10, 10))
     for _res in res:
-        _prefix, _in, _rm, _re = _res
+        _example, _in, _rm, _re, _fname = _res
         solutions = _in[-1]
         measurements = list(solutions.keys())
         regression_mean, slope_mean, \
             regression_vars, slope_vars, \
             means, variances = _rm
         plt.plot(measurements[:len(regression_mean)], regression_mean,
-                 label=f"{_prefix:4s} slope: {slope_mean:1.4f}",
+                 label=f"{_example.upper()} slope: {slope_mean:1.4f}",
                  lw=linewidth)
         plt.scatter(measurements[:len(means)], means, marker='x', lw=20)
     plt.xscale('log')
     plt.yscale('log')
     plt.Axes.set_aspect(plt.gca(), 1)
-    plt.ylim(0.9 * min(means), 1.3 * max(means))
-    plt.ylim(2E-3, 2E-1)
+    # plt.ylim(0.9 * min(means), 1.3 * max(means))
+    # plt.ylim(2E-3, 2E-1)
     plt.xlabel(xlabel, fontsize=fsize)
     if legend:
         plt.legend(fontsize=fsize * 0.8)
@@ -152,20 +160,20 @@ def plot_experiment_measurements(res, prefix,
     plt.title(title, fontsize=1.25 * fsize)
     if save:
         fdir = '/'.join(prefix.split('/')[:-1])
-        check_dir(fdir)
+        check_dir(f'figures/{_fname}/{fdir}')
         _logger.info("Saving measurement experiments: mean convergence.")
-        plt.savefig(f'{prefix}_convergence_obs_mean.png', bbox_inches='tight')
+        plt.savefig(f'figures/{_fname}/{prefix}_convergence_obs_mean.png', bbox_inches='tight')
     else:
         plt.show()
 
     plt.figure(figsize=(10, 10))
     for _res in res:
-        _prefix, _in, _rm, _re = _res
+        _example, _in, _rm, _re, _fname = _res
         regression_mean, slope_mean, \
             regression_vars, slope_vars, \
             means, variances = _rm
         plt.plot(measurements[:len(regression_vars)], regression_vars,
-                 label=f"{_prefix:4s} slope: {slope_vars:1.4f}",
+                 label=f"{_example.upper()} slope: {slope_vars:1.4f}",
                  lw=linewidth)
         plt.scatter(measurements[:len(variances)], variances,
                     marker='x', lw=20)
@@ -174,7 +182,7 @@ def plot_experiment_measurements(res, prefix,
     plt.Axes.set_aspect(plt.gca(), 1)
 #     if not len(np.unique(variances)) == 1:
 #         plt.ylim(0.9 * min(variances), 1.3 * max(variances))
-    plt.ylim(5E-6, 5E-4)
+    # plt.ylim(5E-6, 5E-4)
     plt.xlabel(xlabel, fontsize=fsize)
     if legend:
         plt.legend(fontsize=fsize * 0.8)
@@ -183,6 +191,6 @@ def plot_experiment_measurements(res, prefix,
               fontsize=1.25 * fsize)
     if save:
         _logger.info("Saving measurement experiments: variance convergence.")
-        plt.savefig(f'{prefix}_convergence_obs_var.png', bbox_inches='tight')
+        plt.savefig(f'figures/{_fname}/{prefix}_convergence_obs_var.png', bbox_inches='tight')
     else:
         plt.show()

src/mud_examples/ode.py

@@ -91,11 +91,13 @@ def wrapper(num_obs, sd):
             raise ValueError("Unknown example type")
 
         _logger.info("Increasing Measurements Quantity Study")
-        experiments, solutions = experiment_measurements(num_measurements=measurements,
-                                                 sd=sigma,
-                                                 num_trials=num_trials,
-                                                 seed=seed,
-                                                 fun=wrapper)
+        experiments, solutions = experiment_measurements(
+            num_measurements=measurements,
+            sd=sigma,
+            num_trials=num_trials,
+            seed=seed,
+            fun=wrapper,
+            )
 
         means, variances = extract_statistics(solutions, lam_true)
         regression_mean, slope_mean = fit_log_linear_regression(measurements, means)
@@ -105,14 +107,16 @@ def wrapper(num_obs, sd):
 
         num_sensors = num_measure
         if len(tolerances) > 1:
-            
             _logger.info("Increasing Measurement Precision Study")
-            sd_means, sd_vars = experiment_equipment(num_trials=num_trials,
-                                                  num_measure=num_sensors,
-                                                  sd_vals=sd_vals,
-                                                  reference_value=lam_true,
-                                                  fun=wrapper)
-
+            experiments, solutions = experiment_equipment(
+                sd_vals=sd_vals,
+                num_measure=num_sensors,
+                num_trials=num_trials,
+                seed=seed,
+                fun=wrapper,
+                )
+
+            sd_means, sd_vars = extract_statistics(solutions, lam_true)
             regression_err_mean, slope_err_mean = fit_log_linear_regression(tolerances, sd_means)
             regression_err_vars, slope_err_vars = fit_log_linear_regression(tolerances, sd_vars)
             _re = (regression_err_mean, slope_err_mean,
@@ -125,7 +129,7 @@ def wrapper(num_obs, sd):
         _in = (lam, qoi, sensors, qoi_true, experiments, solutions)
         _rm = (regression_mean, slope_mean, regression_vars, slope_vars, means, variances)
         example_name = example.upper()
-        res.append((example_name, _in, _rm, _re))
+        res.append((example_name, _in, _rm, _re, ''))
 
         # TODO check for existence of save directory, grab subset of measurements properly.
         plot_decay_solution(solutions, generate_decay_model, fsize=fsize,