From dcb2fc63fba230f65272d4d9404b58abd8977861 Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Sat, 2 Dec 2023 15:48:51 +0100 Subject: [PATCH 1/6] Replaced interp2d in data folder for issue #74, checking diff --- 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) From 8b8e4675c04e3b2805dafdc25a3e591e4531fe94 Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Mon, 4 Dec 2023 17:35:56 +0100 Subject: [PATCH 2/6] Replaced interp2d in pointlens.py, issue #74 --- data/interpolate_elliptic_integral_3.py | 21 --------------------- source/MulensModel/pointlens.py | 8 ++++++-- 2 files changed, 6 insertions(+), 23 deletions(-) diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index f6a96fe9..26d0dd79 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/data/interpolate_elliptic_integral_3.py @@ -49,27 +49,6 @@ def get_ellip(x, y): # 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): check_ = np.logspace(np.log10(x[i]), np.log10(x[i+1]), n_divide) diff --git a/source/MulensModel/pointlens.py b/source/MulensModel/pointlens.py index a187193b..ec7e36ce 100644 --- a/source/MulensModel/pointlens.py +++ b/source/MulensModel/pointlens.py @@ -4,6 +4,7 @@ from math import sin, cos, sqrt, log10 from scipy import integrate from scipy.interpolate import interp1d, interp2d +from scipy.interpolate import RegularGridInterpolator as RGI from scipy.special import ellipk, ellipe # These are complete elliptic integrals of the first and the second kind. from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3 @@ -92,7 +93,9 @@ def _read_elliptic_files(self): if line[:3] == "# Y": yy = np.array([float(t) for t in line.split()[2:]]) pp = np.loadtxt(file_3) - PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic') + # PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic') + PointLens._interpolate_3 = RGI((xx, yy), pp.T, method='cubic', + bounds_error=False) PointLens._interpolate_3_min_x = np.min(xx) PointLens._interpolate_3_max_x = np.max(xx) PointLens._interpolate_3_min_y = np.min(yy) @@ -595,7 +598,8 @@ def _get_ellip3(self, n, k): cond_4 = (k <= PointLens._interpolate_3_max_y) if cond_1 and cond_2 and cond_3 and cond_4: - return PointLens._interpolate_3(n, k)[0] + # return PointLens._interpolate_3(n, k)[0] + return float(PointLens._interpolate_3((n,k)).T) return ellip3(n, k) def get_point_lens_large_LD_integrated_magnification(self, u, gamma): From fe58801d1d729b0d60cfb0b4fa21e6bd5c73a6ce Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Thu, 7 Dec 2023 16:50:48 +0100 Subject: [PATCH 3/6] Fixed small bug and added scipy version, issue #74 --- data/interpolate_elliptic_integral_3.py | 29 ++++++++++++++----------- source/MulensModel/pointlens.py | 16 +++++++++----- 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index 26d0dd79..1eb8d105 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/data/interpolate_elliptic_integral_3.py @@ -1,11 +1,12 @@ """ Calculates interpolation tables for elliptical integral of the third kind. """ -import math +# import math import numpy as np -from math import sin, cos, sqrt -from scipy import integrate -from scipy.interpolate import interp1d, interp2d +# from math import sin, cos, sqrt +import scipy +# from scipy import integrate +from scipy.interpolate import interp2d # , interp1d from scipy.interpolate import RegularGridInterpolator as RGI from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3 @@ -46,8 +47,10 @@ def get_ellip(x, y): add_y = [] p = get_ellip(x, y) - # interp_p = interp2d(x, y, p.T, kind='cubic') - interp_rgi = RGI((x, y), p, method='cubic', bounds_error=False) + if scipy.__version__ >= "1.9.0": + interp_p = RGI((x, y), p, method='cubic', bounds_error=False) + else: + interp_p = interp2d(x, y, p.T, kind='cubic') check_x = [] for i in range(len(x)-1): @@ -60,13 +63,13 @@ def get_ellip(x, y): check_true_p = get_ellip(check_x, check_y) check_p = np.zeros((len(check_x), len(check_y))) for (ix, cx) in enumerate(check_x): - for (iy, cy) in enumerate(check_y): - if cy > cx: - 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] = float(interp_rgi((cx, cy)).T) + cond = np.array(check_y) > cx + check_p[ix, cond] = 1. + check_true_p[ix, cond] = 1. + if scipy.__version__ >= "1.9.0": + check_p[ix, ~cond] = interp_p((cx, np.array(check_y)[~cond])) + else: + check_p[ix, ~cond] = interp_p(cx, np.array(check_y)[~cond]).T[0] 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) diff --git a/source/MulensModel/pointlens.py b/source/MulensModel/pointlens.py index ec7e36ce..b084e336 100644 --- a/source/MulensModel/pointlens.py +++ b/source/MulensModel/pointlens.py @@ -2,6 +2,7 @@ import warnings import numpy as np from math import sin, cos, sqrt, log10 +import scipy from scipy import integrate from scipy.interpolate import interp1d, interp2d from scipy.interpolate import RegularGridInterpolator as RGI @@ -93,9 +94,12 @@ def _read_elliptic_files(self): if line[:3] == "# Y": yy = np.array([float(t) for t in line.split()[2:]]) pp = np.loadtxt(file_3) - # PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic') - PointLens._interpolate_3 = RGI((xx, yy), pp.T, method='cubic', - bounds_error=False) + + if scipy.__version__ >= "1.9.0": + PointLens._interpolate_3 = RGI((xx, yy), pp, method='cubic', + bounds_error=False) + else: + PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic') PointLens._interpolate_3_min_x = np.min(xx) PointLens._interpolate_3_max_x = np.max(xx) PointLens._interpolate_3_min_y = np.min(yy) @@ -598,8 +602,10 @@ def _get_ellip3(self, n, k): cond_4 = (k <= PointLens._interpolate_3_max_y) if cond_1 and cond_2 and cond_3 and cond_4: - # return PointLens._interpolate_3(n, k)[0] - return float(PointLens._interpolate_3((n,k)).T) + if scipy.__version__ >= "1.9.0": + return float(PointLens._interpolate_3((n, k)).T) + else: + return PointLens._interpolate_3(n, k)[0] return ellip3(n, k) def get_point_lens_large_LD_integrated_magnification(self, u, gamma): From 8e879bbe9e336cdb8081c095c993979cb8b6f837 Mon Sep 17 00:00:00 2001 From: "Raphael A. P. Oliveira" Date: Thu, 14 Dec 2023 19:54:16 +0100 Subject: [PATCH 4/6] Changed scipy version check to try/except, issue #74 --- data/interpolate_elliptic_integral_3.py | 10 +++++----- source/MulensModel/pointlens.py | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index 1eb8d105..86e77c8c 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/data/interpolate_elliptic_integral_3.py @@ -47,10 +47,10 @@ def get_ellip(x, y): add_y = [] p = get_ellip(x, y) - if scipy.__version__ >= "1.9.0": - interp_p = RGI((x, y), p, method='cubic', bounds_error=False) - else: - interp_p = interp2d(x, y, p.T, kind='cubic') + try: + interp_p = RGI((x, y), p, method="cubic", bounds_error=False) + except ValueError: + interp_p = interp2d(x, y, p.T, kind="cubic") check_x = [] for i in range(len(x)-1): @@ -66,7 +66,7 @@ def get_ellip(x, y): cond = np.array(check_y) > cx check_p[ix, cond] = 1. check_true_p[ix, cond] = 1. - if scipy.__version__ >= "1.9.0": + if isinstance(interp_p, RGI): check_p[ix, ~cond] = interp_p((cx, np.array(check_y)[~cond])) else: check_p[ix, ~cond] = interp_p(cx, np.array(check_y)[~cond]).T[0] diff --git a/source/MulensModel/pointlens.py b/source/MulensModel/pointlens.py index b084e336..eb541465 100644 --- a/source/MulensModel/pointlens.py +++ b/source/MulensModel/pointlens.py @@ -95,10 +95,10 @@ def _read_elliptic_files(self): yy = np.array([float(t) for t in line.split()[2:]]) pp = np.loadtxt(file_3) - if scipy.__version__ >= "1.9.0": + try: PointLens._interpolate_3 = RGI((xx, yy), pp, method='cubic', bounds_error=False) - else: + except ValueError: PointLens._interpolate_3 = interp2d(xx, yy, pp.T, kind='cubic') PointLens._interpolate_3_min_x = np.min(xx) PointLens._interpolate_3_max_x = np.max(xx) @@ -602,7 +602,7 @@ def _get_ellip3(self, n, k): cond_4 = (k <= PointLens._interpolate_3_max_y) if cond_1 and cond_2 and cond_3 and cond_4: - if scipy.__version__ >= "1.9.0": + if isinstance(PointLens._interpolate_3, RGI): return float(PointLens._interpolate_3((n, k)).T) else: return PointLens._interpolate_3(n, k)[0] From 25fd6ee2250fb80a3d4bb186ef73b8f3fd419741 Mon Sep 17 00:00:00 2001 From: radek_poleski Date: Thu, 21 Dec 2023 16:00:26 +0100 Subject: [PATCH 5/6] version update before merge --- source/MulensModel/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/MulensModel/version.py b/source/MulensModel/version.py index 326b86ff..71fae1bc 100644 --- a/source/MulensModel/version.py +++ b/source/MulensModel/version.py @@ -1 +1 @@ -__version__ = "2.19.0" +__version__ = "2.19.3" From 9db266817da169217964856467f4cc898238558e Mon Sep 17 00:00:00 2001 From: radek_poleski Date: Thu, 21 Dec 2023 16:04:06 +0100 Subject: [PATCH 6/6] flake8 cleanup in data/*.py --- data/interpolate_elliptic_integral_1_2.py | 3 --- data/interpolate_elliptic_integral_3.py | 5 ----- 2 files changed, 8 deletions(-) diff --git a/data/interpolate_elliptic_integral_1_2.py b/data/interpolate_elliptic_integral_1_2.py index f74ee9f6..61263163 100644 --- a/data/interpolate_elliptic_integral_1_2.py +++ b/data/interpolate_elliptic_integral_1_2.py @@ -2,10 +2,7 @@ Calculates interpolation tables for elliptical integral of the first and second kind. """ -import math import numpy as np -from math import sin, cos, sqrt -from scipy import integrate from scipy.interpolate import interp1d from scipy.special import ellipk, ellipe # These are complete elliptic integrals of the first and the second kind. diff --git a/data/interpolate_elliptic_integral_3.py b/data/interpolate_elliptic_integral_3.py index 86e77c8c..f38276ca 100644 --- a/data/interpolate_elliptic_integral_3.py +++ b/data/interpolate_elliptic_integral_3.py @@ -1,11 +1,7 @@ """ Calculates interpolation tables for elliptical integral of the third kind. """ -# import math import numpy as np -# from math import sin, cos, sqrt -import scipy -# from scipy import integrate from scipy.interpolate import interp2d # , interp1d from scipy.interpolate import RegularGridInterpolator as RGI from sympy.functions.special.elliptic_integrals import elliptic_pi as ellip3 @@ -24,7 +20,6 @@ def get_ellip(x, y): - p = [] z = np.zeros((len(x), len(y))) for (i, x_) in enumerate(x): for (j, y_) in enumerate(y):