RFC: FITPACK analog in scipy.interpolate

--===============4418117799611547978==
Content-Type: multipart/alternative; boundary="0000000000007dbfac06107d6c2d"

--0000000000007dbfac06107d6c2d
Content-Type: text/plain; charset="UTF-8"

Hi,

I'd be grateful for suggestions/opinions about the API of
scipy.interpolate. We currently rely on the FITPACK library to construct
smoothing splines. The user-facing API is `splrep` and `UnivariateSpline`.
There's a project to replicate/rework what FITPACK is doing, with two goals:
- bring the equivalent functionality to CuPy
- break the FITPACK monolith to make it easier to add enhancements to ---
see, for instance https://github.com/scipy/scipy/issues/2579 for feature
request going into its second decade of shelf life.

Specifically, here's a "stack" of three PRs aiming to replicate the spline
smoothing of `splrep` and `UnivariateSpline`. The split is mainly to keep
diff size manageable:

- https://github.com/scipy/scipy/pull/19753 brings in linalg primitives (a
version of the banded QR factorization)
- https://github.com/scipy/scipy/pull/19873 adds a generator to construct
successive knot vectors
- https://github.com/scipy/scipy/pull/19970 includes the previous two, and
adds an equivalent of `splrep`.

Detailed code reviews of these PRs would of course be greatly appreciated.
What I'd like to discuss here is the end API. The current state is this. We
are adding two objects:

1. `generate_knots(x, y, k, s)` is a python generator which yields
successive knot vectors for a given smoothing parameter `s`. The generator
is finite, and the last value is what `splrep` would return (*). The use of
the generator comes from https://github.com/scipy/scipy/issues/2579, only
here we yield knot vectors, not splines. This way, we can decouple
constructing the knots from smoothing criteria.

The generator has two stopping criteria: either required smoothness has
been reached (the deviation of the LSQ spline from data is less than the
input `s`), or the desired number of knots have been placed. The latter
functionality is present in the Fortran FITPACK, but effectively disabled
in the python-to-fortran wrappers.

(*) fine print: the knots are *almost* equivalent to what FITPACK returns.
The core procedure is to split the data into some intervals, and place new
knots into intervals with the maximum deviation. Now, if you have two
intervals with nearly identical deviation, which of them is "maximal" boils
down to numerical noise --- and this one may differ in Fortran and Python.

2. `make_splrep(x, y, k, s, t=None) -> BSpline` function. It accepts the
same parameters as `splrep`, and returns a BSpline instance. If the knots,
t, are not provided, it calls `generate_knots`, so we get the same result
as `splrep(x, y, k=k, s=s)`.
If knots are specified, it just uses them to construct the smoothing spline
with the same procedure that FITPACK uses.

So the first question is the new `make_splrep` name. Because of (*) above,
just calling the new function `splrep` has backwards compat issues.
Therefore, a new function, and the plan is to declare `splrep` and
`UnivariateSpline` legacy and leave them be, at least for the time being.

An alternative is, of course, to make a gradual transition along the lines
of https://github.com/scipy/scipy/issues/19896#issuecomment-1898608008 ---
hide the change behind a variable or an env var to switch "backends", and
make the change from the legacy Fortran backend to the new one over several
releases.

Thoughts?

Evgeni

P.S. If you are using `splrep` or `UnivariateSpline`, it'd be super helpful
if you could try `make_splrep` on your data. I would be very very
interested in how it fares, whether the results are identical, and if they
are not, how bad/serious the difference is.

--0000000000007dbfac06107d6c2d
Content-Type: text/html; charset="UTF-8"
Content-Transfer-Encoding: quoted-printable

<div dir=3D"ltr">Hi,<div><br></div><div>I&#39;d be grateful for suggestions=
/opinions about the API of scipy.interpolate. We currently rely on the FITP=
ACK library to construct smoothing splines. The user-facing API is `splrep`=
 and `UnivariateSpline`. There&#39;s a project to replicate/rework what FIT=
PACK is doing, with two goals:</div><div>- bring the equivalent functionali=
ty to CuPy</div><div>- break the FITPACK monolith to make it easier to add =
enhancements to --- see, for instance=C2=A0<a href=3D"https://github.com/sc=
ipy/scipy/issues/2579">https://github.com/scipy/scipy/issues/2579</a> for f=
eature request going into its second decade of shelf life.</div><div><br></=
div><div>Specifically, here&#39;s a &quot;stack&quot; of three PRs aiming t=
o replicate the spline smoothing of `splrep` and `UnivariateSpline`. The sp=
lit is mainly to keep diff size manageable:<br></div><div><br></div><div>-=
=C2=A0<a href=3D"https://github.com/scipy/scipy/pull/19753">https://github.=
com/scipy/scipy/pull/19753</a> brings in linalg primitives (a version of th=
e banded QR factorization)<br></div><div>-=C2=A0<a href=3D"https://github.c=
om/scipy/scipy/pull/19873">https://github.com/scipy/scipy/pull/19873</a> ad=
ds a generator to construct successive knot vectors</div><div>-=C2=A0<a hre=
f=3D"https://github.com/scipy/scipy/pull/19970">https://github.com/scipy/sc=
ipy/pull/19970</a> includes the previous two, and adds an equivalent of `sp=
lrep`.</div><div><br></div><div>Detailed code reviews of these PRs would of=
 course be greatly appreciated. What I&#39;d like to discuss here is the en=
d API. The current state is this. We are adding two objects:</div><div><br>=
</div><div>1. `generate_knots(x, y, k, s)` is a python generator which yiel=
ds successive knot vectors for a given smoothing parameter `s`. The generat=
or is finite, and the last value is what `splrep` would return (*). The use=
 of the generator comes from=C2=A0<a href=3D"https://github.com/scipy/scipy=
/issues/2579">https://github.com/scipy/scipy/issues/2579</a>, only here we =
yield knot vectors, not splines. This way, we can decouple constructing the=
 knots from smoothing criteria.</div><div><br></div><div>The generator has =
two stopping criteria: either required smoothness has been reached (the dev=
iation of the LSQ spline from data is less than the input `s`), or the desi=
red number of knots have been placed. The latter functionality is present i=
n the Fortran FITPACK, but effectively disabled in the python-to-fortran wr=
appers.</div><div><br></div><div><div>(*) fine print: the knots are *almost=
* equivalent to what FITPACK returns. The core procedure is to split the da=
ta into some intervals, and place new knots into intervals with the maximum=
 deviation. Now, if you have two intervals with nearly identical deviation,=
 which of them is &quot;maximal&quot; boils down to numerical noise --- and=
 this one may differ in Fortran and Python.</div><div><br></div></div><div>=
2. `make_splrep(x, y, k, s, t=3DNone) -&gt; BSpline` function. It accepts t=
he same parameters as `splrep`, and returns a BSpline instance. If the knot=
s, t, are not provided, it calls `generate_knots`, so we get the same resul=
t as `splrep(x, y, k=3Dk, s=3Ds)`.=C2=A0</div><div>If knots are specified, =
it just uses them to construct the smoothing spline with the same procedure=
 that FITPACK uses.</div><div><br></div><div>So the first question is the n=
ew `make_splrep` name. Because of (*) above, just calling the new function =
`splrep` has backwards compat issues. Therefore, a new function, and the pl=
an is to declare `splrep` and `UnivariateSpline` legacy and leave them be, =
at least for the time being.</div><div><br></div><div>An alternative is, of=
 course, to make a gradual transition along the lines of=C2=A0<a href=3D"ht=
tps://github.com/scipy/scipy/issues/19896#issuecomment-1898608008">https://=
github.com/scipy/scipy/issues/19896#issuecomment-1898608008</a> --- hide th=
e change behind a variable or an env var to switch &quot;backends&quot;, an=
d make the change from the legacy Fortran backend to the new one over sever=
al releases.</div><div><br></div><div>Thoughts?</div><div><br></div><div>Ev=
geni</div><div><br></div><div>P.S. If you are using `splrep` or `Univariate=
Spline`, it&#39;d be super helpful if you could try `make_splrep` on your d=
ata. I would be very very interested in how it fares, whether the results a=
re identical, and if they are not, how bad/serious the difference is.</div>=
<div><br></div></div>

--0000000000007dbfac06107d6c2d--

--===============4418117799611547978==
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]

--===============4418117799611547978==--