From f36952afb346f9d70527dc24a3ffecd03b75e3e3 Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Sat, 2 Dec 2023 15:48:51 +0100 Subject: [PATCH] Replaced interp2d in data folder for issue #74, still checking difference --- data/interpolate_elliptic_integral_3.py | 28 +++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index 58a02653..f6a96fe9 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/data/interpolate_elliptic_integral_3.py @@ -6,6 +6,7 @@ from math import sin, cos, sqrt from scipy import integrate from scipy.interpolate import interp1d, interp2d +from scipy.interpolate import RegularGridInterpolator as RGI from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3 @@ -45,7 +46,29 @@ def get_ellip(x, y): add_y = [] p = get_ellip(x, y) - interp_p = interp2d(x, y, p.T, kind='cubic') + # interp_p = interp2d(x, y, p.T, kind='cubic') + interp_rgi = RGI((x, y), p, method='cubic', bounds_error=False) + + # Testing the change between scipy's interp2d and RGI + # xnew = np.arange(min(x), max(x), 1e-2) + # ynew = np.arange(min(y), max(y), 1e-2) + # xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True) + # znew_interp2d = interp_p(xnew, ynew) + # znew_rgi = interp_rgi((xxnew, yynew)).T + # diff = abs(znew_interp2d - znew_rgi).flatten() + # diff_max, diff_mean = diff.max(), diff.mean() # maximum: 1e-10, 1e-12 + + # Plotting to check difference + # import matplotlib.pyplot as plt + # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) + # ax1.plot(x, p.T[0, :], 'ro-', xnew, znew_interp2d[0, :], 'b-') + # im = ax2.imshow(znew) + # plt.colorbar(im, ax=ax2) + # fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) + # ax1.plot(x, p.T[0, :], 'ro-', xnew, znew_rgi[0, :], 'b-') + # im = ax2.imshow(znew_rgi) + # plt.colorbar(im, ax=ax2) + # plt.show() check_x = [] for i in range(len(x)-1): @@ -63,7 +86,8 @@ def get_ellip(x, y): check_p[ix, iy] = 1. check_true_p[ix, iy] = 1. else: - check_p[ix, iy] = interp_p(cx, cy)[0] + # check_p[ix, iy] = interp_p(cx, cy)[0] + check_p[ix, iy] = float(interp_rgi((cx, cy)).T) relative_diff_p = np.abs(check_p - check_true_p) / check_true_p index = np.unravel_index(relative_diff_p.argmax(), relative_diff_p.shape)