Python linalg.hankel函数代码示例

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)