import numpy as np
cimport numpy as np
cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def _xmvar(np.ndarray[np.double_t, ndim=1] x, int m):

    cdef double c
    cdef double v
    cdef int n = x.shape[0]
    cdef unsigned int j
    cdef unsigned int i

    c = 0
    for j in range(n - 3 * m + 1):
        v = 0
        for i in range(j, j + m):
            v += (x[i + 2 * m] - 2 * x[i + m] + x[i]) / m
        c += v * v
    return c / (2.0 * m*m * (n - 3 * m + 1))


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def _xmvar2(np.ndarray[np.double_t, ndim=1] x, np.ndarray[np.int_t, ndim=1] mm):

    cdef double z
    cdef double zi
    cdef unsigned int n = x.shape[0]
    cdef unsigned int i
    cdef unsigned int k
    cdef unsigned int m
    cdef np.ndarray[np.double_t, ndim=1] w = np.empty([n + 1, ], np.double)
    cdef np.ndarray[np.double_t, ndim=1] mvar = np.empty([mm.shape[0], ], np.double)

    mvar.fill(np.nan)

    w[0] = 0.0
    for i in range(n):
        w[i + 1] = w[i] + x[i]

    for k, m in enumerate(mm):
        if 3 * m > n:
            continue
        z = 0
        for i in range(n - 3 * m + 1):
            zi = (w[i + 3 * m] - 3 * w[i + 2 * m] + 3 * w[i + m] - w[i]) / m
            z += zi * zi
        mvar[k] = z / (2.0 * m*m * (n - 3 * m + 1))
    
    return mvar


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def _xavar(double[:] x, int[:] mm):

    cdef double z
    cdef double zi
    cdef unsigned int n = x.shape[0]
    cdef unsigned int i
    cdef unsigned int k
    cdef unsigned int m
    cdef int N

    r = np.empty([mm.shape[0], ], np.double)
    r.fill(np.nan)
    cdef double[::1] avar = r

    for k, m in enumerate(mm):    	
        if m < 1:
            continue
        N = (n // m) - 2
        if N < 1:
            continue
        z = 0
        for i in range(N):
            zi = x[i * m] - 2 * x[(i + 1) * m] + x[(i + 2) * m]
            z += (zi * zi)
        avar[k] = z / (2.0 * m*m * N)

    return r