Skip to content

Commit

Permalink
convergence to projection for PDE, closes #36
Browse files Browse the repository at this point in the history
  • Loading branch information
@mathematicalmichael
mathematicalmichael authored Mar 23, 2021
2 parents dbb70b9 + 9690439 commit 41d41c1
Show file tree
Hide file tree
Showing 9 changed files with 273 additions and 147 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,3 +62,4 @@ pde_*D/
*-data/
*.pkl

data
4 changes: 2 additions & 2 deletions setup.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
72 changes: 40 additions & 32 deletions src/mud_examples/experiments.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
"""
Expand All @@ -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()
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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()
30 changes: 17 additions & 13 deletions src/mud_examples/ode.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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,
Expand All @@ -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,
Expand Down
Loading

0 comments on commit 41d41c1

Please sign in to comment.