From 127b89c35801ee7f7e4ab3ace2649b97eb6b7091 Mon Sep 17 00:00:00 2001 From: Elijah Andrews Date: Wed, 25 Mar 2020 14:53:59 +0000 Subject: [PATCH] Various curve fit testing improvements. --- numerical/testing/error_simulation.py | 85 ++++++++++++------- numerical/testing/simulated_peak_errors.py | 93 +++++++++++++++----- numerical/testing/std_vs_N.py | 99 ++++++++++++++++++++++ 3 files changed, 223 insertions(+), 54 deletions(-) create mode 100644 numerical/testing/std_vs_N.py diff --git a/numerical/testing/error_simulation.py b/numerical/testing/error_simulation.py index 19c3f28..e6ecbf1 100644 --- a/numerical/testing/error_simulation.py +++ b/numerical/testing/error_simulation.py @@ -8,7 +8,7 @@ from numerical.models.slot.slot_opt import find_slot_peak from experimental.plotting.analyse_slot import SweepData m_0 = 1 -n = 5000 +n = 10000 w_thresh = 15 density_ratio = 0.15 H = 2 @@ -20,64 +20,87 @@ print("Requested n = {0}, using n = {1}.".format(n, len(centroids))) R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32) R_inv = scipy.linalg.inv(R_matrix) -xs = np.linspace(-5, 3, 16) +xs = np.linspace(-12, 8, 32) Xs = xs * 0.5 * W -Y = W +Y = W * 1.5 points = np.empty((len(Xs), 3)) points[:, 0] = Xs points[:, 1] = Y points[:, 2] = 0 vels = bem.get_jet_dirs(points, centroids, normals, areas, m_0, R_inv) -theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi +thetas = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi _, x_star, theta_j_star, _ = find_slot_peak(W, Y, H, n, centroids=centroids, normals=normals, areas=areas, R_inv=R_inv) fake_sweep = SweepData('test', Y, W, H) std = 0.015085955056793596 -N = 5 +N = 10 all_rand_xs = [] all_rand_theta_js = [] -for p_bar, theta_j in zip(xs, theta_js): +for p_bar, theta_j in zip(xs, thetas): rand_theta_js = np.random.normal(theta_j, std, N) for rand_theta_j in rand_theta_js: all_rand_xs.append(p_bar) all_rand_theta_js.append(rand_theta_j) - fake_sweep.add_point(p_bar, p_bar, rand_theta_js, Y) - -(fit_max_peak_x, fit_max_peak_theta_j, max_poly_coeffs, max_theta_j_std), \ -(fit_min_peak_x, fit_min_peak_theta_j, min_poly_coeffs, min_theta_j_std) = fake_sweep.get_curve_fits(range_fact=1.5) - + fake_sweep.add_point(p_bar, p_bar, rand_theta_j, Y) + +peak_fit_range = 1.5 +(fit_max_peak_x, fit_max_peak_theta_j, max_poly_coeffs), \ +(fit_min_peak_x, fit_min_peak_theta_j, min_poly_coeffs), \ +combined_tj_std, _ \ + = fake_sweep.get_curve_fits(range_fact=peak_fit_range) +# fit_x_star = (fit_max_peak_x - fit_min_peak_x) / 2 -fit_theta_j_star = (fit_max_peak_theta_j - fit_min_peak_theta_j) / 2 +fit_theta_star = (fit_max_peak_theta_j - fit_min_peak_theta_j) / 2 +# +# fit_x_star, fit_theta_star, x_offset, theta_offset = num_prediction_fit(xs, thetas, std) +# num_interp = get_num_prediction_function() +# +# print(x_offset, theta_offset) -print(max_theta_j_std, min_theta_j_std) +# print(max_theta_j_std, min_theta_j_std) -peak_fit_range = 2.5 -max_peak_x = xs[np.argmax(theta_js)] +max_peak_x = xs[np.argmax(thetas)] max_fit_xs = np.linspace(0, peak_fit_range * max_peak_x, 100) max_fit_ys = np.polyval(max_poly_coeffs, max_fit_xs) -min_peak_x = xs[np.argmin(theta_js)] -min_fit_xs = np.linspace(0, peak_fit_range * min_peak_x, 100) -min_fit_ys = np.polyval(min_poly_coeffs, min_fit_xs) +# min_peak_x = xs[np.argmin(thetas)] +# min_fit_xs = np.linspace(0, peak_fit_range * min_peak_x, 100) +# min_fit_ys = np.polyval(min_poly_coeffs, min_fit_xs) + +# fit_xs = np.linspace(min(xs), max(xs), 500) +# fit_ys = num_interp(fit_xs, fit_x_star, fit_theta_star, x_offset, theta_offset) plt.scatter(all_rand_xs, all_rand_theta_js, marker='.', label="Randomized") -plt.scatter(xs, theta_js, label="Numerical", marker='x') +plt.scatter(xs, thetas, label="Numerical", marker='x') plt.scatter([x_star, -x_star], [theta_j_star, -theta_j_star], label="Numerical Peaks", color='g', marker='+') -plt.plot(max_fit_xs, max_fit_ys, label="Curve Fits", color='purple') -plt.plot(min_fit_xs, min_fit_ys, color='purple') - -plt.axvline(x_star, color='g', linestyle='--') -plt.axvline(fit_x_star, color='r', linestyle='--') -plt.axhline(theta_j_star, color='g', linestyle='--') -plt.axhline(fit_theta_j_star, color='r', linestyle='--') - -plt.scatter([fit_max_peak_x, fit_min_peak_x], [fit_max_peak_theta_j, fit_min_peak_theta_j], label="Fitted Peaks", color='k', - marker='+', s=100) - +plt.plot(max_fit_xs, max_fit_ys, label="Curve Fit", color='purple') +confidence_interval = 0.99 +# plt.plot(max_fit_xs, np.add(max_fit_ys, stats.norm.interval(confidence_interval, loc=0, scale=max_theta_j_std)[0]), +# color='gray', linestyle='--') +# plt.plot(max_fit_xs, np.add(max_fit_ys, stats.norm.interval(confidence_interval, loc=0, scale=max_theta_j_std)[1]), +# color='gray', linestyle='--') +# plt.plot(np.add(max_fit_xs, stats.norm.interval(confidence_interval, loc=0, scale=max_x_star_std)[0]), max_fit_ys, +# color='gray', linestyle='--') +# plt.plot(np.add(max_fit_xs, stats.norm.interval(confidence_interval, loc=0, scale=max_x_star_std)[1]), max_fit_ys, +# color='gray', linestyle='--') +# plt.plot(min_fit_xs, min_fit_ys, color='purple') + +# plt.axvline(x_star, color='g', linestyle='--') +# plt.axvline(fit_x_star, color='r', linestyle='--') +# plt.axhline(theta_j_star, color='g', linestyle='--') +# plt.axhline(fit_theta_j_star, color='r', linestyle='--') + +# plt.errorbar([fit_max_peak_x, fit_min_peak_x], [fit_max_peak_theta_j, fit_min_peak_theta_j], label="Fitted Peaks", +# yerr=[stats.norm.interval(confidence_interval, loc=0, scale=max_theta_j_std)[0], +# stats.norm.interval(confidence_interval, loc=0, scale=min_theta_j_std)[0]], +# xerr=[stats.norm.interval(confidence_interval, loc=0, scale=max_x_star_std)[0], +# stats.norm.interval(confidence_interval, loc=0, scale=min_x_star_std)[0]], +# color='k', capsize=3, fmt='.') +plt.axvline(0, color='k', linestyle='--') plt.xlabel("$x$") -plt.ylabel("$\\theta_j$") +plt.ylabel("$\\theta_j$ (rad)") plt.legend() plt.show() diff --git a/numerical/testing/simulated_peak_errors.py b/numerical/testing/simulated_peak_errors.py index c953761..b58f1a0 100644 --- a/numerical/testing/simulated_peak_errors.py +++ b/numerical/testing/simulated_peak_errors.py @@ -8,7 +8,7 @@ from numerical.models.slot.slot_opt import find_slot_peak from experimental.plotting.analyse_slot import SweepData m_0 = 1 -n = 5000 +n = 10000 w_thresh = 15 density_ratio = 0.15 H = 2 @@ -20,9 +20,9 @@ print("Requested n = {0}, using n = {1}.".format(n, len(centroids))) R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32) R_inv = scipy.linalg.inv(R_matrix) -xs = np.linspace(-5, 3, 16) +xs = np.linspace(-5, 3, 32) Xs = xs * 0.5 * W -Y = W +Y = W * 0.25 points = np.empty((len(Xs), 3)) points[:, 0] = Xs points[:, 1] = Y @@ -34,52 +34,99 @@ _, x_star, theta_j_star, _ = find_slot_peak(W, Y, H, n, centroids=centroids, nor R_inv=R_inv) std = 0.015085955056793596 -N = 5 - -num_sweeps = 10000 +N = 10 +num_sweeps = 100000 fit_x_stars = [] -fit_theta_j_stars = [] +fit_theta_stars = [] +fit_theta_stds = [] +max_mean_xs = [] for i in range(num_sweeps): + if i % 100 == 0: + print(f"{100 * i / num_sweeps}%") fake_sweep = SweepData('test', Y, W, H) all_rand_xs = [] all_rand_theta_js = [] + max_mean_x = 0 + max_mean_theta_j = 0 for x, theta_j in zip(xs, theta_js): rand_theta_js = np.random.normal(theta_j, std, N) + if np.mean(rand_theta_js) > max_mean_theta_j: + max_mean_theta_j = np.mean(rand_theta_js) + max_mean_x = x for rand_theta_j in rand_theta_js: all_rand_xs.append(x) all_rand_theta_js.append(rand_theta_j) - fake_sweep.add_point(x, x, rand_theta_js, Y) - + fake_sweep.add_point(x, x, rand_theta_j, Y) try: - (max_peak_x, max_peak_theta_j, max_poly_coeffs, max_theta_j_std), \ - (min_peak_x, min_peak_theta_j, min_poly_coeffs, min_theta_j_std) = fake_sweep.get_curve_fits(range_fact=0.5) + # fit_x_star, fit_theta_star, _, _ = num_prediction_fit(all_rand_xs, all_rand_theta_js, std) + (max_peak_x, max_peak_theta_j, max_poly_coeffs), \ + (min_peak_x, min_peak_theta_j, min_poly_coeffs), \ + combined_theta_j_std, _ = fake_sweep.get_curve_fits(range_fact=1.5) + # fake_sweep.get_curve_fits(range_fact=1.75 * 1.117647058823529 / max_mean_x) except ValueError: continue - fit_p_bar_star = (max_peak_x - min_peak_x) / 2 - fit_theta_j_star = (max_peak_theta_j - min_peak_theta_j) / 2 + if np.abs((max_peak_x - min_peak_x) / 2) > 2 * Y: + continue + + max_mean_xs.append(max_mean_x) + fit_x_star = (max_peak_x - min_peak_x) / 2 + fit_theta_star = (max_peak_theta_j - min_peak_theta_j) / 2 - fit_x_stars.append(fit_p_bar_star) - fit_theta_j_stars.append(fit_theta_j_star) + fit_x_stars.append(fit_x_star) + fit_theta_stars.append(fit_theta_star) -plt.scatter(fit_x_stars, fit_theta_j_stars, marker='x', color='k', label="Randomized Fitted Peaks") + fit_theta_stds.append(combined_theta_j_std) + # fit_theta_stds.append(max_theta_j_std) + # fit_x_star_stds.append(max_x_star_std) +plt.scatter(fit_x_stars, fit_theta_stars, marker='.', c=max_mean_xs, label="Randomized Fitted Peaks") plt.scatter([x_star], [theta_j_star], label="Numerical Peak", color='r', marker='o') +print(f"Mean error: {np.mean(fit_x_stars) - x_star}") + plt.xlabel("$x$") -plt.ylabel("$\\theta_j$") +plt.ylabel("$\\theta_j$ (rad)") plt.legend() +confidence_interval = 0.99 + fig, axes = plt.subplots(nrows=2, ncols=1) axes[0].hist(fit_x_stars, color='k', bins=int(np.ceil(num_sweeps / 10))) axes[0].axvline(x_star, color='r') +# x_star_int = stats.norm.interval(confidence_interval, loc=np.mean(fit_x_stars), scale=np.mean(fit_x_star_stds)) +# axes[0].axvline(x_star_int[0], color='g', label="Estimated") +# axes[0].axvline(x_star_int[1], color='g') +axes[0].axvline(np.mean(fit_x_stars) + 0.075, color='g', label="Estimated") +axes[0].axvline(np.mean(fit_x_stars) - 0.075, color='g') +# x_star_int = stats.norm.interval(confidence_interval, loc=np.mean(fit_x_stars), scale=np.std(fit_x_stars)) +# axes[0].axvline(x_star_int[0], color='b', label="Measured") +# axes[0].axvline(x_star_int[1], color='b') +axes[0].axvline(np.mean(fit_x_stars) + np.std(fit_x_stars), color='b', label="Measured") +axes[0].axvline(np.mean(fit_x_stars) - np.std(fit_x_stars), color='b') axes[0].set_xlabel("$x^\\star$") +axes[0].legend() -axes[1].hist(fit_theta_j_stars, color='k', bins=int(np.ceil(num_sweeps / 10))) +axes[1].hist(fit_theta_stars, color='k', bins=int(np.ceil(num_sweeps / 10))) axes[1].axvline(theta_j_star, color='r') -axes[1].set_xlabel("$\\theta_j^\\star$") - -print(np.std(fit_x_stars)) -print(np.std(fit_theta_j_stars)) - +# theta_j_star_int = stats.norm.interval(confidence_interval, loc=np.mean(fit_theta_j_stars), +# scale=np.mean(fit_theta_j_stds)) +# axes[1].axvline(theta_j_star_int[0], color='g') +# axes[1].axvline(theta_j_star_int[1], color='g') +axes[1].axvline(np.mean(fit_theta_stars) + np.mean(fit_theta_stds), color='g') +axes[1].axvline(np.mean(fit_theta_stars) - np.mean(fit_theta_stds), color='g') +# theta_j_star_int = stats.norm.interval(confidence_interval, loc=np.mean(fit_theta_j_stars), +# scale=np.std(fit_theta_j_stars)) +# axes[1].axvline(theta_j_star_int[0], color='b') +# axes[1].axvline(theta_j_star_int[1], color='b') +axes[1].axvline(np.mean(fit_theta_stars) + np.std(fit_theta_stars), color='b') +axes[1].axvline(np.mean(fit_theta_stars - np.std(fit_theta_stars)), color='b') +axes[1].set_xlabel("$\\theta_j^\\star$ (rad)") + +# print(f"x_star: measured = {np.std(fit_x_stars)}, estimated = {np.mean(fit_x_star_stds)}") +# print(np.std(fit_x_star_stds)) +# print(f"theta_j_star: measured = {np.std(fit_theta_stars)}, estimated = {np.mean(fit_theta_j_stds)}") +# print(np.std(fit_theta_j_stds)) + +plt.tight_layout() plt.show() diff --git a/numerical/testing/std_vs_N.py b/numerical/testing/std_vs_N.py new file mode 100644 index 0000000..3d1d83b --- /dev/null +++ b/numerical/testing/std_vs_N.py @@ -0,0 +1,99 @@ +import scipy +import numpy as np +import matplotlib.pyplot as plt + +import numerical.util.gen_utils as gen +import numerical.bem as bem +from numerical.models.slot.slot_opt import find_slot_peak +from experimental.plotting.analyse_slot import SweepData + +m_0 = 1 +n = 5000 +w_thresh = 15 +density_ratio = 0.15 +H = 2 +W = 2 + +centroids, normals, areas = gen.gen_varied_slot(n=n, H=H, W=W, length=100, depth=50, w_thresh=w_thresh, + density_ratio=density_ratio) +print("Requested n = {0}, using n = {1}.".format(n, len(centroids))) +R_matrix = bem.get_R_matrix(centroids, normals, areas, dtype=np.float32) +R_inv = scipy.linalg.inv(R_matrix) + +xs = np.linspace(-5, 3, 16) +Xs = xs * 0.5 * W +Y = W +points = np.empty((len(Xs), 3)) +points[:, 0] = Xs +points[:, 1] = Y +points[:, 2] = 0 +vels = bem.get_jet_dirs(points, centroids, normals, areas, m_0, R_inv) +theta_js = np.arctan2(vels[:, 1], vels[:, 0]) + 0.5 * np.pi + +_, x_star, theta_j_star, _ = find_slot_peak(W, Y, H, n, centroids=centroids, normals=normals, areas=areas, + R_inv=R_inv) + +std = 0.015085955056793596 +num_sweeps = 50000 +Ns = range(10, 50, 5) +theta_j_star_stds = [] +x_star_stds = [] +mean_Ns = [] +for N in Ns: + print(N) + fit_x_stars = [] + fit_theta_j_stars = [] + fit_theta_j_stds = [] + fit_Ns = [] + for i in range(num_sweeps): + fake_sweep = SweepData('test', Y, W, H) + all_rand_xs = [] + all_rand_theta_js = [] + for x, theta_j in zip(xs, theta_js): + rand_theta_js = np.random.normal(theta_j, std, N) + for rand_theta_j in rand_theta_js: + all_rand_xs.append(x) + all_rand_theta_js.append(rand_theta_j) + fake_sweep.add_point(x, x, rand_theta_j, Y) + + try: + (max_peak_x, max_peak_theta_j, max_poly_coeffs), \ + (min_peak_x, min_peak_theta_j, min_poly_coeffs), \ + combined_theta_j_star_std, N_points = fake_sweep.get_curve_fits(range_fact=1.5) + except ValueError: + continue + + fit_Ns.append(N_points) + fit_x_star = (max_peak_x - min_peak_x) / 2 + fit_theta_j_star = (max_peak_theta_j - min_peak_theta_j) / 2 + # fit_p_bar_star = max_peak_x + # fit_theta_j_star = max_peak_theta_j + + fit_x_stars.append(fit_x_star) + fit_theta_j_stars.append(fit_theta_j_star) + + fit_theta_j_stds.append(combined_theta_j_star_std) + # fit_theta_j_stds.append((max_theta_j_std + min_theta_j_std) / 2) + # fit_x_star_stds.append((max_x_star_std + min_x_star_std) / 2) + # fit_theta_j_stds.append(max_theta_j_std) + # fit_x_star_stds.append(max_x_star_std) + mean_Ns.append(np.mean(fit_Ns)) + theta_j_star_stds.append(np.std(fit_theta_j_stars)) + x_star_stds.append(np.std(fit_x_stars)) + +print(np.polyfit(np.log(mean_Ns), np.log(np.divide(theta_j_star_stds, std)), 1)) +print(np.polyfit(np.log(mean_Ns), np.log(np.divide(x_star_stds, std)), 1)) + +plt.plot(mean_Ns, theta_j_star_stds, 'C0') +# plt.plot(np.multiply(Ns, 4), np.multiply(theta_j_star_stds, np.sqrt(np.multiply(Ns, 4))) / std) +plt.ylabel("$\\theta_j^\\star \\sigma$", color='C0') +plt.xlabel("$N$ - Samples per position") +plt.xscale('log') +plt.yscale('log') +plt.twinx() +plt.plot(mean_Ns, x_star_stds, 'C1') +# plt.plot(np.multiply(Ns, 4), np.multiply(x_star_stds, np.sqrt(np.multiply(Ns, 4))) / (25 * std)) +plt.ylabel("$x^\\star \\sigma$", color='C1') +plt.yscale('log') +plt.tight_layout() +plt.show() -- GitLab