Skip to content

Commit

Permalink
Fixed small bug and added scipy version, issue #74
Browse files Browse the repository at this point in the history
  • Loading branch information
rapoliveira committed Dec 7, 2023
1 parent 8b8e467 commit fe58801
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 18 deletions.
29 changes: 16 additions & 13 deletions data/interpolate_elliptic_integral_3.py
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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):
Expand All @@ -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)

Expand Down
16 changes: 11 additions & 5 deletions source/MulensModel/pointlens.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down

4 comments on commit fe58801

@rpoleski
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should not check the version number this way: if scipy.__version__ >= "1.9.0":. I suggest distutils.version - link

@rapoliveira
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The first answer says it is deprecated and undocumented. Isn't it better to use packaging.version.parse?
Let me know if I also should update MulensModel version number to 2.19.1.

@rpoleski
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I knew you will read more about it, then I did. The problem is that packaging is a third party package. In that case it's better to use try/except and use RGI as a first choice.
The version number is decided just before the merge.

@rapoliveira
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done!

Please sign in to comment.