FIX: Clean up ``scipy.signal.cubic`` deprecation

[oesteban]

Resolves: #370.

[codecov[bot]]

Codecov Report

Patch coverage: 100.00% and project coverage change: -0.06% :warning:

Comparison is base (a608927) 80.53% compared to head (2bf307f) 80.48%.

Additional details and impacted files

```diff @@ Coverage Diff @@ ## master #388 +/- ## ========================================== - Coverage 80.53% 80.48% -0.06% ========================================== Files 26 26 Lines 2296 2290 -6 Branches 286 286 ========================================== - Hits 1849 1843 -6 Misses 402 402 Partials 45 45 ``` | [Flag](https://app.codecov.io/gh/nipreps/sdcflows/pull/388/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) | Coverage Δ | | |---|---|---| | [unittests](https://app.codecov.io/gh/nipreps/sdcflows/pull/388/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) | `∅ <ø> (∅)` | | Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps#carryforward-flags-in-the-pull-request-comment) to find out more. | [Files Changed](https://app.codecov.io/gh/nipreps/sdcflows/pull/388?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps) | Coverage Δ | | |---|---|---| | [sdcflows/transform.py](https://app.codecov.io/gh/nipreps/sdcflows/pull/388?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nipreps#diff-c2RjZmxvd3MvdHJhbnNmb3JtLnB5) | `76.53% <100.00%> (-0.77%)` | :arrow_down: |

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

[effigies]

I've read and re-read the BSpline docs, and I cannot for the life of me figure out how we're supposed to get the derivatives to generate the Jacobian. Do you know how that's supposed to work?

[oesteban]

I guess https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.derivative.html

Something like BSpline.derivative().design_matrix(...) ?

[oesteban]

The coefficients are the same. It's just the weights that change. Pretty convenient :)

[effigies]

design_matrix() is a classmethod. If we find a way to write generate_design_matrix(BSpline()), then we could do generate_design_matrix(BSpline().derivative()).

[oesteban]

I see. We can initialize a BSpline because t, c and k are known when we want to calculate the jacobian. However, we will need to see how the design matrix is generated for an object because we need to calculate the tensor product.

[oesteban]

Now I'm hoping my Cython implementation is still in my local git repository...

[effigies]

Well, they do have an equivalence in https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.design_matrix.html#scipy.interpolate.BSpline.design_matrix:

from scipy.interpolate import make_interp_spline, BSpline
import numpy as np

k = 2
t = [-1, 0, 1, 2, 3, 4, 5, 6]
x = [1, 2, 3, 4]
c = np.eye(len(t) - k - 1)

design_matrix = BSpline.design_matrix(x, t, k)
design_matrix_gh = BSpline(t, c, k)(x)
np.allclose(design_matrix.toarray(), design_matrix_gh, atol=1e-14)

So we could try:

orig = BSpline(knots, np.eye(len(knots) - 3 - 1), 3)
wd0.append(scipy.sparse.csr_matrix(orig(locs).T))

if jacobian:
    deriv = orig.derivative()
    wd1.append(scipy.sparse.csr_matrix(deriv(locs).T))

[oesteban]

Well, they do have an equivalence in https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.design_matrix.html#scipy.interpolate.BSpline.design_matrix:

from scipy.interpolate import make_interp_spline, BSpline
import numpy as np

k = 2
t = [-1, 0, 1, 2, 3, 4, 5, 6]
x = [1, 2, 3, 4]
c = np.eye(len(t) - k - 1)

design_matrix = BSpline.design_matrix(x, t, k)
design_matrix_gh = BSpline(t, c, k)(x)
np.allclose(design_matrix.toarray(), design_matrix_gh, atol=1e-14)

So we could try:

orig = BSpline(knots, np.eye(len(knots) - 3 - 1), 3)
wd0.append(scipy.sparse.csr_matrix(orig(locs).T))

if jacobian:
    deriv = orig.derivative()
    wd1.append(scipy.sparse.csr_matrix(deriv(locs).T))

I don't see why this would not work.

[oesteban]

I think it's the opposite, if I understood your message correctly. If we don't drop weights, we will need to pad here.

On Fri, Aug 18, 2023, 18:42 Chris Markiewicz @.***> wrote:

@.**** commented on this pull request.

In sdcflows/transform.py https://github.com/nipreps/sdcflows/pull/388#discussion_r1298655014:

• coeffs_data.append(np.pad(
• level.get_fdata(dtype="float32"),
• ((1, 1), (1, 1), (1, 1)),
• ).reshape(-1))

If we're not removing weights when fitting, the coefficients should match the size of the design matrix, so we should not need to pad them here, correct?

— Reply to this email directly, view it on GitHub https://github.com/nipreps/sdcflows/pull/388#pullrequestreview-1584891183, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAESDRUIIPOVWY7BY2TYX7LXV6LQBANCNFSM6AAAAAA3TVE5YU . You are receiving this because you authored the thread.Message ID: @.***>

[effigies]

Math is hard, and I should probably trust the tests. If you're comfortable that we're testing everything thoroughly, I'm okay merging. I'm working on a branch off of here to add Jacobians based on my above comment.

[effigies]

Replaced by #393.