def coffee_pred(step):
    n=1290-step
    fea=20
    cof=1
    from sklearn import svm
    from scipy.linalg import hankel        
    import matplotlib.pyplot as plt
    from sklearn import metrics
    import math
    import numpy as np
    #input data from csv file
    path='/Users/royyang/Desktop/time_series_forecasting/csv_files/coffee_ls.txt'
    path1 = '/Users/royyang/Desktop/time_series_forecasting/csv_files/coffee_ls_nor.txt'
    result_tem=[]
    with open(path) as f:
        next(f)
        for line in f:
            item=line.replace('\n','').split(' ')
            result_tem.append(float(item[1]))
    mean = np.mean(result_tem)
    sd = np.std(result_tem)
    result = (result_tem-mean)/sd
    #form hankel matrix
    X=hankel(result[0:-fea-step+1], result[-1-fea:-1])
    y=result[fea+step-1:]
    #split data into training and testing
    Xtrain=X[:n]
    ytrain=y[:n]
    Xtest=X[n:]
    ytest=y[n:]
    # linear kernal
    svr_lin1 = svm.SVR(kernel='linear', C=cof)
    y_lin1 = svr_lin1.fit(Xtrain, ytrain)
    
    return int(y_lin1.predict(result[-1-fea:-1])*sd+mean)

def firls(m, bands, desired, weight=None):

    if weight==None : weight = ones(len(bands)/2)
    bands, desired, weight = array(bands), array(desired), array(weight)

    #if not desired[-1] == 0 and bands[-1] == 1 and m % 2 == 1:
    if m % 2 == 1:
        m = m + 1

    M = m/2
    w = kron(weight, [-1,1])
    omega = bands * pi
    i1 = arange(1,M+1)

    # generate the matrix q
    # as illustrated in the above-cited reference, the matrix can be
    # expressed as the sum of a hankel and toeplitz matrix. a factor of
    # 1/2 has been dropped and the final filter hficients multiplied
    # by 2 to compensate.
    cos_ints = append(omega, sin(mat(arange(1,m+1)).T*mat(omega))).reshape((-1,omega.shape[0]))
    q = append(1, 1.0/arange(1.0,m+1)) * array(mat(cos_ints) * mat(w).T).T[0]
    q = toeplitz(q[:M+1]) + hankel(q[:M+1], q[M : ])

    # the vector b is derived from solving the integral:
    #
    #           _ w
    #          /   2
    #  b  =   /       w(w) d(w) cos(kw) dw
    #   k    /    w
    #       -      1
    #
    # since we assume that w(w) is constant over each band (if not, the
    # computation of q above would be considerably more complex), but
    # d(w) is allowed to be a linear function, in general the function
    # w(w) d(w) is linear. the computations below are derived from the
    # fact that:
    #     _
    #    /                          a              ax + b
    #   /   (ax + b) cos(nx) dx =  --- cos (nx) +  ------ sin(nx)
    #  /                             2                n
    # -                             n
    #


    enum = append(omega[::2]**2 - omega[1::2]**2, cos(mat(i1).T * mat(omega[1::2])) - cos(mat(i1).T * mat(omega[::2]))).flatten()
    deno = mat(append(2, i1)).T * mat(omega[1::2] - omega[::2])
    cos_ints2 = enum.reshape(deno.shape)/array(deno)

    d = zeros_like(desired)
    d[::2]  = -weight * desired[::2]
    d[1::2] =  weight * desired[1::2]

    b = append(1, 1.0/i1) * array(mat(kron (cos_ints2, [1, 1]) + cos_ints[:M+1,:]) * mat(d).T)[:,0]

    # having computed the components q and b of the  matrix equation,
    # solve for the filter hficients.
    a = (array(inv(q)*mat(b).T).T)[0]
    h = append( a[:0:-1], append(2*a[0],  a[1:]))

    return h

 def update(self, x_n, d_n):
     
     # Update the internal buffers
     self.n += 1
     slot = self.L - ((self.n-1) % self.L) - 1
     self.x[slot] = x_n
     self.d[slot] = d_n
     
     # Block update
     if self.n % self.L == 0:
         
         # block-update parameters
         X = la.hankel(self.x[:self.L],r=self.x[self.L-1:])
         Lambda_diag = (self.lmbd**np.arange(self.L))
         lmbd_L = self.lmbd**self.L
         
         alpha = self.d - np.dot(X, self.w)
         
         pi = np.dot(self.P, X.T)
         g = np.linalg.solve((lmbd_L*np.diag(1/Lambda_diag) + np.dot(X,pi)).T, pi.T).T
         
         self.w += np.dot(g, alpha)
         
         self.P = self.P/lmbd_L - np.dot(g, pi.T)/lmbd_L
         
         # Remember a few values
         self.x[-self.length+1:] = self.x[0:self.length-1]

    def THinv5(self,phi,K,M,eps):
        """
         function G=THinv5(phi,K,M,eps)
         %
         %%%% Implements fast (complexity O(M*K^2))
         %%%% computation of the following piece of code:
         %
         %C=[];
         %for im=1:M 
         %  A=toeplitz(phi(1:K,im),phi(1:K,im)')+hankel(phi(1:K,im),phi(K:2*K-1,im)')+eps(im)*eye(K);
         %  C=[C inv(A)];
         %end  
         %
         % DEFAULT PARAMETERS: M=2; phi=randn(2*K-1,M); eps=randn(1,2);
         %   SIZE of phi SHOULD BE (2*K-1,M).
         %   SIZE of eps SHOULD BE (1,M).
        """

        # %C=[];
        C = []
        # %for im=1:M 
        # %  A=toeplitz(phi(1:K,im),phi(1:K,im)')+hankel(phi(1:K,im),phi(K:2*K-1,im)')+eps(im)*eye(K);
        # %  C=[C inv(A)];
        # %end
        for im in xrange(M):
            A = (toeplitz(phi[:K,im],phi[:K,im].T) + 
                 hankel(phi[:K,im],phi[K-1:2*K,im].T) +
                 eps[im]*np.eye(K))
            C.append(np.linalg.inv(A))

        return np.concatenate(C,axis=1)

def _firls_even(numtaps, bands, desired):
    
    # This function implements an algorithm similar to the one of the
    # SciPy `firls` function, but for even-length filters rather than
    # odd-length ones. See paper notes entitled "Least squares FIR
    # filter design for even N" for derivation. The derivation is
    # similar to that of Ivan Selesnick's "Linear-Phase FIR Filter
    # Design By Least Squares" (available online at
    # http://cnx.org/contents/[email protected]),
    # with due alteration of detail for even filter lengths.
    
    bands.shape = (-1, 2)
    desired.shape = (-1, 2)
    weights = np.ones(len(desired))
    M = int(numtaps / 2)
    
    # Compute M x M matrix Q (actually twice Q).
    n = np.arange(numtaps)[:, np.newaxis, np.newaxis]
    q = np.dot(np.diff(np.sinc(bands * n) * bands, axis=2)[:, :, 0], weights)
    Q1 = linalg.toeplitz(q[:M])
    Q2 = linalg.hankel(q[1:M + 1], q[M:])
    Q = Q1 + Q2
    
    # Compute M-vector b.
    k = np.arange(M) + .5
    e = bands[1]
    b = np.diff(e * np.sinc(e * k[:, np.newaxis])).reshape(-1)
    
    # Compute a (actually half a).
    a = np.dot(linalg.pinv(Q), b)
    
    # Compute h.
    h = np.concatenate((np.flipud(a), a))
    
    return h

def MomsFromHankelMoms (hm):
    """
    Returns the raw moments given the Hankel moments.
    
    The raw moments are: `m_i=E(\mathcal{X}^i)`
    
    The ith Hankel moment is the determinant of matrix 
    `\Delta_{i/2}`, if i is even, 
    and it is the determinant of `\Delta^{(1)}_{(i+1)/2}`, 
    if i is odd. For the definition of matrices `\Delta`
    and `\Delta^{(1)}` see [1]_.
       
    Parameters
    ----------
    hm : vector of doubles
        The list of Hankel moments (starting with the first
        moment)
        
    Returns
    -------
    m : vector of doubles
        The list of raw moments
    
    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Stieltjes_moment_problem
    """
    m = [hm[0]]
    for i in range(1,len(hm)):
        if i%2 == 0:
            N = i//2 + 1
            H = la.hankel(m[0:N], m[N-1:2*N-2] + [0])
        else:
            N = (i+1) // 2 + 1
            H = la.hankel([1] + m[0:N-1], m[N-2:2*N-3] + [0])
        h = hm[i]
        rH = np.delete(H, N-1, 0)
        for j in range(N):
            cofactor = (-1)**(N+j-1) * la.det(np.delete(rH, j, 1))
            if j<N-1:
                h -= cofactor * H[N-1,j]
            else:
                m.append(h/cofactor)
    return m

def CheckMoments (m, prec=1e-14):
    """
    Checks if the given moment sequence is valid in the sense
    that it belongs to a distribution with support (0,inf).
    
    This procedure checks the determinant of `\Delta_n`
    and `\Delta_n^{(1)}` according to [1]_.
    
    Parameters
    ----------
    m : list of doubles, length 2N+1
        The (raw) moments to check 
        (starts with the first moment).
        Its length must be odd.
    prec : double, optional
        Entries with absolute value less than prec are 
        considered to be zeros. The default value is 1e-14.
        
    Returns
    -------
    r : bool
        The result of the check.
    
    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Stieltjes_moment_problem
    """

    if butools.checkInput and len(m)%2==0:
        raise Exception("CheckMoments: the number of moments must be odd!")
    
    m = [1.0] + m
    N = math.floor(len(m)/2)-1
    
    for n in range(N+1):
        H = la.hankel(m[0:n+1], m[n:2*n+1])
        H0 = la.hankel(m[1:n+2], m[n+1:2*n+2])
        if la.det(H)<-prec or la.det(H0)<-butools.checkPrecision:
            return False
    return True

    def THinv5(self,phi,K,M,eps):
        # %C=[];
        C = []
        # %for im=1:M 
        # %  A=toeplitz(phi(1:K,im),phi(1:K,im)')+hankel(phi(1:K,im),phi(K:2*K-1,im)')+eps(im)*eye(K);
        # %  C=[C inv(A)];
        # %end
        for im in range(M):
            A = (toeplitz(phi[:K,im],phi[:K,im].T) + 
                 hankel(phi[:K,im],phi[K-1:2*K,im].T) +
                 eps[im]*np.eye(K))
            C.append(np.linalg.inv(A))

        return np.concatenate(C,axis=1)

def HankelMomsFromMoms (m):
    """
    Returns the Hankel moments given the raw moments.
    
    The raw moments are: `m_i=E(\mathcal{X}^i)`
    
    The ith Hankel moment is the determinant of matrix 
    `\Delta_{i/2}`, if i is even, 
    and it is the determinant of `\Delta^{(1)}_{(i+1)/2}`, 
    if i is odd. For the definition of matrices `\Delta`
    and `\Delta^{(1)}` see [1]_.
       
    Parameters
    ----------
    m : vector of doubles
        The list of raw moments (starting with the first
        moment)
        
    Returns
    -------
    hm : vector of doubles
        The list of Hankel moments
    
    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Stieltjes_moment_problem
    """
    hm = []
    for i in range(len(m)):
        if i%2 == 0:
            N = i//2 + 1
            H = la.hankel(m[0:N], m[N-1:2*N-1])
        else:
            N = (i+1) // 2 + 1
            H = la.hankel([1] + m[0:N-1], m[N-2:2*N-2])
        hm.append (la.det(H))
    return hm

 def PronyExpFit(self, deg, x, y):
     '''
     Construct a sum of exponentials fit to a given time-sequence y by
     using prony's method
     
     input:
         [deg]: int, number of exponentials
         [x]: numpy array, sequence of regularly spaced points at which y is evaluated
         [y]: numpy array, sequence
         
     output:
         [a]: numpy array, exponential coefficient
         [c]: numpy array, exponentail magnitudes
         [rms]: float, root mean square error of the data
     '''
     # stepsize
     h = x[1] - x[0]
     #Build matrix
     A = la.hankel(y[:-deg],y[-deg-1:])
     a = -A[:,:deg]    
     b = A[:,deg]
     #Solve it
     s = la.lstsq(a,b)[0]
     #Solve polynomial
     p = np.flipud(np.hstack((s,1)))
     u = np.roots(p)
     #Only keep roots in unit circle
     inds = np.where(np.logical_and((np.abs(u) < 1.), \
             np.logical_not(np.logical_and(np.imag(u) == 0., np.real(u) <= 0.))))[0]
     u = u[inds]
     #Calc exponential factors
     a = np.log(u)/h
     #Build power matrix
     A = np.power((np.ones((len(y),1))*u[:,None].T),np.arange(len(y))[:,None]*np.ones((1,len(inds))))
     #solve it 
     f = la.lstsq(A,y)[0]
     #calc amplitudes 
     c = f/np.exp(a*x[0])
     #build x, approx and calc rms
     approx = self.sumExp(x, a, c).real
     rms = np.sqrt(((approx-y)**2).sum() / len(y))
     return a, c, rms

def fit_direct(x, y, F=0, weighted=True, _weights=None):
    """Fit a Gaussian to the given data.

    Returns a fit so that y ~ gauss(x, A, mu, sigma)

    Parameters
    ----------
    x : ndarray
        Sampling positions.
    y : ndarray
        Sampled values.
    F : float
        Ignore values of y <= F.
    weighted : bool
        Whether to use weighted least squares.  If True, weigh
        the error function by y, ensuring that small values
        has less influence on the outcome.

    Additional Parameters
    ---------------------
    _weights : ndarray
        Weights used in weighted least squares.  For internal use
        by fit_iterative.

    Returns
    -------
    A : float
        Amplitude.
    mu : float
        Mean.
    std : float
        Standard deviation.

    """
    mask = (y > F)
    x = x[mask]
    y = y[mask]

    if _weights is None:
        _weights = y
    else:
        _weights = _weights[mask]

    # We do not want to risk working with negative values
    np.clip(y, 1e-10, np.inf, y)

    e = np.ones(len(x))
    if weighted:
        e = e * (_weights**2)
    
    v = (np.sum(np.vander(x, 5) * e[:, None], axis=0))[::-1]
    A = v[sl.hankel([0, 1, 2], [2, 3, 4])]

    ly = e * np.log(y)
    ls = np.sum(ly)
    x_ls = np.sum(ly * x)
    xx_ls = np.sum(ly * x**2)
    B = np.array([ls, x_ls, xx_ls])

    (a, b, c), res, rank, s = np.linalg.lstsq(A, B)

    A = np.exp(a - (b**2 / (4 * c)))
    mu = -b / (2 * c)
    sigma = sp.sqrt(-1 / (2 * c))

    return A, mu, sigma

    def calc_ss_radiation(self, max_order=10, r2_thresh=0.95):
        '''Function to calculate the state space reailization of the wave
        radiation IRF.

        Args:
            max_order : int
                Maximum order of the state space reailization fit
            r2_thresh : float
                The R^2 threshold used for the fit. THe value must be between 0
                and 1

        Returns:
            No variables are directily returned by thi function. The following
            internal variables are calculated:

            Ass :
                time-invariant state matrix
            Bss :
                time-invariant input matrix
            Css :
                time-invariant output matrix
            Dss :
                time-invariant feedthrough matrix
            k_ss_est :
                Impusle response function as cacluated from state space
                approximation
            status :
                status of the realization, 0 - zero hydrodynamic oefficients, 1
                - state space realization meets R2 thresholdm 2 - state space
                realization does not meet R2 threshold and at ss_max limit

        Examples:
        '''
        dt = self.rd.irf.t[2] - self.rd.irf.t[1]
        r2bt = np.zeros([self.am.inf.shape[0], self.am.inf.shape[0], self.rd.irf.t.size])
        k_ss_est = np.zeros(self.rd.irf.t.size)
        self.rd.ss.irk_bss = np.zeros([self.am.inf.shape[0], self.am.inf.shape[0], self.rd.irf.t.size])
        self.rd.ss.A = np.zeros([6, self.am.inf.shape[1], max_order, max_order])
        self.rd.ss.B = np.zeros([6, self.am.inf.shape[1], max_order, 1])
        self.rd.ss.C = np.zeros([6, self.am.inf.shape[1], 1, max_order])
        self.rd.ss.D = np.zeros([6, self.am.inf.shape[1], 1])
        self.rd.ss.irk_bss = np.zeros([6, self.am.inf.shape[1], self.rd.irf.t.size])
        self.rd.ss.rad_conv = np.zeros([6, self.am.inf.shape[1]])
        self.rd.ss.it = np.zeros([6, self.am.inf.shape[1]])
        self.rd.ss.r2t = np.zeros([6, self.am.inf.shape[1]])

        pbar = ProgressBar(widgets=['Radiation damping state space realization for ' + self.name + ':', Percentage(), Bar()], maxval=self.am.inf.shape[0] * self.am.inf.shape[1]).start()
        count = 0
        for i in xrange(self.am.inf.shape[0]):

          for j in xrange(self.am.inf.shape[1]):

            r2bt = np.linalg.norm(
                self.rd.irf.K[i, j, :] - self.rd.irf.K.mean(axis=2)[i, j])

            ss = 2  # Initial state space order

            if r2bt != 0.0:
              while True:

                # Perform Hankel Singular Value Decomposition
                y = dt * self.rd.irf.K[i, j, :]
                h = hankel(y[1::])
                u, svh, v = np.linalg.svd(h)

                u1 = u[0:self.rd.irf.t.size - 2, 0:ss]
                v1 = v.T[0:self.rd.irf.t.size - 2, 0:ss]
                u2 = u[1:self.rd.irf.t.size - 1, 0:ss]
                sqs = np.sqrt(svh[0:ss].reshape(ss, 1))
                invss = 1 / sqs
                ubar = np.dot(u1.T, u2)

                a = ubar * np.dot(invss, sqs.T)
                b = v1[0, :].reshape(ss, 1) * sqs
                c = u1[0, :].reshape(1, ss) * sqs.T
                d = y[0]

                CoeA = dt / 2
                CoeB = 1
                CoeC = -CoeA
                CoeD = 1

                # (T/2*I + T/2*A)^{-1}         = 2/T(I + A)^{-1}
                iidd = np.linalg.inv(CoeA * np.eye(ss) - CoeC * a)

                # (A-I)2/T(I + A)^{-1}         = 2/T(A-I)(I + A)^{-1}
                ac = np.dot(CoeB * a - CoeD * np.eye(ss), iidd)
                # (T/2+T/2)*2/T(I + A)^{-1}B   = 2(I + A)^{-1}B
                bc = (CoeA * CoeB - CoeC * CoeD) * np.dot(iidd, b)
                # C * 2/T(I + A)^{-1}          = 2/T(I + A)^{-1}
                cc = np.dot(c, iidd)
                # D - T/2C (2/T(I + A)^{-1})B  = D - C(I + A)^{-1})B
                dc = d + CoeC * np.dot(np.dot(c, iidd), b)

                for jj in xrange(self.rd.irf.t.size):

                  # Calculate impulse response function from state space
                  # approximation
                  k_ss_est[jj] = np.dot(np.dot(cc, expm(ac * dt * jj)), bc)

#.........这里部分代码省略.........

 def time_hankel(self, size):
     sl.hankel(self.x)

def bicoherencex(w, x, y, nfft=None, wind=None, nsamp=None, overlap=None):
  """
  Direct (FD) method for estimating cross-bicoherence
  Parameters:
    w,x,y - data vector or time-series
          - should have identical dimensions
    nfft - fft length [default = power of two > nsamp]
           actual size used is power of two greater than 'nsamp'
    wind - specifies the time-domain window to be applied to each
           data segment; should be of length 'segsamp' (see below);
      otherwise, the default Hanning window is used.
    segsamp - samples per segment [default: such that we have 8 segments]
            - if x is a matrix, segsamp is set to the number of rows
    overlap - percentage overlap, 0 to 99  [default = 50]
            - if y is a matrix, overlap is set to 0.

  Output:
    bic     - estimated cross-bicoherence: an nfft x nfft array, with
              origin at center, and axes pointing down and to the right.
    waxis   - vector of frequencies associated with the rows and columns
              of bic;  sampling frequency is assumed to be 1.
  """

  if w.shape != x.shape or x.shape != y.shape:
    raise ValueError('w, x and y should have identical dimentions')

  (ly, nrecs) = y.shape
  if ly == 1:
    ly = nrecs
    nrecs = 1
    w = w.reshape(1,-1)
    x = x.reshape(1,-1)
    y = y.reshape(1,-1)

  if not nfft:
    nfft = 128

  if not overlap: overlap = 50
  overlap = max(0,min(overlap,99))
  if nrecs > 1: overlap = 0
  if not nsamp: nsamp = 0
  if nrecs > 1: nsamp = ly
  if nrecs == 1 and nsamp <= 0:
    nsamp = np.fix(ly/ (8 - 7 * overlap/100))
  if nfft < nsamp:
    nfft = 2**nextpow2(nsamp)

  overlap = np.fix(overlap/100 * nsamp)
  nadvance = nsamp - overlap
  nrecs = np.fix((ly*nrecs - overlap) / nadvance)

  if not wind:
    wind = np.hanning(nsamp)

  try:
    (rw, cw) = wind.shape
  except ValueError:
    (rw,) = wind.shape
    cw = 1

  if min(rw, cw) != 1 or max(rw, cw) != nsamp:
    print "Segment size is " + str(nsamp)
    print "Wind array is " + str(rw) + " by " + str(cw)
    print "Using default Hanning window"
    wind = np.hanning(nsamp)

  wind = wind.reshape(1,-1)


  # Accumulate triple products
  bic = np.zeros([nfft, nfft])
  Pyy = np.zeros([nfft,1])
  Pww = np.zeros([nfft,1])
  Pxx = np.zeros([nfft,1])

  mask = hankel(np.arange(nfft),np.array([nfft-1]+range(nfft-1)))
  Yf12 = np.zeros([nfft,nfft])
  ind  = np.transpose(np.arange(nsamp))
  w = w.ravel(order='F')
  x = x.ravel(order='F')
  y = y.ravel(order='F')

  for k in xrange(nrecs):
    ws = w[ind]
    ws = (ws - np.mean(ws)) * wind
    Wf = np.fft.fft(ws, nfft) / nsamp
    CWf = np.conjugate(Wf)
    Pww = Pww + flat_eq(Pww, (Wf*CWf))

    xs = x[ind]
    xs = (xs - np.mean(xs)) * wind
    Xf = np.fft.fft(xs, nfft) / nsamp
    CXf = np.conjugate(Xf)
    Pxx = Pxx + flat_eq(Pxx, (Xf*CXf))

    ys = y[ind]
    ys = (ys - np.mean(ys)) * wind
    Yf = np.fft.fft(ys, nfft) / nsamp
    CYf = np.conjugate(Yf)
    Pyy = Pyy + flat_eq(Pyy, (Yf*CYf))
#.........这里部分代码省略.........

def ssa(x, r): # x is nt x 1 vector, r is embedding dimension
    #nt = len(x); mean = sum(x)/nt
    #x = x - ones(nt)*mean
    H = linalg.hankel(x, zeros(r)) #Construct Hankel matrix
    return linalg.svd(H, compute_uv=False)

#.........这里部分代码省略.........
    ...         ax.semilogy(band, np.maximum(gains, 1e-7), 'k--', linewidth=2)
    ...     if bi == 0:
    ...         ax.legend(hs, ('firls', 'remez', 'firwin2'),
    ...                   loc='lower center', frameon=False)
    ...     else:
    ...         ax.set_xlabel('Frequency (Hz)')
    ...     ax.grid(True)
    ...     ax.set(title='Band-pass %d-%d Hz' % bands[2:4], ylabel='Magnitude')
    ...
    >>> fig.tight_layout()
    >>> plt.show()

    """  # noqa
    nyq = 0.5 * _get_fs(fs, nyq)

    numtaps = int(numtaps)
    if numtaps % 2 == 0 or numtaps < 1:
        raise ValueError("numtaps must be odd and >= 1")
    M = (numtaps-1) // 2

    # normalize bands 0->1 and make it 2 columns
    nyq = float(nyq)
    if nyq <= 0:
        raise ValueError('nyq must be positive, got %s <= 0.' % nyq)
    bands = np.asarray(bands).flatten() / nyq
    if len(bands) % 2 != 0:
        raise ValueError("bands must contain frequency pairs.")
    bands.shape = (-1, 2)

    # check remaining params
    desired = np.asarray(desired).flatten()
    if bands.size != desired.size:
        raise ValueError("desired must have one entry per frequency, got %s "
                         "gains for %s frequencies."
                         % (desired.size, bands.size))
    desired.shape = (-1, 2)
    if (np.diff(bands) <= 0).any() or (np.diff(bands[:, 0]) < 0).any():
        raise ValueError("bands must be monotonically nondecreasing and have "
                         "width > 0.")
    if (bands[:-1, 1] > bands[1:, 0]).any():
        raise ValueError("bands must not overlap.")
    if (desired < 0).any():
        raise ValueError("desired must be non-negative.")
    if weight is None:
        weight = np.ones(len(desired))
    weight = np.asarray(weight).flatten()
    if len(weight) != len(desired):
        raise ValueError("weight must be the same size as the number of "
                         "band pairs (%s)." % (len(bands),))
    if (weight < 0).any():
        raise ValueError("weight must be non-negative.")

    # Set up the linear matrix equation to be solved, Qa = b

    # We can express Q(k,n) = 0.5 Q1(k,n) + 0.5 Q2(k,n)
    # where Q1(k,n)=q(k−n) and Q2(k,n)=q(k+n), i.e. a Toeplitz plus Hankel.

    # We omit the factor of 0.5 above, instead adding it during coefficient
    # calculation.

    # We also omit the 1/π from both Q and b equations, as they cancel
    # during solving.

    # We have that:
    #     q(n) = 1/π ∫W(ω)cos(nω)dω (over 0->π)
    # Using our nomalization ω=πf and with a constant weight W over each
    # interval f1->f2 we get:
    #     q(n) = W∫cos(πnf)df (0->1) = Wf sin(πnf)/πnf
    # integrated over each f1->f2 pair (i.e., value at f2 - value at f1).
    n = np.arange(numtaps)[:, np.newaxis, np.newaxis]
    q = np.dot(np.diff(np.sinc(bands * n) * bands, axis=2)[:, :, 0], weight)

    # Now we assemble our sum of Toeplitz and Hankel
    Q1 = toeplitz(q[:M+1])
    Q2 = hankel(q[:M+1], q[M:])
    Q = Q1 + Q2

    # Now for b(n) we have that:
    #     b(n) = 1/π ∫ W(ω)D(ω)cos(nω)dω (over 0->π)
    # Using our normalization ω=πf and with a constant weight W over each
    # interval and a linear term for D(ω) we get (over each f1->f2 interval):
    #     b(n) = W ∫ (mf+c)cos(πnf)df
    #          = f(mf+c)sin(πnf)/πnf + mf**2 cos(nπf)/(πnf)**2
    # integrated over each f1->f2 pair (i.e., value at f2 - value at f1).
    n = n[:M + 1]  # only need this many coefficients here
    # Choose m and c such that we are at the start and end weights
    m = (np.diff(desired, axis=1) / np.diff(bands, axis=1))
    c = desired[:, [0]] - bands[:, [0]] * m
    b = bands * (m*bands + c) * np.sinc(bands * n)
    # Use L'Hospital's rule here for cos(nπf)/(πnf)**2 @ n=0
    b[0] -= m * bands * bands / 2.
    b[1:] += m * np.cos(n[1:] * np.pi * bands) / (np.pi * n[1:]) ** 2
    b = np.dot(np.diff(b, axis=2)[:, :, 0], weight)

    # Now we can solve the equation (use pinv because Q can be rank deficient)
    a = np.dot(pinv(Q), b)

    # make coefficients symmetric (linear phase)
    coeffs = np.hstack((a[:0:-1], 2 * a[0], a[1:]))
    return coeffs

def get_hankel_matrix(y, L=None):
    N = len(y)
    if L is None:
        L = N // 2 + 1
    M = N - L + 1
    return la.hankel(y)[:L,:M]

    def update(self, x_n, d_n):

        self.n += 1
        self.update_buf(x_n, d_n)
        
        # get this randomness
        pr = np.random.uniform(size=(self.N)) 
        I = np.arange(self.N)[pr < self.q]
        
        if I.shape[0] > 0:
            
            # extract data from buffer
            L = self.buf_ind
            x = self.x[(self.length-1)+self.buf_ind-1::-1]
            d = self.d[self.buf_ind-1::-1]
            
            # compute more power of lambda if needed
            if L >= self.lmbd_pwr.shape[0]:
                self.lmbd_pwr = np.concatenate((self.lmbd_pwr, self.lmbd**np.arange(self.lmbd_pwr.shape[0], L+1)))
            
            # block-update parameters
            X = la.hankel(x[:L],r=x[L-1:])
            Lambda_diag = self.lmbd_pwr[:L]
            lmbd_L = self.lmbd_pwr[L]
            
            """ As of now, using BLAS dot is more efficient than the FFT based method toeplitz_multiplication.
            rhs = (self.lmbd**np.arange(L) * X.T).T
            self.R = (self.lmbd**L)*self.R + hankel_multiplication(x[:self.length], x[self.length-1:], rhs)
            """ # Use thus dot instead
            
            # block-update the auto-correlation matrix 
            rhs = (Lambda_diag * X.T).T
            self.R = lmbd_L*self.R + np.dot(X.T, rhs)
            
            # block-update of the cross-correlation vector here
            self.r_dx = lmbd_L*self.r_dx + np.dot(X.T, Lambda_diag*d)
            
            # update the sketches as needed
            for i in np.arange(self.N):

                if pr[i] < self.q:

                    self.last_update[i] = self.n


                    buf = np.zeros(self.length)
                    if self.sketch == 'RowSample':
                        buf = x[:self.length] 
                    elif self.sketch == 'Binary':
                        data = (np.sqrt(Lambda_diag) * X.T)
                        buf = np.dot(data, np.random.choice([-1,1], size=self.L))*np.sqrt(self.L)
                    x_up = buf.copy()

                    if self.mem_buf_size > 0:
                        for row in self.mem_buf:
                            x_up += np.random.choice([-1,1])*row*np.sqrt(lmbd_L)

                        # update the sketch buffer
                        self.mem_buf *= np.sqrt(lmbd_L)
                        self.mem_buf[1:,:] = self.mem_buf[0:-1,:]
                        self.mem_buf[0,:] = buf

                    # update the inverse correlation matrix
                    mu = 1/lmbd_L
                    pi = np.dot(self.P[i,:,:], x_up)

                    denom = (self.q*lmbd_L + np.inner(x_up, pi))
                    np.outer(pi, pi/denom, out=self.out_buf)
                    self.P[i,:,:] = self.P[i,:,:] - self.out_buf
                    self.P[i,:,:] *= mu
                    
                # do IHS
                self.w = self.w + \
                        np.dot(self.P[i,:,:], self.r_dx - np.dot(self.R, self.w))/self.n

            # Flush the buffers
            self.x[:self.length-1] = self.x[self.buf_ind:self.buf_ind+self.length-1]
            self.x[self.length-1:] = 0
            self.d[:] = 0
            self.buf_ind = 0

#how many data points do we have?
len(result) 

for i,item in enumerate(result):
    if i % 52 == 0:
        print i

# use 10% to test

# feature # of training set
fre=52

#form hankel matrix
from scipy.linalg import hankel        
X=hankel(result[0:-fre], result[-1-fre:-1])
y=result[fre:]

#use n datapoints
n=430

#split data into training and testing
Xtrain=X[:n]
ytrain=result[:n]
Xtest=X[n:]
ytest=y[n:]

#check if the split is correct
len(X)
len(y)
len(X[0])

 def construct_Hankel_matrices(self, y):
     ind0 = int(len(y)/2)
     # original and shifted hankel matrix
     H0 = la.hankel(y[0:ind0], y[ind0-1:2*ind0-1])
     H1 = la.hankel(y[1:ind0+1], y[ind0:2*ind0])
     return H0, H1