1 files changed, 84 insertions, 13 deletions
diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 90843b1..9c8f3ee 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -37,11 +37,42 @@ def find_shift(z, target):
     return zs, factor
 
 
-def laguerre_gamma_shift(z, x=None, w=None, n=8, target=11):
+def find_optimal_shift(z, n):
+    mhat = 1.34093 * n + 0.854093
+    steps = int(np.ceil(mhat - np.real(z))) - 1
+    return steps
+
+
+def get_shifting_factor(z, steps):
+    if steps > 0:
+        factor = 1 / pochhammer(z, steps)
+    elif steps < 0:
+        factor = pochhammer(z + steps, -steps)
+    return factor
+
+
+def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
     x, w = _prep_zeros_and_weights(x, w, n)
+    n = len(x)
 
     z += 0j
     z_shifted, correction_factor = find_shift(z, target)
+
+    res = np.sum(x ** (z_shifted - 1) * w)
+    res *= correction_factor
+    res = drop_imag(res)
+    return res
+
+
+def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
+    x, w = _prep_zeros_and_weights(x, w, n)
+    n = len(x)
+
+    z += 0j
+    opt_shift = find_optimal_shift(z, n)
+    correction_factor = get_shifting_factor(z, opt_shift)
+    z_shifted = z + opt_shift
+
     res = np.sum(x ** (z_shifted - 1) * w)
     res *= correction_factor
     res = drop_imag(res)
@@ -72,8 +103,10 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
         f = laguerre_gamma_simple
     elif func == "mirror":
         f = laguerre_gamma_mirror
+    elif func == "optimal_shifted":
+        f = laguerre_gamma_opt_shifted
     else:
-        f = laguerre_gamma_shift
+        f = laguerre_gamma_shifted
     return np.array([f(zi, x, w, n, **kwargs) for zi in z])
 
 
@@ -81,11 +114,10 @@ def calc_rel_error(x, y):
     return (y - x) / x
 
 
-ns = np.arange(2, 12, 2)
-
 # Simple / naive
 xmin = -5
 xmax = 30
+ns = np.arange(2, 12, 2)
 ylim = np.array([-11, 6])
 x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400)
 gamma = scipy.special.gamma(x)
@@ -104,7 +136,7 @@ ax.set_xlabel(r"$z$")
 ax.set_ylabel("Relativer Fehler")
 ax.legend(ncol=3, fontsize="small")
 ax.grid(1, "both")
-fig.savefig(f"{img_path}/rel_error_simple.pgf")
+# fig.savefig(f"{img_path}/rel_error_simple.pgf")
 
 
 # Mirrored
@@ -130,7 +162,7 @@ ax2.set_xlabel(r"$z$")
 ax2.set_ylabel("Relativer Fehler")
 ax2.legend(ncol=1, loc="upper left", fontsize="small")
 ax2.grid(1, "both")
-fig2.savefig(f"{img_path}/rel_error_mirror.pgf")
+# fig2.savefig(f"{img_path}/rel_error_mirror.pgf")
 
 
 # Move to target
@@ -163,12 +195,14 @@ v = ax3.imshow(bests, cmap="inferno", aspect="auto", interpolation="nearest")
 plt.colorbar(v, ax=ax3, label=r"$m$")
 ticks = np.arange(0, N + 1, N // 5)
 ax3.set_xlim(0, 1)
-ax3.set_xticks(ticks, [f"{v:.2f}" for v in ticks / N])
+ax3.set_xticks(ticks)
+ax3.set_xticklabels([f"{v:.2f}" for v in ticks / N])
 ax3.set_xticks(np.arange(0, N + 1, N // 20), minor=True)
-ax3.set_yticks(np.arange(len(ns)), ns)
+ax3.set_yticks(np.arange(len(ns)))
+ax3.set_yticklabels(ns)
 ax3.set_xlabel(r"$z$")
 ax3.set_ylabel(r"$n$")
-fig3.savefig(f"{img_path}/targets.pdf")
+# fig3.savefig(f"{img_path}/targets.pdf")
 
 targets = np.mean(bests, -1)
 intercept, bias = np.polyfit(ns, targets, 1)
@@ -178,9 +212,7 @@ fig4, axs4 = plt.subplots(
 xl = np.array([ns[0] - 0.5, ns[-1] + 0.5])
 axs4[0].plot(xl, intercept * xl + bias, label=r"$\hat{m}$")
 axs4[0].plot(ns, targets, "x", label=r"$\bar{m}$")
-axs4[1].plot(
-    ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \bar{m}$"
-)
+axs4[1].plot(ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \bar{m}$")
 axs4[0].set_xlim(*xl)
 # axs4[0].set_title("Schätzung von Mittelwert")
 # axs4[1].set_title("Fehler")
@@ -188,10 +220,49 @@ axs4[-1].set_xlabel(r"$z$")
 for ax in axs4:
     ax.grid(1)
     ax.legend()
-fig4.savefig(f"{img_path}/schaetzung.pgf")
+# fig4.savefig(f"{img_path}/schaetzung.pgf")
 
 print(f"Intercept={intercept:.6g}, Bias={bias:.6g}")
 predicts = np.ceil(intercept * ns[:, None] + bias - x)
 print(f"Error: {int(np.sum(np.abs(bests-predicts)))}")
 
+# Comparison relative error between methods
+N = 200
+step = 1 / (N - 1)
+x = np.linspace(step, 1 - step, N + 1)
+gamma = scipy.special.gamma(x)[:, None]
+n = 8
+targets = np.arange(10, 14)
+gamma = scipy.special.gamma(x)
+fig5, ax5 = plt.subplots(num=1, clear=True, constrained_layout=True)
+for target in targets:
+    gamma_lag = eval_laguerre_gamma(x, target=target, n=n, func="shifted")
+    rel_error = np.abs(calc_rel_error(gamma, gamma_lag))
+    ax5.semilogy(x, rel_error, label=f"$m={target}$")
+gamma_lgo = eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+rel_error = np.abs(calc_rel_error(gamma, gamma_lgo))
+ax5.semilogy(x, rel_error, label="$m^*$")
+ax5.set_xlim(x[0], x[-1])
+ax5.set_ylim(5e-9, 5e-8)
+ax5.set_xlabel(r"$z$")
+ax5.grid(1, "both")
+ax5.legend()
+fig5.savefig(f"{img_path}/rel_error_shifted.pgf")
+
+N = 200
+x = np.linspace(-5+ EPSILON, 5-EPSILON, N)
+gamma = scipy.special.gamma(x)[:, None]
+n = 8
+gamma = scipy.special.gamma(x)
+fig6, ax6 = plt.subplots(num=1, clear=True, constrained_layout=True)
+gamma_lgo = eval_laguerre_gamma(x, n=n, func="optimal_shifted")
+rel_error = np.abs(calc_rel_error(gamma, gamma_lgo))
+ax6.semilogy(x, rel_error, label="$m^*$")
+ax6.set_xlim(x[0], x[-1])
+ax6.set_ylim(5e-9, 5e-8)
+ax6.set_xlabel(r"$z$")
+ax6.grid(1, "both")
+ax6.legend()
+fig6.savefig(f"{img_path}/rel_error_range.pgf")
+
 # plt.show()