search

rpoleski / MulensModel

Microlensing Modelling package
https://rpoleski.github.io/MulensModel/
Other
56 stars 15 forks source link

replace interp2d #74

Open rpoleski opened 1 year ago

rpoleski commented 1 year ago

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

rapoliveira commented 10 months ago

@rpoleski I can do that, if it's still useful

rapoliveira commented 9 months ago

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. interp2d_after

@rpoleski please let me know if I can proceed with the RegularGridInterpolator function.

rpoleski commented 9 months ago

What changes have you made to data/interpolate_elliptic_integral_3.py?

rapoliveira commented 9 months ago

@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.

rpoleski commented 9 months ago

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.