MulensModel
MulensModel copied to clipboard
rpoleski
replace interp2d
In SciPy v1.10 (repeased Jan 3, 2023) interp2d() (used only in pointlens.py) is deprecated and will be removed in v1.12 (probably Jan 2024).
We have to replace it with LinearNDInterpolator or CloughTocher2DInterpolator.
Also see https://gist.github.com/ev-br/8544371b40f414b7eaf3fe6217209bff
interp2d is used in pointlens.py when reading the table with values of elliptical integrals of the 3rd kind. In order to test changes without affecting the source code, I edited data/interpolate_elliptic_integral_3.py, which generates this table.
Following the tutorial in the link, LinearNDInterpolator and CloughTocher2DInterpolator are specific for using in scattered 2D linear interpolations, i.e. for irregular grids not regularly spaced (len(x) == len(y) == len(z)).
In our case, since the grid is regularly spaced with len(x) = len(y) = 10 and len(z) = len(x)*len(y), where z is a (10,10) array obtained from the function get_ellip(x, y), the scipy function RegularGridInterpolator is more adequate.
Changing from:
from scipy.interpolate import interp1d, interp2d
interp_p = interp2d(x, y, p.T, kind='cubic')
xnew = np.arange(min(x), max(x), 1e-2)
ynew = np.arange(min(y), max(y), 1e-2)
znew = interp_p(xnew, ynew)
To:
from scipy.interpolate import RegularGridInterpolator as RGI
r = RGI((x, y), p, method='cubic', bounds_error=False)
xxnew, yynew = np.meshgrid(xnew, ynew, indexing='ij', sparse=True)
znew = r((xxnew, yynew)).T
leads to the same results, as shown in this figure.
@rpoleski please let me know if I can proceed with the RegularGridInterpolator function.
@rpoleski please check the changes in the last commit.
Some lines to check the difference between the interpolations and to make the plot above are commented.
The new function RegularGridInterpolator takes twice as long as interp2d and the interpolated values are equal up to the 10th decimal place, at minimum.
Ok, you just replaced the old function with the new one. I thought there were more changes. In that case, this should be implemented. Please make a pull request with the changes in the main code.