Various curve fit testing improvements.

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()

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()

numerical/testing/std_vs_N.py

+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()