import math
#given an array of Y values at consecutive integral x abscissas,
#return array of corresponding derivatives to make a natural cubic spline
def naturalSpline(ys):
    vs = [0.0] * len(ys)
    if (len(ys) < 2):
        return vs
    
    DECAY = math.sqrt(3)-2;
    endi = len(ys)-1

    # make convolutional spline
    S = 0.0;E = 0.0
    for i in range(len(Y)):
        vs[i]+=S;vs[endi-i]+=E;
        S=(S+3.0*ys[i])*DECAY;
        E=(E-3.0*ys[endi-i])*DECAY;

    #Natural Boundaries
    S2 = 6.0*(ys[1]-ys[0]) - 4.0*vs[0] - 2.0*vs[1]
    E2 = 6.0*(ys[endi-1]-ys[endi]) + 4.0*vs[endi] + 2.0*vs[endi-1]
    # A = dE2/dE = -dS2/dS,  B = dE2/dS = -dS2/dS
    A = 4.0+2.0*DECAY
    B = (4.0*DECAY+2.0)*(DECAY**(len(ys)-2))
    DEN = A*A - B*B
    S = (A*S2 + B*E2) / DEN
    E = (-A*E2 - B*S2) / DEN
    for i in range(len(ys)):
        vs[i]+=S;vs[endi-i]+=E
        S*=DECAY;E*=DECAY
    return vs

#
#Plot a different natural spline, along with its 1st and 2nd derivatives, each time you run this
#
%run plothelp.py
%matplotlib inline
import random
import numpy
Y = [random.random()*10.0+2 for _ in range(5)]
V = naturalSpline(Y)
xs = numpy.linspace(0,len(Y)-1, 1000)
plt.figure(0, figsize=(12.0,4.0))
plt.plot(xs,[hermite_interp(Y,V,x) for x in xs])
plt.plot(range(0,len(Y)),[Y[x] for x in range(0,len(Y))], "bo")
plt.figure(1, figsize=(12.0,4.0));plt.grid(True)
plt.plot(xs,[hermite_interp1(Y,V,x) for x in xs])
plt.plot(xs,[hermite_interp2(Y,V,x) for x in xs])

"""
Explanation: Convolutional Cubic Splines
$C^2$-continuous cubic splines through evenly spaced data points can be created by convolving the data points with a $C^2$-continuous piecewise cubic kernel, characterized as follows:
\begin{split}
 y(0) &= 1\
 y(x) &= 0,\text{ for all integer }x \neq 0\
 y'(x) &= sgn(x) * 3(\sqrt3-2)^{|x|},\text{ for all integer }x\
\end{split}
which implies
$$y''(x) = -6\sqrt 3(\sqrt3-2)^{|x|},\text{ for all integer }x \neq 0$$
The double-sided exponential allows the convolution to be performed extremely efficiently.
Any of the standard boundary conditions can then be applied by adjusting the start and end derivatives appropritately, and propagating the change to the rest of the curve.
The following function calculates "natural" cubic splines (with $y'' = 0$ at start and end) using this technique, which is much easiser than the way everyone is taught!
Use it however you like.
Cheers,
Matt Timmermans
End of explanation
"""


#
# Plot the kernel
#
DECAY = math.sqrt(3)-2;
vs = [3*(DECAY**x) for x in range(1,7)]
ys = [0]*len(vs) + [1] + [0]*len(vs)
vs = [-v for v in vs[::-1]] + [0.0] + vs
xs = numpy.linspace(0,len(ys)-1, 1000)
plt.figure(0, figsize=(12.0,4.0));plt.grid(True);plt.ylim([-0.2,1.1]);plt.xticks(range(-5,6))
plt.plot([x-6.0 for x in xs],[hermite_interp(ys,vs,x) for x in xs])

"""
Explanation: The Kernel
The kernel decays quickly around $x=0$, which is why cubic splines suffer from very little "ringing" -- moving one point doesn't significantly affect the curve at points far away.
End of explanation
"""


