cvxopt spline example a bit more flexible

spline1.py

from cvxopt import lapack, solvers, matrix, mul
from cvxopt.modeling import op, variable, max, sum


#def const_spline(n=12):
   # """ Adapted from example on CVXOPT site -- Changing as little as possible """

#t = matrix(v.index.values[::4])
#y = matrix(v.values[::4])
t, y = data['t'], data['y']
m = len(t)


lbound = -1
rbound = 1
nknots = 5
order = 4 #num coeffs, so 4 for "cubic"

n = order * (nknots+1)
knots = np.linspace(lbound,rbound,nknots+1,endpoint=False)[1:]
u1, u2 = -1.0/3, 1.0/3

#get array indices by knot locations
Is = []
Is.append( [ k for k in range(m) if lbound  <= t[k] < knots[0] ])
prev_knot = knots[0]
for i,knot in enumerate(knots[1:]):
    Is.append( [ k for k in range(m) if prev_knot  <= t[k] < knot ])
    prev_knot = knot
Is.append( [ k for k in range(m) if prev_knot  <= t[k] <= rbound ])

# build system matrix
A = matrix(0.0, (m,n))
for k in range(order):
    offset = 0
    for I in Is:
        A[I,k+offset] = t[I]**k
        offset+=order

# build equality constraint

num_constraints = 3 * nknots
G = matrix(0.0, (num_constraints,n))

#1st order
for i,knot in enumerate(knots):
    G[i,list(range(order*i,order*i+8))] = 1.0,knot,knot**2,knot**3,-1.0,-knot,-knot**2,-knot**3

#2nd order
for i,knot in enumerate(knots):
    G[i+nknots,list(range(order*i,order*i+8))] = 0, 1.0,2*knot,3*knot**2, 0, -1.0, -2*knot,-3*knot**2

#3rd order
for i,knot in enumerate(knots):
    G[i+nknots*2,list(range(order*i,order*i+8))] = 0, 0, 2, 6*knot,0,0,-2,-6*knot
    


xcheb = variable(n)
c1 = G*xcheb == 0
op( max(abs(A*xcheb - y)), [c1]).solve()
xcheb = xcheb.value

nopts = 100
ts = lbound + (rbound - 1.0*lbound)/nopts * matrix(list(range(nopts)), tc='d')

#get array indices by knot locations
Is = []
Is.append( [ k for k in range(nopts) if lbound  <= ts[k] < knots[0] ])
prev_knot = knots[0]
for i,knot in enumerate(knots[1:]):
    Is.append( [ k for k in range(nopts) if prev_knot  <= ts[k] < knot ])
    prev_knot = knot
Is.append( [ k for k in range(nopts) if prev_knot  <= ts[k] <= rbound ])

ycheb = matrix(0.0, (nopts,1))

offset = 0
for I in Is:
    ycheb[I] = sum( xcheb[k+offset]*ts[I]**k for k in range(order) )
    offset+=order


pylab.plot(t, y, 'bo')
pylab.plot(ts, ycheb, '-r')