From 6149839224755c21225d2decddeae12207c2cbab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Tue, 31 May 2022 16:31:25 +0200
Subject: Add rule of thumb, analyse integrand, correct mistake in integration
 SLP<->LP

---
 buch/papers/laguerre/scripts/gamma_approx.py | 197 +++++++++++++++++++++++++++
 1 file changed, 197 insertions(+)
 create mode 100644 buch/papers/laguerre/scripts/gamma_approx.py

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
new file mode 100644
index 0000000..90843b1
--- /dev/null
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -0,0 +1,197 @@
+from pathlib import Path
+
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+import scipy.special
+
+EPSILON = 1e-7
+root = str(Path(__file__).parent)
+img_path = f"{root}/../images"
+
+
+def _prep_zeros_and_weights(x, w, n):
+    if x is None or w is None:
+        return np.polynomial.laguerre.laggauss(n)
+    return x, w
+
+
+def drop_imag(z):
+    if abs(z.imag) <= EPSILON:
+        z = z.real
+    return z
+
+
+def pochhammer(z, n):
+    return np.prod(z + np.arange(n))
+
+
+def find_shift(z, target):
+    factor = 1.0
+    steps = int(np.floor(target - np.real(z)))
+    zs = z + steps
+    if steps > 0:
+        factor = 1 / pochhammer(z, steps)
+    elif steps < 0:
+        factor = pochhammer(zs, -steps)
+    return zs, factor
+
+
+def laguerre_gamma_shift(z, x=None, w=None, n=8, target=11):
+    x, w = _prep_zeros_and_weights(x, w, n)
+
+    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_simple(z, x=None, w=None, n=8):
+    x, w = _prep_zeros_and_weights(x, w, n)
+    z += 0j
+    res = np.sum(x ** (z - 1) * w)
+    res = drop_imag(res)
+    return res
+
+
+def laguerre_gamma_mirror(z, x=None, w=None, n=8):
+    x, w = _prep_zeros_and_weights(x, w, n)
+    z += 0j
+    if z.real < 1e-3:
+        return np.pi / (
+            np.sin(np.pi * z) * laguerre_gamma_simple(1 - z, x, w)
+        )  # Reflection formula
+    return laguerre_gamma_simple(z, x, w)
+
+
+def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
+    x, w = _prep_zeros_and_weights(x, w, n)
+    if func == "simple":
+        f = laguerre_gamma_simple
+    elif func == "mirror":
+        f = laguerre_gamma_mirror
+    else:
+        f = laguerre_gamma_shift
+    return np.array([f(zi, x, w, n, **kwargs) for zi in z])
+
+
+def calc_rel_error(x, y):
+    return (y - x) / x
+
+
+ns = np.arange(2, 12, 2)
+
+# Simple / naive
+xmin = -5
+xmax = 30
+ylim = np.array([-11, 6])
+x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400)
+gamma = scipy.special.gamma(x)
+fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
+for n in ns:
+    gamma_lag = eval_laguerre_gamma(x, n=n)
+    rel_err = calc_rel_error(gamma, gamma_lag)
+    ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+ax.set_xlim(x[0], x[-1])
+ax.set_ylim(*(10.0 ** ylim))
+ax.set_xticks(np.arange(xmin, xmax + EPSILON, 5))
+ax.set_xticks(np.arange(xmin, xmax), minor=True)
+ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+ax.set_yticks(10.0 ** np.arange(*ylim, 2))
+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")
+
+
+# Mirrored
+xmin = -15
+xmax = 15
+ylim = np.array([-11, 1])
+x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400)
+gamma = scipy.special.gamma(x)
+fig2, ax2 = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(5, 2.5))
+for n in ns:
+    gamma_lag = eval_laguerre_gamma(x, n=n, func="mirror")
+    rel_err = calc_rel_error(gamma, gamma_lag)
+    ax2.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
+ax2.set_xlim(x[0], x[-1])
+ax2.set_ylim(*(10.0 ** ylim))
+ax2.set_xticks(np.arange(xmin, xmax + EPSILON, 5))
+ax2.set_xticks(np.arange(xmin, xmax), minor=True)
+ax2.set_yticks(10.0 ** np.arange(*ylim, 2))
+# locmin = mpl.ticker.LogLocator(base=10.0,subs=0.1*np.arange(1,10),numticks=100)
+# ax2.yaxis.set_minor_locator(locmin)
+# ax2.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
+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")
+
+
+# Move to target
+bests = []
+N = 200
+step = 1 / (N - 1)
+a = 11 / 8
+b = 1 / 2
+x = np.linspace(step, 1 - step, N + 1)
+gamma = scipy.special.gamma(x)[:, None]
+ns = np.arange(2, 13)
+for n in ns:
+    zeros, weights = np.polynomial.laguerre.laggauss(n)
+    est = np.ceil(b + a * n)
+    targets = np.arange(max(est - 2, 0), est + 3)
+    gamma_lag = np.stack(
+        [
+            eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted")
+            for target in targets
+        ],
+        -1,
+    )
+    rel_error = np.abs(calc_rel_error(gamma, gamma_lag))
+    best = np.argmin(rel_error, -1) + targets[0]
+    bests.append(best)
+bests = np.stack(bests, 0)
+
+fig3, ax3 = plt.subplots(num=3, clear=True, constrained_layout=True, figsize=(5, 3))
+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(np.arange(0, N + 1, N // 20), minor=True)
+ax3.set_yticks(np.arange(len(ns)), ns)
+ax3.set_xlabel(r"$z$")
+ax3.set_ylabel(r"$n$")
+fig3.savefig(f"{img_path}/targets.pdf")
+
+targets = np.mean(bests, -1)
+intercept, bias = np.polyfit(ns, targets, 1)
+fig4, axs4 = plt.subplots(
+    2, num=4, sharex=True, clear=True, constrained_layout=True, figsize=(5, 4)
+)
+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[0].set_xlim(*xl)
+# axs4[0].set_title("Schätzung von Mittelwert")
+# axs4[1].set_title("Fehler")
+axs4[-1].set_xlabel(r"$z$")
+for ax in axs4:
+    ax.grid(1)
+    ax.legend()
+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)))}")
+
+# plt.show()
-- 
cgit v1.2.1


From b2f6c58490cd3517d1a813e1b51b4e2aafc945a0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Thu, 2 Jun 2022 08:09:02 +0200
Subject: Add relative error plots with shift

---
 buch/papers/laguerre/scripts/gamma_approx.py | 97 ++++++++++++++++++++++++----
 1 file changed, 84 insertions(+), 13 deletions(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

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()
-- 
cgit v1.2.1


From 85e7d741f78ca0874b42db5cfbd18f4c28a933b3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Thu, 2 Jun 2022 15:23:21 +0200
Subject: Add presentation

---
 buch/papers/laguerre/scripts/gamma_approx.py | 36 ++++++++++++++++++----------
 1 file changed, 23 insertions(+), 13 deletions(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 9c8f3ee..dd50d92 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -39,7 +39,7 @@ def find_shift(z, target):
 
 def find_optimal_shift(z, n):
     mhat = 1.34093 * n + 0.854093
-    steps = int(np.ceil(mhat - np.real(z))) - 1
+    steps = int(np.floor(mhat - np.real(z)))
     return steps
 
 
@@ -136,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
@@ -162,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
@@ -202,7 +202,7 @@ 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)
@@ -211,16 +211,16 @@ 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[0].plot(ns, targets, "x", label=r"$\overline{m}$")
+axs4[1].plot(ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \overline{m}$")
 axs4[0].set_xlim(*xl)
 # axs4[0].set_title("Schätzung von Mittelwert")
 # axs4[1].set_title("Fehler")
-axs4[-1].set_xlabel(r"$z$")
+axs4[-1].set_xlabel(r"$n$")
 for ax in axs4:
     ax.grid(1)
     ax.legend()
-# fig4.savefig(f"{img_path}/schaetzung.pgf")
+fig4.savefig(f"{img_path}/estimate.pgf")
 
 print(f"Intercept={intercept:.6g}, Bias={bias:.6g}")
 predicts = np.ceil(intercept * ns[:, None] + bias - x)
@@ -234,14 +234,14 @@ 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)
+fig5, ax5 = plt.subplots(num=5, 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}$")
+    ax5.semilogy(x, rel_error, label=f"$m={target}$", linewidth=3)
 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.semilogy(x, rel_error, "c", linestyle="dotted", label="$m^*$", linewidth=3)
 ax5.set_xlim(x[0], x[-1])
 ax5.set_ylim(5e-9, 5e-8)
 ax5.set_xlabel(r"$z$")
@@ -254,10 +254,10 @@ 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)
+fig6, ax6 = plt.subplots(num=6, 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.semilogy(x, rel_error, label="$m^*$", linewidth=3)
 ax6.set_xlim(x[0], x[-1])
 ax6.set_ylim(5e-9, 5e-8)
 ax6.set_xlabel(r"$z$")
@@ -265,4 +265,14 @@ ax6.grid(1, "both")
 ax6.legend()
 fig6.savefig(f"{img_path}/rel_error_range.pgf")
 
+N = 2001
+x = np.linspace(-5, 5, N)
+gamma = scipy.special.gamma(x)
+fig7, ax7 = plt.subplots(num=7, clear=True, constrained_layout=True)
+ax7.plot(x, gamma)
+ax7.set_xlim(x[0], x[-1])
+ax7.set_ylim(-7.5, 25)
+ax7.grid(1, "both")
+fig7.savefig(f"{img_path}/gamma.pgf")
+
 # plt.show()
-- 
cgit v1.2.1


From fac45f54d4cee5018c063b4a720695cbf3040fa9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Thu, 2 Jun 2022 15:53:49 +0200
Subject: Correct typos in presentation

---
 buch/papers/laguerre/scripts/gamma_approx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index dd50d92..857c735 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -192,7 +192,7 @@ bests = np.stack(bests, 0)
 
 fig3, ax3 = plt.subplots(num=3, clear=True, constrained_layout=True, figsize=(5, 3))
 v = ax3.imshow(bests, cmap="inferno", aspect="auto", interpolation="nearest")
-plt.colorbar(v, ax=ax3, label=r"$m$")
+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)
-- 
cgit v1.2.1


From 8fb46098cb8e42a94b8e01ecc809f536d5c7efaf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Fri, 3 Jun 2022 07:23:21 +0200
Subject: Minor tweaks of presentation

---
 buch/papers/laguerre/scripts/gamma_approx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 857c735..53ba76b 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -241,7 +241,7 @@ for target in targets:
     ax5.semilogy(x, rel_error, label=f"$m={target}$", linewidth=3)
 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, "c", linestyle="dotted", label="$m^*$", linewidth=3)
+ax5.semilogy(x, rel_error, "m", linestyle="dotted", label="$m^*$", linewidth=3)
 ax5.set_xlim(x[0], x[-1])
 ax5.set_ylim(5e-9, 5e-8)
 ax5.set_xlabel(r"$z$")
-- 
cgit v1.2.1


From fde57297b3efbef28d09a532e1b3895d2b2ad917 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Thu, 14 Jul 2022 15:03:28 +0200
Subject: Correct Makefile, add text to gamma.tex, separate python-scripts for
 each image

---
 buch/papers/laguerre/scripts/gamma_approx.py | 181 +--------------------------
 1 file changed, 5 insertions(+), 176 deletions(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 53ba76b..208f770 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -8,6 +8,7 @@ import scipy.special
 EPSILON = 1e-7
 root = str(Path(__file__).parent)
 img_path = f"{root}/../images"
+fontsize = "medium"
 
 
 def _prep_zeros_and_weights(x, w, n):
@@ -26,17 +27,6 @@ def pochhammer(z, n):
     return np.prod(z + np.arange(n))
 
 
-def find_shift(z, target):
-    factor = 1.0
-    steps = int(np.floor(target - np.real(z)))
-    zs = z + steps
-    if steps > 0:
-        factor = 1 / pochhammer(z, steps)
-    elif steps < 0:
-        factor = pochhammer(zs, -steps)
-    return zs, factor
-
-
 def find_optimal_shift(z, n):
     mhat = 1.34093 * n + 0.854093
     steps = int(np.floor(mhat - np.real(z)))
@@ -44,6 +34,7 @@ def find_optimal_shift(z, n):
 
 
 def get_shifting_factor(z, steps):
+    factor = 1.0
     if steps > 0:
         factor = 1 / pochhammer(z, steps)
     elif steps < 0:
@@ -56,7 +47,9 @@ def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
     n = len(x)
 
     z += 0j
-    z_shifted, correction_factor = find_shift(z, target)
+    steps = int(np.floor(target - np.real(z)))
+    z_shifted = z + steps
+    correction_factor = get_shifting_factor(z, steps)
 
     res = np.sum(x ** (z_shifted - 1) * w)
     res *= correction_factor
@@ -112,167 +105,3 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
 
 def calc_rel_error(x, y):
     return (y - x) / x
-
-
-# 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)
-fig, ax = plt.subplots(num=1, clear=True, constrained_layout=True, figsize=(5, 2.5))
-for n in ns:
-    gamma_lag = eval_laguerre_gamma(x, n=n)
-    rel_err = calc_rel_error(gamma, gamma_lag)
-    ax.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
-ax.set_xlim(x[0], x[-1])
-ax.set_ylim(*(10.0 ** ylim))
-ax.set_xticks(np.arange(xmin, xmax + EPSILON, 5))
-ax.set_xticks(np.arange(xmin, xmax), minor=True)
-ax.set_yticks(10.0 ** np.arange(*ylim, 2))
-ax.set_yticks(10.0 ** np.arange(*ylim, 2))
-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")
-
-
-# Mirrored
-xmin = -15
-xmax = 15
-ylim = np.array([-11, 1])
-x = np.linspace(xmin + EPSILON, xmax - EPSILON, 400)
-gamma = scipy.special.gamma(x)
-fig2, ax2 = plt.subplots(num=2, clear=True, constrained_layout=True, figsize=(5, 2.5))
-for n in ns:
-    gamma_lag = eval_laguerre_gamma(x, n=n, func="mirror")
-    rel_err = calc_rel_error(gamma, gamma_lag)
-    ax2.semilogy(x, np.abs(rel_err), label=f"$n={n}$")
-ax2.set_xlim(x[0], x[-1])
-ax2.set_ylim(*(10.0 ** ylim))
-ax2.set_xticks(np.arange(xmin, xmax + EPSILON, 5))
-ax2.set_xticks(np.arange(xmin, xmax), minor=True)
-ax2.set_yticks(10.0 ** np.arange(*ylim, 2))
-# locmin = mpl.ticker.LogLocator(base=10.0,subs=0.1*np.arange(1,10),numticks=100)
-# ax2.yaxis.set_minor_locator(locmin)
-# ax2.yaxis.set_minor_formatter(mpl.ticker.NullFormatter())
-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")
-
-
-# Move to target
-bests = []
-N = 200
-step = 1 / (N - 1)
-a = 11 / 8
-b = 1 / 2
-x = np.linspace(step, 1 - step, N + 1)
-gamma = scipy.special.gamma(x)[:, None]
-ns = np.arange(2, 13)
-for n in ns:
-    zeros, weights = np.polynomial.laguerre.laggauss(n)
-    est = np.ceil(b + a * n)
-    targets = np.arange(max(est - 2, 0), est + 3)
-    gamma_lag = np.stack(
-        [
-            eval_laguerre_gamma(x, target=target, x=zeros, w=weights, func="shifted")
-            for target in targets
-        ],
-        -1,
-    )
-    rel_error = np.abs(calc_rel_error(gamma, gamma_lag))
-    best = np.argmin(rel_error, -1) + targets[0]
-    bests.append(best)
-bests = np.stack(bests, 0)
-
-fig3, ax3 = plt.subplots(num=3, clear=True, constrained_layout=True, figsize=(5, 3))
-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)
-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)))
-ax3.set_yticklabels(ns)
-ax3.set_xlabel(r"$z$")
-ax3.set_ylabel(r"$n$")
-fig3.savefig(f"{img_path}/targets.pdf")
-
-targets = np.mean(bests, -1)
-intercept, bias = np.polyfit(ns, targets, 1)
-fig4, axs4 = plt.subplots(
-    2, num=4, sharex=True, clear=True, constrained_layout=True, figsize=(5, 4)
-)
-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"$\overline{m}$")
-axs4[1].plot(ns, ((intercept * ns + bias) - targets), "-x", label=r"$\hat{m} - \overline{m}$")
-axs4[0].set_xlim(*xl)
-# axs4[0].set_title("Schätzung von Mittelwert")
-# axs4[1].set_title("Fehler")
-axs4[-1].set_xlabel(r"$n$")
-for ax in axs4:
-    ax.grid(1)
-    ax.legend()
-fig4.savefig(f"{img_path}/estimate.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=5, 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}$", linewidth=3)
-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, "m", linestyle="dotted", label="$m^*$", linewidth=3)
-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=6, 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^*$", linewidth=3)
-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")
-
-N = 2001
-x = np.linspace(-5, 5, N)
-gamma = scipy.special.gamma(x)
-fig7, ax7 = plt.subplots(num=7, clear=True, constrained_layout=True)
-ax7.plot(x, gamma)
-ax7.set_xlim(x[0], x[-1])
-ax7.set_ylim(-7.5, 25)
-ax7.grid(1, "both")
-fig7.savefig(f"{img_path}/gamma.pgf")
-
-# plt.show()
-- 
cgit v1.2.1


From 3cb2fa354f814fa98474610dac744281285dafc6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Fri, 15 Jul 2022 11:40:55 +0200
Subject: First version of section 'Gauss Quadratur', fix to gamma_approx.py
 when z=0

---
 buch/papers/laguerre/scripts/gamma_approx.py | 18 +++++++++++++-----
 1 file changed, 13 insertions(+), 5 deletions(-)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 208f770..9f9dae7 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -1,7 +1,5 @@
 from pathlib import Path
 
-import matplotlib as mpl
-import matplotlib.pyplot as plt
 import numpy as np
 import scipy.special
 
@@ -58,6 +56,8 @@ def laguerre_gamma_shifted(z, x=None, w=None, n=8, target=11):
 
 
 def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
     x, w = _prep_zeros_and_weights(x, w, n)
     n = len(x)
 
@@ -73,6 +73,8 @@ def laguerre_gamma_opt_shifted(z, x=None, w=None, n=8):
 
 
 def laguerre_gamma_simple(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
     x, w = _prep_zeros_and_weights(x, w, n)
     z += 0j
     res = np.sum(x ** (z - 1) * w)
@@ -81,6 +83,8 @@ def laguerre_gamma_simple(z, x=None, w=None, n=8):
 
 
 def laguerre_gamma_mirror(z, x=None, w=None, n=8):
+    if z == 0.0:
+        return np.infty
     x, w = _prep_zeros_and_weights(x, w, n)
     z += 0j
     if z.real < 1e-3:
@@ -90,8 +94,8 @@ def laguerre_gamma_mirror(z, x=None, w=None, n=8):
     return laguerre_gamma_simple(z, x, w)
 
 
-def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
-    x, w = _prep_zeros_and_weights(x, w, n)
+def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):    
+    x, w = _prep_zeros_and_weights(x, w, n)    
     if func == "simple":
         f = laguerre_gamma_simple
     elif func == "mirror":
@@ -104,4 +108,8 @@ def eval_laguerre_gamma(z, x=None, w=None, n=8, func="simple", **kwargs):
 
 
 def calc_rel_error(x, y):
-    return (y - x) / x
+    mask = np.abs(x) != np.infty
+    rel_error = np.zeros_like(y)
+    rel_error[mask] = (y[mask] - x[mask]) / x[mask]
+    rel_error[~mask] = 0.0
+    return rel_error
-- 
cgit v1.2.1


From e1f5d6267540ea8dc758696fb08cb7540362cf8f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Patrik=20M=C3=BCller?= <patrik.mueller@ost.ch>
Date: Mon, 18 Jul 2022 17:34:37 +0200
Subject: First complete draft of Laguerre chapter

---
 buch/papers/laguerre/scripts/gamma_approx.py | 1 +
 1 file changed, 1 insertion(+)

(limited to 'buch/papers/laguerre/scripts/gamma_approx.py')

diff --git a/buch/papers/laguerre/scripts/gamma_approx.py b/buch/papers/laguerre/scripts/gamma_approx.py
index 9f9dae7..5b09e59 100644
--- a/buch/papers/laguerre/scripts/gamma_approx.py
+++ b/buch/papers/laguerre/scripts/gamma_approx.py
@@ -7,6 +7,7 @@ EPSILON = 1e-7
 root = str(Path(__file__).parent)
 img_path = f"{root}/../images"
 fontsize = "medium"
+cmap = "plasma"
 
 
 def _prep_zeros_and_weights(x, w, n):
-- 
cgit v1.2.1