global continuous least square parabolic smoothing

Example

x = np.linspace(0.,7,20000)
y0 = np.sin(x)
y = sin(x) + 2*randn(20000)

plot(x,y,'x',alpha=0.8)
plot(x, y0, lw=3,alpha=0.2)
s = global_parabolic_smooth(y,2000)
plot(x,s,lw=2)
#ylim(-1,1)

global_parabolic_smooth.py

def global_parabolic_smooth(y, window=20, w=None):
    #global least square fit of piecewise parabola
    #with continuous first and zeroth derivative
    #this works but it can be optimized a lot more
    #this requires solving a huge system of linear equation
    #truncate y so that it fits
    orglen=len(y)
    y = y.copy()
    if w is None:
        w = np.ones(len(y))
        
    #do the padding and reshape y in to chunks
    q, r = divmod(len(y), window) #r is also the length of the last chunk
    nchunk = q + int(bool(r)) # round up integer division
    ncoeff = nchunk+2 # +2 for b and c; n for a's    
    
    #reshape y and w into chunks
    y.resize(window*nchunk)
    y = y.reshape(nchunk, window)
    
    w.resize(window*nchunk)#this relies on zero padding
    w=w.reshape(nchunk, window)
    
    x = np.linspace(1.0/window, 1, window)
    #1.0/window is important since the 0 is actually a but further behind
    #think about the boundary between windows
    
    #array of size nchunk each one of them is sum_i (y*w*x) for each chunk
    tmp = y*w
    yx0w = np.sum(tmp, axis=1)
    tmp = tmp*x
    yx1w = np.sum(tmp, axis=1)
    tmp = tmp*x
    yx2w = np.sum(tmp, axis=1)
    #print 'yx2w='+str(yx2w)
    #array of size nchunk each one of them is sum_i(w*x^n) for each chunk
    tmp = w
    wx0 = np.sum(w, axis=1)
    tmp = tmp*x
    wx1 = np.sum(tmp, axis=1)
    tmp = tmp*x
    wx2 = np.sum(tmp, axis=1)
    tmp = tmp*x
    wx3 = np.sum(tmp, axis=1)
    tmp = tmp*x
    wx4 = np.sum(tmp, axis=1)

    #print'wx0='+str(wx0)    
    #print'wx1='+str(wx1)
    #print'wx2='+str(wx2)
    #print'wx3='+str(wx3)
    #print'wx4='+str(wx4)

    npower = 2*2 + 1 #the partial derivative can gives x**2 which multiply by former x**2 (+1 for zeroth)
    
    aqp = np.identity(nchunk) # d/ d a_q of a_p
    
    bqp = 2*np.triu(np.ones((nchunk,nchunk)), k=1) #d/ d_a_q of b_p
    #array([[ 0.,  2.,  2.,  2.,  2.],
    #   [ 0.,  0.,  2.,  2.,  2.],
    #   [ 0.,  0.,  0.,  2.,  2.],
    #   [ 0.,  0.,  0.,  0.,  2.],
    #   [ 0.,  0.,  0.,  0.,  0.]])

    tmp1 = np.arange(nchunk)
    tmp2 = np.arange(nchunk)+1
    cqp = np.triu(2*np.add.outer(-tmp1, +tmp2)-3, k=1)
    #array([[0, 1, 3, 5, 7],
    #   [0, 0, 1, 3, 5],
    #   [0, 0, 0, 1, 3],
    #   [0, 0, 0, 0, 1],
    #   [0, 0, 0, 0, 0]])
    
    apb = np.zeros(nchunk)# d/db of a_p
    bpb = np.ones(nchunk)#d/db of b_p
    cpb = np.arange(nchunk)#d/db of c_p
    
    apc = np.zeros(nchunk)#d/dc of a_p
    bpc = np.zeros(nchunk)#d/dc of b_p
    cpc = np.ones(nchunk)#d/dc of c_p
    
    #defind S => S[0] = c, S[1] = b, S[2]=a[0], ... S[n]=a[n-2]
    dads = np.concatenate(([apc], [apb], aqp),axis=0)
    
    dbds = np.concatenate(([bpc], [bpb], bqp),axis=0)
    
    dcds = np.concatenate(([cpc], [cpb], cqp), axis=0)
    #print 'dads='+str(dads)
    #print 'dbds='+str(dbds)
    #print 'dcds='+str(dcds)
    
    #constant terms from d/da F
    B = dads.dot(yx2w) + dbds.dot(yx1w) + dcds.dot(yx0w)
    
    A = (dads*wx4).dot(dads.T) + \
        (dbds*wx3).dot(dads.T) + (dads*wx3).dot(dbds.T) + \
        (dcds*wx2).dot(dads.T) + (dbds*wx2).dot(dbds.T) + (dads*wx2).dot(dcds.T) + \
        (dcds*wx1).dot(dbds.T) + (dbds*wx1).dot(dcds.T) +\
        (dcds*wx0).dot(dcds.T)
    #print 'here'
    #print (dads*wx4).dot(dads.T)
    #print (dbds*wx3).dot(dads.T) + (dads*wx3).dot(dbds.T)
    #print (dcds*wx2).dot(dads.T) 
    #print dbds.dot(dbds.T)
    #print (dcds*wx0).dot(dcds.T)
    
    #print 'A='+str(A)
    #print 'B='+str(B)
    S = np.linalg.solve(A,B)
    #print 'S='+str(S)
    #print 'd='+str(A.dot(S)-B)
    
    a = S[2:]
    b = np.zeros(nchunk)
    b[0] = S[1]
    c = np.zeros(nchunk)
    c[0] = S[0]
    for i in range(1,nchunk):
        b[i] = 2*a[i-1] + b[i-1]
        c[i] = a[i-1]+b[i-1]+c[i-1]
    
    
    a = a.reshape((nchunk,1))
    b = b.reshape((nchunk,1))
    c = c.reshape((nchunk,1))
    #print 'a='+str(a)
    #print 'b='+str(b)
    #print 'c='+str(c)
    xx = np.tile(x,(nchunk, 1))
    ret = xx*xx*a + xx*b + c
    ret = ret.reshape(nchunk*window)
    return ret[:orglen]