scipy.interpolate: RGI design question
Evgeni Burovski <[email protected]> Tue, 2 Jan 2024 16:33:47 +0300
| Newsgroups | gmane.comp.python.scientific.devel |
|---|---|
| Message-ID | <CAMRo0is6jxQ13PBa2V2THt+n4+6x2xsfEikFqhqrcEwZZYW-4g@mail.gmail.com> |
--===============5024049724730947840== Content-Type: multipart/alternative; boundary="000000000000fba6d9060df68c10" --000000000000fba6d9060df68c10 Content-Type: text/plain; charset="UTF-8" Hi, I'd like to solicit some opinions about a PR in scipy.interpolate. The PR itself is https://github.com/scipy/scipy/pull/19633 --- it has several moving parts; here on list I'd be most thankful for input on its (relatively) high-level design. (Of course, reviews of internals in the PR itself would be very welcome, too, if you've spare cycles!). To briefly recap the current state in the main branch --------------------------------------------------------------------- - Our current recommendation for N-dimensional interpolation on a grid is the RegularGridInterpolator class. - Its spline modes (method="cubic" etc) are very slow: they contain a python level loop of the size `num_data_points * num_dimensions` with a lot of python overhead on each iteration. As a result, it is two to three orders of magnitude slower than now-deprecated interp2d et al (see e.g. https://github.com/scipy/scipy/issues/18010). - The PR, https://github.com/scipy/scipy/pull/19633, brings performance to within a factor of 1-3 of specialized 2D `interp2d` and `*BivariateSpline`s without the loss of N-D generality, brings features, and allows further enhancements. Under the hood, it pre-constructs an instance of an NdBSpline and forwards the evaluations to it. Construction of NdBSpline is what is finicky. The current state of the PR ------------------------------------ Currently, https://github.com/scipy/scipy/pull/19633, adds new interpolation modes, (e.g. "cubic_" --- note the trailing underscore!) to existing modes (e.g. "cubic"). The existing methods are very quick to set up and slow to evaluate; the new ones are fast to evaluate, but the construction is potentially slow and might require some user control (there is large sparse linear algebra involved, so a user may need to choose between a direct or iterative solver etc). Here's how it looks: >>> from scipy.interpolate import RegularGridInterpolator as RGI >>> RGI((x, y), values, method="cubic_") # trailing underscore! is the classic fast-to-construct, slow-to-evaluate method, and >>> RGI((x, y), values, method="cubic") is a new slow-to-construct, fast-to-evaluate method. Additionally, the RGI constructor needs optional keyword arguments for sparse linalg solvers --- the PR offers reasonable defaults, but they will not cover all cases. This way, call signatures can be >>> RGI((x, y), values, method="cubic", solver=gmres, atol=1e-6) # atol is for gmres but this will raise (old methods do not accept solver kwargs): >>> RGI((x, y), values, method="cubic_", solver=gmres, atol=1e-6) # cubic_ What I'd like to ask the input for ------------------------------------------ The story with method="cubic" and "cubic_" and method-dependent kwargs is a bit messy :-). Possible options include: 1. Keep the PR as is. 2. Keep method="cubic" to be what it is in main, make new methods have underscores, "cubic_" etc. 3. As above, plus invent better names for new methods (am open to suggestions!) 4. Remove the old methods completely. The backwards compat effect is not known, I cannot confidently guarantee the behavior of the new methods in all cases. So would prefer to keep the old ones as fallbacks in some form. 5. Keep the RGI class intact, offer a parallel API for the new methods. The most straightforward way would be to have a factory function, `make_nd_spline` (a made-up name) to construct and return an NdBSpline instance. This would parallel the 1D `make_interp_spline` factory. My preference would be option 1. Pros: - improves performance in common cases with no user intervention; - offers multiple knobs to tweak things for use cases where defaults are not adequate - does not bring new names to an already large API surface. Cons: - not fully backwards compatible: some users might need to change their code to either keep strict backwards compat (add an undescore to the method invocation) or experiment with the solver/solver arguments. So, I'd greatly appreciate inputs on these possible options --- or if you've other alternatives, great, I'm all ears! If you're an RGI user, test-driving https://github.com/scipy/scipy/pull/19633 would be extremely valuable, too! Cheers, Evgeni --000000000000fba6d9060df68c10 Content-Type: text/html; charset="UTF-8" Content-Transfer-Encoding: quoted-printable <div dir=3D"ltr">Hi,<div><br></div><div>I'd like to solicit some opinio= ns about a PR in scipy.interpolate. The PR itself is <a href=3D"https://git= hub.com/scipy/scipy/pull/19633">https://github.com/scipy/scipy/pull/19633</= a> --- it has several moving parts; here on list I'd be most thankful f= or input on its (relatively) high-level design. (Of course, reviews of inte= rnals in the PR itself would be very welcome, too, if you've spare cycl= es!).</div><div><br></div><div>To briefly recap the current state in the ma= in branch</div><div>-------------------------------------------------------= --------------</div><div><br></div><div>- Our current recommendation for N-= dimensional interpolation on a grid is the RegularGridInterpolator class.</= div><div>- Its spline modes (method=3D"cubic" etc) are very slow:= =C2=A0 they contain a python level loop of the size `num_data_points * num_= dimensions` with a lot of python overhead on each iteration.</div><div>As a= result, it is two to three orders of magnitude slower than now-deprecated = interp2d et al (see e.g. <a href=3D"https://github.com/scipy/scipy/issues/1= 8010">https://github.com/scipy/scipy/issues/18010</a>).</div><div>- The PR,= =C2=A0=C2=A0<a href=3D"https://github.com/scipy/scipy/pull/19633">https://g= ithub.com/scipy/scipy/pull/19633</a>, brings performance to within a factor= of 1-3 of specialized 2D `interp2d` and `*BivariateSpline`s without the lo= ss of N-D generality, brings features, and allows further enhancements. Und= er the hood, it pre-constructs an instance of an NdBSpline and forwards the= evaluations to it. Construction of NdBSpline is what is finicky.</div><div= ><br></div><div>The current state of the PR</div><div>---------------------= ---------------</div><div><br></div><div>Currently,=C2=A0<a href=3D"https:/= /github.com/scipy/scipy/pull/19633">https://github.com/scipy/scipy/pull/196= 33</a>, adds new interpolation modes, (e.g. "cubic_" --- note the= trailing underscore!) to existing modes (e.g. "cubic").</div><di= v><br></div><div>The existing methods are very quick to set up and slow to = evaluate; the new ones are fast to evaluate, but the construction is potent= ially slow and might require some user control</div><div>(there is large sp= arse linear algebra involved, so a user may need to choose between a direct= or iterative solver etc).</div><div><br></div><div>Here's how it looks= :</div><div><br></div><div>>>> from scipy.interpolate import Regul= arGridInterpolator as RGI</div><div>>>> RGI((x, y), values, method= =3D"cubic_")=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0# trailing undersco= re!<br></div><div><br></div><div><div>is the classic fast-to-construct, slo= w-to-evaluate method, and=C2=A0</div><div><br></div></div><div>>>>= RGI((x, y), values, method=3D"cubic")</div><div><br></div><div>i= s a new slow-to-construct, fast-to-evaluate method.=C2=A0</div><div><br></d= iv><div>Additionally, the RGI constructor needs optional keyword arguments = for sparse linalg solvers --- the PR offers reasonable defaults, but they w= ill not cover all cases.</div><div>This way, call signatures can be<br></di= v><div><br></div><div>>>> RGI((x, y), values, method=3D"cubic= ", solver=3Dgmres, atol=3D1e-6)=C2=A0 =C2=A0# atol is for gmres</div><= div><br></div><div>but this will raise (old methods do not accept solver kw= args):</div><div><br></div><div>>>> RGI((x, y), values, method=3D&= quot;cubic_", solver=3Dgmres, atol=3D1e-6)=C2=A0 =C2=A0# cubic_</div><= div><br></div><div><div><br></div><div><div>What I'd like to ask the in= put for</div><div>------------------------------------------</div><div><br>= </div></div><div><div>The story with method=3D"cubic" and "c= ubic_" and method-dependent kwargs is a bit messy :-).</div><div><br><= /div></div><div>Possible options include:</div><div><br></div><div>1. Keep = the PR as is.</div><div><br></div><div>2. Keep method=3D"cubic" t= o be what it is in main, make new methods have underscores, "cubic_&qu= ot; etc.</div><div><br></div><div>3. As above, plus invent better names for= new methods (am open to suggestions!)</div><div><br></div><div>4. Remove t= he old methods completely.</div><div><br></div><div>The backwards compat ef= fect is not known, I cannot confidently guarantee the behavior of the new m= ethods in all cases. So would prefer to keep the old ones as fallbacks in s= ome form.</div><div><br></div><div>5. Keep the RGI class intact, offer a pa= rallel API for the new methods.</div><div><br></div><div>The most straightf= orward way would be to have a factory function, `make_nd_spline` (a made-up= name) to construct and return an NdBSpline instance. This would parallel t= he 1D `make_interp_spline` factory.</div><div><br></div><div>My preference = would be option 1.<br></div><div><br></div><div>Pros:</div><div>- improves = performance in common cases with no user intervention;</div><div>- offers m= ultiple knobs to tweak things for use cases where defaults are not adequate= </div><div>- does not bring new names to an already large API surface.</div= ><div><br></div><div>Cons:</div><div>-=C2=A0 not fully backwards compatible= : some users might need to change their code to either keep strict backward= s compat (add an undescore to the method invocation) or experiment with the= solver/solver arguments.<br></div><div><br></div><div>So, I'd greatly = appreciate inputs on these possible options --- or if you've other alte= rnatives, great, I'm all ears!</div><div><br></div><div>If you're a= n RGI user, test-driving=C2=A0=C2=A0<a href=3D"https://github.com/scipy/sci= py/pull/19633">https://github.com/scipy/scipy/pull/19633</a> would be extre= mely valuable, too!</div></div><div><br></div><div>Cheers,</div><div><br></= div><div>Evgeni</div><div><br></div></div> --000000000000fba6d9060df68c10-- --===============5024049724730947840== Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ SciPy-Dev mailing list -- [email protected] To unsubscribe send an email to [email protected] https://mail.python.org/mailman3/lists/scipy-dev.python.org/ Member address: [email protected] --===============5024049724730947840==--