Python linalg.toeplitz函数代码示例

def estimate_time_constant(Y, sn, p=None, lags=5, include_noise=False, pixels=None):
    """ 
    Estimating global time constants for the dataset Y through the autocovariance function (optional).
    The function is no longer used in the standard setting of the algorithm since every trace has its own 
    time constant.
    Inputs:
    Y: np.ndarray (2D)
        input movie data with time in the last axis
    p: positive integer
        order of AR process, default: 2
    lags: positive integer
        number of lags in the past to consider for determining time constants. Default 5
    include_noise: Boolean
        Flag to include pre-estimated noise value when determining time constants. Default: False
    pixels: np.ndarray
        Restrict estimation to these pixels (e.g., remove saturated pixels). Default: All pixels
        
    Output:
    g:  np.ndarray (p x 1) 
        Discrete time constants
    """
    if p is None:
        raise Exception("You need to define p")

    if pixels is None:
        pixels = np.arange(np.size(Y) / np.shape(Y)[-1])

    from scipy.linalg import toeplitz

    npx = len(pixels)
    g = 0
    lags += p
    XC = np.zeros((npx, 2 * lags + 1))
    for j in range(npx):
        XC[j, :] = np.squeeze(axcov(np.squeeze(Y[pixels[j], :]), lags))

    gv = np.zeros(npx * lags)
    if not include_noise:
        XC = XC[:, np.arange(lags - 1, -1, -1)]
        lags -= p

    A = np.zeros((npx * lags, p))
    for i in range(npx):
        if not include_noise:
            A[i * lags + np.arange(lags), :] = toeplitz(
                np.squeeze(XC[i, np.arange(p - 1, p + lags - 1)]), np.squeeze(XC[i, np.arange(p - 1, -1, -1)])
            )
        else:
            A[i * lags + np.arange(lags), :] = toeplitz(
                np.squeeze(XC[i, lags + np.arange(lags)]), np.squeeze(XC[i, lags + np.arange(p)])
            ) - (sn[i] ** 2) * np.eye(lags, p)
            gv[i * lags + np.arange(lags)] = np.squeeze(XC[i, lags + 1 :])

    if not include_noise:
        gv = XC[:, p:].T
        gv = np.squeeze(np.reshape(gv, (np.size(gv), 1), order="F"))

    g = np.dot(np.linalg.pinv(A), gv)

    return g

def X1B(xtrain,btrain,ytest,degre=10):
    """ Excercie 1 partie B"""
    
    # (RXX+RBB) h = rXX         
    rXX=correlate(xtrain,xtrain)
    rBB=correlate(btrain,btrain)
    
    # Ah=b
    A=scl.toeplitz(rXX[:degre])+scl.toeplitz(rBB[:degre])
    b=rXX[:degre]
    h=npl.inv(A).dot(b)
    
    # Estimation de X par filtrage h de Y
    xest=scipy.signal.lfilter(h,[1],ytest)
    
    plt.figure(1)
    plt.title('X1 B')
    plt.subplot(3,1,1)
    plt.ylabel('y')
    plt.plot(ytest,'b-')        
    plt.subplot(3,1,2)
    plt.ylabel('x estime')
    plt.plot(xest,'r-')
    plt.subplot(3,1,3)
    plt.ylabel('y, x estime')
    plt.plot(ytest,'b-')
    plt.plot(xest,'r-')    
    plt.savefig('X1B.pdf')
    plt.close()
    
    if inspection:
        global X1Bvars
        X1Bvars=inspect.currentframe().f_locals         

def simulate_data_v1(nCells = 5*10**4, nPersons = 40, seed = 123456, ratio_P =  [1., 1., 0.8, 0.1]):
	"""
		Simulates the data following the instruction presented in the article
	
	"""

	if seed is not None:
		npr.seed(seed)
		
		
		
	P = [0.49, 0.3, 0.2 , 0.01 ]
	Thetas = [np.array([0.,0, 0]), np.array([0, -2, 1]), np.array([1., 2, 0]), np.array([-2,2,1.5])]
	Z_Sigma  = [np.array([[1.27, 0.25, 0],[0.25, 0.27, -0.001],[0., -0.001, 0.001]]),
				np.array([[0.06, 0.04, -0.03],[0.04, 0.05, 0],[-0.03, 0., 0.09]]),
				np.array([[0.44, 0.08, 0.08],[0.08, 0.16, 0],[0.08, 0., 0.16]]),
				0.01*np.eye(3)]
	Sigmas = [0.1*np.eye(3), 0.1*spl.toeplitz([2.,0.5,0]),0.1* spl.toeplitz([2.,-0.5,1]),
			  0.1*spl.toeplitz([1.,.3,.3]) ] 
	
	nu = 100
		
	Y, act_Class, mus,  x = simulate_data_(Thetas, Z_Sigma, Sigmas, P, nu = nu, ratio_act = ratio_P, n_cells = nCells, n_persons = nPersons,
				seed = seed)
	
	
	return Y, act_Class, mus, Thetas, Sigmas, P

def estimate_time_constant(Y, sn, p = 2, lags = 5, include_noise = False, pixels = None):
        
    if pixels is None:
        pixels = np.arange(np.size(Y)/np.shape(Y)[-1])
    
    from scipy.linalg import toeplitz    
    
    npx = len(pixels)
    g = 0
    lags += p
    XC = np.zeros((npx,2*lags+1))
    for j in range(npx):
        XC[j,:] = np.squeeze(axcov(np.squeeze(Y[pixels[j],:]),lags))
        
    gv = np.zeros(npx*lags)
    if not include_noise:
        XC = XC[:,np.arange(lags-1,-1,-1)]
        lags -= p
        
    A = np.zeros((npx*lags,p))
    for i in range(npx):
        if not include_noise:
            A[i*lags+np.arange(lags),:] = toeplitz(np.squeeze(XC[i,np.arange(p-1,p+lags-1)]),np.squeeze(XC[i,np.arange(p-1,-1,-1)])) 
        else:
            A[i*lags+np.arange(lags),:] = toeplitz(np.squeeze(XC[i,lags+np.arange(lags)]),np.squeeze(XC[i,lags+np.arange(p)])) - (sn[i]**2)*np.eye(lags,p)
            gv[i*lags+np.arange(lags)] = np.squeeze(XC[i,lags+1:])
        
    if not include_noise:
        gv = XC[:,p:].T
        gv = np.squeeze(np.reshape(gv,(np.size(gv),1),order='F'))
        
    g = np.dot(np.linalg.pinv(A),gv)
    
    return g

def test_mv(sim, d, n, mv, string):
	row = np.zeros(d)
	row2 = np.zeros(d)
	row[0] = 4.
	row2[0] = 2.
	if d > 1:
		row[1] = -1
		row2[1] = 1
	
	prior = {'mu':-10 * np.ones(d),'Sigma':spl.toeplitz(row)}
	param = {'Sigma': spl.toeplitz(row2)}
	Y = np.empty((n,d))
	L = np.linalg.cholesky(param['Sigma'])
	for i in range(n):  # @UnusedVariable
		Y[i,:] = np.dot(L,np.random.randn(d,1)).reshape((d,))
	dist = mv(prior = prior, param = param)  
	mu_est = np.zeros((sim,dist.mu_p.shape[0]))
	t0 = time.time()
	dist.set_data(Y)
	for i in range(sim):
		dist.set_parameter(param)
		dist.set_prior(prior)
		dist.set_data(Y)
		mu_est[i,:] = dist.sample()
	
	t1 = time.time()
	string += " %.4f msec/sim (sim, n ,d ) = (%d %d,%d) "%(1000*np.double(t1-t0)/sim, sim, n, d )
	print(string)

def main():

    
    print('Testing SVM class')
    height=3
    width=5
    samples=10
    #Create the D matrix
    D_w=toeplitz(np.hstack([1,np.zeros(width-2)]),np.hstack([1,-1, np.zeros(width-2)]))
    D_h=toeplitz(np.hstack([1,np.zeros(height-2)]),np.hstack([1,-1, np.zeros(height-2)]))
#    D=np.c_[D,np.zeros(width-1)]
#    D[width-2,width-1]=-1
    D2=np.kron(D_w,np.eye(height))
    D3=np.kron(np.eye(width),D_h)
    D=np.r_[D2,D3]
    
    
    w=np.random.randn(height*width)
    X=np.random.randn(samples,height*width)
    Y=np.random.randint(2,size=samples)
    Y[Y==0]=-1
    
    SHuber=SVMHuber(huber=0.5)
    print w
    print X
    print SHuber.gradf(X,Y,w)
    print SHuber.f(X,Y,w)
    
    print('Testing P class')
    KTVS=KtvSVM(D,huber=0.01,lambda1=1,lambda2=0.1,k=1)
    print KTVS.f(X,Y,w)
    print KTVS.gradf(X,Y,w)

def simulate_data(nCells = 5*10**4, nPersons = 40, seed = 123456, ratio_P =  [1., 1., 0.8, 0.1]):
	"""
		Simulates the data following the instruction presented in the article
	
	"""

	if seed != None:
		npr.seed(seed)
		
		
	nClass = 4
	dim    = 3
	P = [0.49, 0.3, 0.2 , 0.01 ]
	Thetas = [np.array([0.,0, 0]), np.array([0, -2, 1]), np.array([1., 2, 0]), np.array([-2,2,1.5])]
	Z_Sigma  = [np.array([[1.27, 0.25, 0],[0.25, 0.27, -0.001],[0., -0.001, 0.001]]),
			    np.array([[0.06, 0.04, -0.03],[0.04, 0.05, 0],[-0.03, 0., 0.09]]),
			    np.array([[0.44, 0.08, 0.08],[0.08, 0.16, 0],[0.08, 0., 0.16]]),
			    0.01*np.eye(3)]
	Sigmas = [0.1*np.eye(3), 0.1*spl.toeplitz([2.,0.5,0]),0.1* spl.toeplitz([2.,-0.5,1]),
			  0.1*spl.toeplitz([1.,.3,.3]) ] 
	

		
	act_Class = np.zeros((nPersons,4))
	for i in range(nClass):
		act_Class[:np.ceil(nPersons*ratio_P[i]),i] = 1.
	Y = []
	
	nu  = 100
	mus = []
	for i in range(nPersons):
		mix_obj = GMM.mixture(K = np.int(np.sum(act_Class[i, :])))
		theta_temp  = []
		sigma_temp  = []
		for j in range(nClass):
			if act_Class[i, j] == 1:
				theta_temp.append(Thetas[j] +  npr.multivariate_normal(np.zeros(3), Z_Sigma[j]))
				sigma_temp.append(wishart.invwishartrand(nu, (nu - dim - 1)* Sigmas[j]))
			else:
				theta_temp.append(np.ones(dim)*np.NAN)
				sigma_temp.append(np.ones((dim,dim))*np.NAN)
		theta_temp_ = [  theta_temp[aC] for aC in np.where(act_Class[i, :] == 1)[0]]
		sigma_temp_ = [  sigma_temp[aC] for aC in np.where(act_Class[i, :] == 1)[0]]

		mix_obj.mu = theta_temp_
		mus.append(theta_temp)
		mix_obj.sigma = sigma_temp_
		
		
		p_ = np.array([ (0.2*np.random.rand()+0.9) * P[aC]  for aC in np.where(act_Class[i, :] == 1)[0]]  )
		p_ /= np.sum(p_)
		mix_obj.p = p_
		mix_obj.d = dim
		Y.append(mix_obj.simulate_data(nCells))
	mus = np.array(mus)
	return Y, act_Class, mus.T, Thetas, Sigmas, P

def maestx(y, pcs, q, norder=3,samp_seg=1,overlap=0):
    """
    MAEST  MA parameter estimation via the GM-RCLS algorithm, with Tugnait's fix
        y  - time-series (vector or matrix)
        q  - MA order
        norder - cumulant-order to use  [default = 3]
        samp_seg - samples per segment for cumulant estimation
                  [default: length of y]
        overlap - percentage overlap of segments  [default = 0]
        flag - 'biased' or 'unbiased'          [default = 'biased']
        Return: estimated MA parameter vector
    """
    assert norder>=2 and norder<=4, "Cumulant order must be 2, 3, or 4!"
    nsamp = len(y)
    overlap = max(0, min(overlap,99))

    c2 = cumx(y, pcs, 2,q, samp_seg, overlap)
    c2 = np.hstack((c2, np.zeros(q)))
    cumd = cumx(y, pcs, norder,q,samp_seg,overlap,0,0)[::-1]
    cumq = cumx(y, pcs, norder,q,samp_seg,overlap,q,q)
    cumd = np.hstack((cumd, np.zeros(q)))
    cumq[:q] = np.zeros(q)

    cmat = toeplitz(cumd, np.hstack((cumd[0],np.zeros(q))))
    rmat = toeplitz(c2,   np.hstack((c2[0],np.zeros(q))))
    amat0 = np.hstack((cmat, -rmat[:,1:q+1]))
    rvec0 = c2

    cumq = np.hstack((cumq[2*q:q-1:-1], np.zeros(q)))
    cmat4 = toeplitz(cumq, np.hstack((cumq[0],np.zeros(q))))
    c3   = cumd[:2*q+1]
    amat0 = np.vstack((np.hstack((amat0, np.zeros((3*q+1,1)))), \
            np.hstack((np.hstack((np.zeros((2*q+1,q+1)), cmat4[:,1:q+1])), \
            np.reshape(-c3,(len(c3),1))))))
    rvec0 = np.hstack((rvec0, -cmat4[:,0]))

    row_sel = range(q)+range(2*q+1,3*q+1)+range(3*q+1,4*q+1)+range(4*q+2,5*q+2)
    amat0 = amat0[row_sel,:]
    rvec0 = rvec0[row_sel]

    bvec = lstsq(amat0, rvec0)[0]
    b1 = bvec[1:q+1]/bvec[0]
    b2 = bvec[q+1:2*q+1]
    if norder == 3:
        if all(b2 > 0):
            b1 = np.sign(b1) * np.sqrt(0.5*(b1**2 + b2))
        else:
            print 'MAEST: alternative solution b1 used'
    else:
        b1 = np.sign(b2)* (abs(b1) + abs(b2)**(1./3))/2
#        if all(np.sign(b2) == np.sign(b1)):
#            b1 = np.sign(b1)* (abs(b1) + abs(b2)**(1./3))/2
#        else:
#            print 'MAEST: alternative solution b1 used'
    return np.hstack(([1], b1))

def kern2mat(H, size):
    """Create a matrix corresponding to the application of the convolution kernel H.
    
    The size argument should be a tuple (output_size, input_size).
    """
    N = H.size
    half = int((N - 1) / 2)
    Nout, Mout = size
    if Nout == Mout:
        return linalg.toeplitz(np.r_[H[half:], np.zeros(Nout - half)], np.r_[H[half:], np.zeros(Mout-half)])
    else:
        return linalg.toeplitz(np.r_[H, np.zeros(Nout - N)], np.r_[H[-1], np.zeros(Mout-1)])

    def _setup_bias_correct(self, model):

        R = np.identity(model.design.shape[0]) - np.dot(model.design, model.calc_beta)
        M = np.zeros((self.p+1,)*2)
        I = np.identity(R.shape[0])

        for i in range(self.p+1):
            Di = np.dot(R, toeplitz(I[i]))
            for j in range(self.p+1):
                Dj = np.dot(R, toeplitz(I[j]))
                M[i,j] = np.diagonal((np.dot(Di, Dj))/(1.+(i>0))).sum()
                    
        self.invM = L.inv(M)
        return

def pacf(x, periodogram=True, lagmax=None):
    """Computes the partial autocorrelation function of series `x` along
    the given axis.

:Parameters:
    x : 1D array
        Time series.
    periodogram : {True, False} optional
        Whether to use a periodogram-like estimate of the ACF or not.
    lagmax : {None, int} optional
        Maximum lag. If None, the maximum lag is set to n/4+1, with n the series
        length.
    """
    acfx = acf(x, periodogram)[:,None]
    #
    if lagmax is None:
        n = len(x) // 4 + 1
    else:
        n = min(lagmax, len(x))
    #
    arkf = np.zeros((n,n),float)
    arkf[1,1] = acfx[1,0]
    for k in range(2,n):
        res = solve(toeplitz(acfx[:k]), acfx[1:k+1]).squeeze()
        arkf[k,1:k+1] = res
    return arkf.diagonal()

def arma_pacf(ar, ma, nobs=10):
    '''partial autocorrelation function of an ARMA process

    Parameters
    ----------
    ar : array_like, 1d
        coefficient for autoregressive lag polynomial, including zero lag
    ma : array_like, 1d
        coefficient for moving-average lag polynomial, including zero lag
    nobs : int
        number of terms (lags plus zero lag) to include in returned pacf

    Returns
    -------
    pacf : array
        partial autocorrelation of ARMA process given by ar, ma

    Notes
    -----
    solves yule-walker equation for each lag order up to nobs lags

    not tested/checked yet
    '''
    apacf = np.zeros(nobs)
    acov = arma_acf(ar,ma, nobs=nobs+1)

    apacf[0] = 1.
    for k in range(2,nobs+1):
        r = acov[:k];
        apacf[k-1] = linalg.solve(linalg.toeplitz(r[:-1]), r[1:])[-1]
    return apacf

def conv_mat(lin_op):
    """Returns the coefficient matrix for CONV linear op.

    Parameters
    ----------
    lin_op : LinOp
        The conv linear op.

    Returns
    -------
    list of NumPy matrix
        The matrix representing the convolution operation.
    """
    constant = const_mat(lin_op.data)
    # Cast to 1D.
    constant = intf.from_2D_to_1D(constant)

    # Create a Toeplitz matrix with constant as columns.
    rows = lin_op.size[0]
    nonzeros = lin_op.data.size[0]
    toeplitz_col = np.zeros(rows)
    toeplitz_col[0:nonzeros] = constant

    cols = lin_op.args[0].size[0]
    toeplitz_row = np.zeros(cols)
    toeplitz_row[0] = constant[0]
    coeff = sp_la.toeplitz(toeplitz_col, toeplitz_row)

    return [np.matrix(coeff)]

 def determineOptimalFilter(self, x):
     # performs the direct solution (Wiener-Hopf solution)
     # to determine optimal filter weights for equalization (optimal w/ respect to MMSE)
     # this function expects data @ 1sps, not 2 sps
     # if the data is available @ a higher samples/sym than 1,
     # it is OK to decimate the data and pass it into this block because
     # the equalization should take care of some of the timing.  Of course,
     # if the sps is 4 or greater, something smarter could probably be done
     # than just throwing data away, so think about that
     numTrainingSyms = len(self.preSyms) 
     x_n = x[0:numTrainingSyms]       # slice the appropriate preamble data
     
     # generate the input correlation matrix
     m = numTrainingSyms       
     X = numpy.fft.fft(x_n, 128)  # the FFT Size is 128 - this will change if the preamble length changes, so beware! 
     X_magSq = numpy.square(numpy.absolute(X))
     rxx = numpy.fft.ifft(X_magSq)
     rxx = rxx/m
     
     toeplitzMatCol = rxx[0:m]
     R = linalg.toeplitz(toeplitzMatCol)
     
     # generate the P vector
     xc = numpy.correlate(self.preSyms, x_n, 'full')
     P = xc[0:m]
             
     w = numpy.linalg.solve(R,P)
     w /= numpy.amax(numpy.absolute(w))      # scale the filter coefficients
     
     return w

def _least_square_evoked(epochs_data, events, tmin, sfreq):
    """Least square estimation of evoked response from epochs data.

    Parameters
    ----------
    epochs_data : array, shape (n_channels, n_times)
        The epochs data to estimate evoked.
    events : array, shape (n_events, 3)
        The events typically returned by the read_events function.
        If some events don't match the events of interest as specified
        by event_id, they will be ignored.
    tmin : float
        Start time before event.
    sfreq : float
        Sampling frequency.

    Returns
    -------
    evokeds : array, shape (n_class, n_components, n_times)
        An concatenated array of evoked data for each event type.
    toeplitz : array, shape (n_class * n_components, n_channels)
        An concatenated array of toeplitz matrix for each event type.
    """

    n_epochs, n_channels, n_times = epochs_data.shape
    tmax = tmin + n_times / float(sfreq)

    # Deal with shuffled epochs
    events = events.copy()
    events[:, 0] -= events[0, 0] + int(tmin * sfreq)

    # Contruct raw signal
    raw = _construct_signal_from_epochs(epochs_data, events, sfreq, tmin)

    # Compute the independent evoked responses per condition, while correcting
    # for event overlaps.
    n_min, n_max = int(tmin * sfreq), int(tmax * sfreq)
    window = n_max - n_min
    n_samples = raw.shape[1]
    toeplitz = list()
    classes = np.unique(events[:, 2])
    for ii, this_class in enumerate(classes):
        # select events by type
        sel = events[:, 2] == this_class

        # build toeplitz matrix
        trig = np.zeros((n_samples, 1))
        ix_trig = (events[sel, 0]) + n_min
        trig[ix_trig] = 1
        toeplitz.append(linalg.toeplitz(trig[0:window], trig))

    # Concatenate toeplitz
    toeplitz = np.array(toeplitz)
    X = np.concatenate(toeplitz)

    # least square estimation
    predictor = np.dot(linalg.pinv(np.dot(X, X.T)), X)
    evokeds = np.dot(predictor, raw.T)
    evokeds = np.transpose(np.vsplit(evokeds, len(classes)), (0, 2, 1))
    return evokeds, toeplitz