'''using scipy signal and numpy correlate to calculate some time series
statistics

original developer notes

see also scikits.timeseries  (movstat is partially inspired by it)
added 2009-08-29
timeseries moving stats are in c, autocorrelation similar to here
I thought I saw moving stats somewhere in python, maybe not)


TODO

moving statistics
- filters don't handle boundary conditions nicely (correctly ?)
e.g. minimum order filter uses 0 for out of bounds value
-> append and prepend with last resp. first value
- enhance for nd arrays, with axis = 0



Note: Equivalence for 1D signals
>>> np.all(signal.correlate(x,[1,1,1],'valid')==np.correlate(x,[1,1,1]))
True
>>> np.all(ndimage.filters.correlate(x,[1,1,1], origin = -1)[:-3+1]==np.correlate(x,[1,1,1]))
True

# multidimensional, but, it looks like it uses common filter across time series, no VAR
ndimage.filters.correlate(np.vstack([x,x]),np.array([[1,1,1],[0,0,0]]), origin = 1)
ndimage.filters.correlate(x,[1,1,1],origin = 1))
ndimage.filters.correlate(np.vstack([x,x]),np.array([[0.5,0.5,0.5],[0.5,0.5,0.5]]), \
origin = 1)

>>> np.all(ndimage.filters.correlate(np.vstack([x,x]),np.array([[1,1,1],[0,0,0]]), origin = 1)[0]==\
ndimage.filters.correlate(x,[1,1,1],origin = 1))
True
>>> np.all(ndimage.filters.correlate(np.vstack([x,x]),np.array([[0.5,0.5,0.5],[0.5,0.5,0.5]]), \
origin = 1)[0]==ndimage.filters.correlate(x,[1,1,1],origin = 1))


update
2009-09-06: cosmetic changes, rearrangements
'''


import numpy as np
from scipy import signal

from numpy.testing import assert_array_equal, assert_array_almost_equal

import scikits.statsmodels.api as sm


def expandarr(x,k):
    #make it work for 2D or nD with axis
    kadd = k
    if np.ndim(x) == 2:
        kadd = (kadd, np.shape(x)[1])
    return np.r_[np.ones(kadd)*x[0],x,np.ones(kadd)*x[-1]]

def movorder(x, order = 'med', windsize=3, lag='lagged'):
    '''moving order statistics
    
    Parameters
    ----------
    x : array
       time series data
    order : float or 'med', 'min', 'max'
       which order statistic to calculate
    windsize : int
       window size
    lag : 'lagged', 'centered', or 'leading'
       location of window relative to current position
       
    Returns
    -------
    filtered array
    
    
    '''

    #if windsize is even should it raise ValueError
    if lag == 'lagged':
        lead = windsize//2
    elif lag == 'centered':
        lead = 0
    elif lag == 'leading':
        lead = -windsize//2 +1
    else:
        raise ValueError
    if np.isfinite(order) == True: #if np.isnumber(order):
        ord = order   # note: ord is a builtin function
    elif order == 'med':
        ord = (windsize - 1)/2
    elif order == 'min':
        ord = 0
    elif order == 'max':
        ord = windsize - 1
    else:
        raise ValueError

    #return signal.order_filter(x,np.ones(windsize),ord)[:-lead]
    xext = expandarr(x, windsize)
    #np.r_[np.ones(windsize)*x[0],x,np.ones(windsize)*x[-1]]
    return signal.order_filter(xext,np.ones(windsize),ord)[windsize-lead:-(windsize+lead)]

def check_movorder():
    '''graphical test for movorder'''
    import matplotlib.pylab as plt
    x = np.arange(1,10)
    xo = movorder(x, order='max')
    assert_array_equal(xo, x)
    x = np.arange(10,1,-1)
    xo = movorder(x, order='min')
    assert_array_equal(xo, x)
    assert_array_equal(movorder(x, order='min', lag='centered')[:-1], x[1:])

    tt = np.linspace(0,2*np.pi,15)
    x = np.sin(tt) + 1
    xo = movorder(x, order='max')
    plt.figure()
    plt.plot(tt,x,'.-',tt,xo,'.-')
    plt.title('moving max lagged')
    xo = movorder(x, order='max', lag='centered')
    plt.figure()
    plt.plot(tt,x,'.-',tt,xo,'.-')
    plt.title('moving max centered')
    xo = movorder(x, order='max', lag='leading')
    plt.figure()
    plt.plot(tt,x,'.-',tt,xo,'.-')
    plt.title('moving max leading')

# identity filter
##>>> signal.order_filter(x,np.ones(1),0)
##array([ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.])
# median filter
##signal.medfilt(np.sin(x), kernel_size=3)
##>>> plt.figure()
##<matplotlib.figure.Figure object at 0x069BBB50>
##>>> x=np.linspace(0,3,100);plt.plot(x,np.sin(x),x,signal.medfilt(np.sin(x), kernel_size=3))

# remove old version
##def movmeanvar(x, windowsize=3, valid='same'):
##    '''
##    this should also work along axis or at least for columns
##    '''
##    n = x.shape[0] 
##    x = expandarr(x, windowsize - 1)
##    takeslice = slice(windowsize-1, n + windowsize-1)
##    avgkern = (np.ones(windowsize)/float(windowsize))
##    m = np.correlate(x, avgkern, 'same')#[takeslice]
##    print m.shape
##    print x.shape
##    xm = x - m
##    v = np.correlate(x*x, avgkern, 'same') - m**2
##    v1 = np.correlate(xm*xm, avgkern, valid) #not correct for var of window
###>>> np.correlate(xm*xm,np.array([1,1,1])/3.0,'valid')-np.correlate(xm*xm,np.array([1,1,1])/3.0,'valid')**2
##    return m[takeslice], v[takeslice], v1

def movmean(x, windowsize=3, lag='lagged'):
    '''moving window mean
    
    
    Parameters
    ----------
    x : array
       time series data
    windsize : int
       window size
    lag : 'lagged', 'centered', or 'leading'
       location of window relative to current position
       
    Returns
    -------
    mk : array
        moving mean, with same shape as x


    Notes
    -----
    for leading and lagging the data array x is extended by the closest value of the array
    
    
    '''    
    return movmoment(x, 1, windowsize=windowsize, lag=lag)

def movvar(x, windowsize=3, lag='lagged'):
    '''moving window variance
    
    
    Parameters
    ----------
    x : array
       time series data
    windsize : int
       window size
    lag : 'lagged', 'centered', or 'leading'
       location of window relative to current position
       
    Returns
    -------
    mk : array
        moving variance, with same shape as x
    
    
    '''    
    m1 = movmoment(x, 1, windowsize=windowsize, lag=lag)
    m2 = movmoment(x, 2, windowsize=windowsize, lag=lag)
    return m2 - m1*m1

def movmoment(x, k, windowsize=3, lag='lagged'):
    '''non-central moment
    
    
    Parameters
    ----------
    x : array
       time series data
    windsize : int
       window size
    lag : 'lagged', 'centered', or 'leading'
       location of window relative to current position
       
    Returns
    -------
    mk : array
        k-th moving non-central moment, with same shape as x


    Notes
    -----
    If data x is 2d, then moving moment is calculated for each
    column.
    
    '''

    windsize = windowsize
    #if windsize is even should it raise ValueError
    if lag == 'lagged':
        #lead = -0 + windsize #windsize//2
        lead = -0# + (windsize-1) + windsize//2
        sl = slice((windsize-1) or None, -2*(windsize-1) or None)
    elif lag == 'centered':
        lead = -windsize//2  #0#-1 #+ #(windsize-1)
        sl = slice((windsize-1)+windsize//2 or None, -(windsize-1)-windsize//2 or None)
    elif lag == 'leading':
        #lead = -windsize +1#+1 #+ (windsize-1)#//2 +1
        lead = -windsize +2 #-windsize//2 +1
        sl = slice(2*(windsize-1)+1+lead or None, -(2*(windsize-1)+lead)+1 or None)
    else:
        raise ValueError

    avgkern = (np.ones(windowsize)/float(windowsize))
    xext = expandarr(x, windsize-1)
    #Note: expandarr increases the array size by 2*(windsize-1)
    
    #sl = slice(2*(windsize-1)+1+lead or None, -(2*(windsize-1)+lead)+1 or None)
    print sl

    if xext.ndim == 1:
        return np.correlate(xext**k, avgkern, 'full')[sl]
        #return np.correlate(xext**k, avgkern, 'same')[windsize-lead:-(windsize+lead)]
    else:
        print xext.shape
        print avgkern[:,None].shape
        
        # try first with 2d along columns, possibly ndim with axis
        return signal.correlate(xext**k, avgkern[:,None], 'full')[sl,:]







#x=0.5**np.arange(10);xm=x-x.mean();a=np.correlate(xm,[1],'full')
#x=0.5**np.arange(3);np.correlate(x,x,'same')
##>>> x=0.5**np.arange(10);xm=x-x.mean();a=np.correlate(xm,xo,'full')
##
##>>> xo=np.ones(10);d=np.correlate(xo,xo,'full')
##>>> xo
##xo=np.ones(10);d=np.correlate(xo,xo,'full')
##>>> x=np.ones(10);xo=x-x.mean();a=np.correlate(xo,xo,'full')
##>>> xo=np.ones(10);d=np.correlate(xo,xo,'full')
##>>> d
##array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,   9.,
##         8.,   7.,   6.,   5.,   4.,   3.,   2.,   1.])


##def ccovf():
##    pass
##    #x=0.5**np.arange(10);xm=x-x.mean();a=np.correlate(xm,xo,'full')

__all__ = ['movorder', 'movmean', 'movvar', 'movmoment']

if __name__ == '__main__':    
    
    print '\ncheckin moving mean and variance'
    nobs = 10
    x = np.arange(nobs)
    ws = 3
    ave = np.array([ 0., 1/3., 1., 2., 3., 4., 5., 6., 7., 8.,
                  26/3., 9])
    va = np.array([[ 0.        ,  0.        ],
                   [ 0.22222222,  0.88888889],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.66666667,  2.66666667],
                   [ 0.22222222,  0.88888889],
                   [ 0.        ,  0.        ]])
    ave2d = np.c_[ave, 2*ave]
    print movmean(x, windowsize=ws, lag='lagged')
    print movvar(x, windowsize=ws, lag='lagged')
    print [np.var(x[i-ws:i]) for i in range(ws, nobs)]
    m1 = movmoment(x, 1, windowsize=3, lag='lagged')
    m2 = movmoment(x, 2, windowsize=3, lag='lagged')
    print m1
    print m2
    print m2 - m1*m1
    
    # this implicitly also tests moment 
    assert_array_almost_equal(va[ws-1:,0], 
                    movvar(x, windowsize=3, lag='leading'))
    assert_array_almost_equal(va[ws//2:-ws//2+1,0], 
                    movvar(x, windowsize=3, lag='centered'))
    assert_array_almost_equal(va[:-ws+1,0], 
                    movvar(x, windowsize=ws, lag='lagged'))
    
    
    
    print '\nchecking moving moment for 2d (columns only)'
    x2d = np.c_[x, 2*x]
    print movmoment(x2d, 1, windowsize=3, lag='centered')
    print movmean(x2d, windowsize=ws, lag='lagged')
    print movvar(x2d, windowsize=ws, lag='lagged')
    assert_array_almost_equal(va[ws-1:,:], 
                    movvar(x2d, windowsize=3, lag='leading'))
    assert_array_almost_equal(va[ws//2:-ws//2+1,:], 
                    movvar(x2d, windowsize=3, lag='centered'))
    assert_array_almost_equal(va[:-ws+1,:], 
                    movvar(x2d, windowsize=ws, lag='lagged'))

    assert_array_almost_equal(ave2d[ws-1:], 
                    movmoment(x2d, 1, windowsize=3, lag='leading'))
    assert_array_almost_equal(ave2d[ws//2:-ws//2+1], 
                    movmoment(x2d, 1, windowsize=3, lag='centered'))
    assert_array_almost_equal(ave2d[:-ws+1], 
                    movmean(x2d, windowsize=ws, lag='lagged'))
    
    from scipy import ndimage
    print ndimage.filters.correlate1d(x2d, np.array([1,1,1])/3., axis=0)


    #regression test check
    
    xg = np.array([  0. ,   0.1,   0.3,   0.6,   1. ,   1.5,   2.1,   2.8,   3.6,
                 4.5,   5.5,   6.5,   7.5,   8.5,   9.5,  10.5,  11.5,  12.5,
                13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,  20.5,  21.5,
                22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,
                31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,
                40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,
                49.5,  50.5,  51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,
                58.5,  59.5,  60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,
                67.5,  68.5,  69.5,  70.5,  71.5,  72.5,  73.5,  74.5,  75.5,
                76.5,  77.5,  78.5,  79.5,  80.5,  81.5,  82.5,  83.5,  84.5,
                85.5,  86.5,  87.5,  88.5,  89.5,  90.5,  91.5,  92.5,  93.5,  
                94.5])

    assert_array_almost_equal(xg, movmean(np.arange(100), 10,'lagged'))

    xd = np.array([  0.3,   0.6,   1. ,   1.5,   2.1,   2.8,   3.6,   4.5,   5.5,
                 6.5,   7.5,   8.5,   9.5,  10.5,  11.5,  12.5,  13.5,  14.5,
                15.5,  16.5,  17.5,  18.5,  19.5,  20.5,  21.5,  22.5,  23.5,
                24.5,  25.5,  26.5,  27.5,  28.5,  29.5,  30.5,  31.5,  32.5,
                33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,  40.5,  41.5,
                42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5,  50.5,
                51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,  58.5,  59.5,
                60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,  67.5,  68.5,
                69.5,  70.5,  71.5,  72.5,  73.5,  74.5,  75.5,  76.5,  77.5,
                78.5,  79.5,  80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,
                87.5,  88.5,  89.5,  90.5,  91.5,  92.5,  93.5,  94.5,  95.4,
                96.2,  96.9,  97.5,  98. ,  98.4,  98.7,  98.9,  99. ])
    assert_array_almost_equal(xd, movmean(np.arange(100), 10,'leading'))

    xc = np.array([ 1.36363636,   1.90909091,   2.54545455,   3.27272727,
                 4.09090909,   5.        ,   6.        ,   7.        ,
                 8.        ,   9.        ,  10.        ,  11.        ,
                12.        ,  13.        ,  14.        ,  15.        ,
                16.        ,  17.        ,  18.        ,  19.        ,
                20.        ,  21.        ,  22.        ,  23.        ,
                24.        ,  25.        ,  26.        ,  27.        ,
                28.        ,  29.        ,  30.        ,  31.        ,
                32.        ,  33.        ,  34.        ,  35.        ,
                36.        ,  37.        ,  38.        ,  39.        ,
                40.        ,  41.        ,  42.        ,  43.        ,
                44.        ,  45.        ,  46.        ,  47.        ,
                48.        ,  49.        ,  50.        ,  51.        ,
                52.        ,  53.        ,  54.        ,  55.        ,
                56.        ,  57.        ,  58.        ,  59.        ,
                60.        ,  61.        ,  62.        ,  63.        ,
                64.        ,  65.        ,  66.        ,  67.        ,
                68.        ,  69.        ,  70.        ,  71.        ,
                72.        ,  73.        ,  74.        ,  75.        ,
                76.        ,  77.        ,  78.        ,  79.        ,
                80.        ,  81.        ,  82.        ,  83.        ,
                84.        ,  85.        ,  86.        ,  87.        ,
                88.        ,  89.        ,  90.        ,  91.        ,
                92.        ,  93.        ,  94.        ,  94.90909091,
                95.72727273,  96.45454545,  97.09090909,  97.63636364])
    assert_array_almost_equal(xc, movmean(np.arange(100), 11,'centered'))
