Inserting a knot in a spline

[FabricioS]

scipy.interpolate.insert allows for inserting a knot (with an arbitrary multiplicity) in a B-spline between two knots (say u_i and u_{i+1}) without modifying the others knots, as described in "Inserting new knots into B-spline", W. Boehm, Comp.-Aided Design 12(4), pp.199-201, 1980.
I have code that relies on that property and worked well for the past few years (since 2010).

Some change may have been committed that did not respect this non-modification of knots, see attached minimal example where for a simple knot sequence (obtained with splprep) some knot is deleted (6), another one is inserted (5 with multiplicity going from 1 to 2).

Has anybody experienced that already ? A clue on the commit that introduced that ?

[FabricioS]

import numpy as np
import scipy.interpolate as si

k=3

plot = True
if plot:
    import matplotlib.pyplot as plt
    plt.cla()

# Some sample values
x = np.arange(10.)
y = np.cos(x/x.max() * 2. * np.pi).reshape(1, -1)
if plot:
    plt.plot(x, y[0], 'ro', label='Samples')


tmp = si.splprep(y, u=x, k=k, full_output=True, quiet=0, s=1e-6)
tck, u = tmp[0]
print tck[0]

if plot:
    x2 = u2 = np.linspace(u.min(), u.max(), 1024)
    y2 = si.splev(u2, tck)[0]
    plt.plot(x2, y2, 'k', lw=.5, label='Spline')
    plt.legend()

# Insert knot at u=3.5 with several multiplicity
# (until k+1=4, i.e. inserting a discontinuity possibility)
for mult in range(1, k+2):
    tck2 = si.insert(3.5, tck, mult)
    if plot:
        y3 = si.splev(u2, tck2)[0]
        plt.plot(x2, y3, label='+insert %d' % mult)
        plt.legend()

    print mult, tck2[0]

[pv]

There have been almost no code changes in splev, splprep and insert.
The last code change in insert was in 2007 (see git blame -M -C -w scipy/interpolate/fitpack/insert.f and git blame -M -C -w scipy/interpolate/fitpack.py).

Your problem may be the same as in gh-2911. The insert Fortran code is the same, but produces different results on some platforms. The issue has not been traced down, since I'm not able to reproduce it.

[pv]

You should as a first step run the Scipy test suite and see if you get the failure in gh-2911.

If not, find out which of the functions produces different results compared to before.

Moreover, try to compile older Scipy version and test with them. Do some of them produce different results? (In light of the above, this is unlikely --- if it's the test_kink problem, then it appears to be either caused by a bug in the original insert fortran routine, or a bug in some recent versions of gfortran.)

[pv]

In case it's a compiler bug, compiling with different optimization levels usually reveals this. rm -rf build; python runtests.py -g -s interpolate vs rm -rf build; python runtests.py -s interpolate. The insert.f.165t.optimized output from gfortran -Wall -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops -fdump-tree-optimized -c -o /dev/null insert.f may also be interesting.

[pv]

Closing as duplicate of gh-2911, which is fixed by gh-3673.
In case this does not fix it for you, please reopen.

[FabricioS]

Thanks Pauli (and Julian) for your responsiveness,
I hadn't the time to set up properly my git repository to perfom the test on the maintenance branch that you already solved the problem!
Last question (but maybe I should ask debian maintainers) : as it touches compiled code, as a user of the debian packaged scipy, how should I proceed to get that bugfix before the release/distribution cycle ?

[pv]

The fix is not yet merged, so to install that version, do

sudo apt-get build-dep python-scipy
sudo apt-get install python-pip
git clone https://github.com/pv/scipy-work scipy
cd scipy
git checkout origin/fitpack-alias-fix
pip install --user .

As usual, it installs under ~/.local/lib/