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

    def covariance_matrix_grid(self, endog_expval, index):

        from scipy.linalg import toeplitz
        r = np.zeros(len(endog_expval))
        r[0] = 1
        r[1:self.max_lag + 1] = self.dep_params
        return toeplitz(r), True

def X1A(xtrain,ytrain,ytest,degre=10):
    """ Excercie 1 partie A"""
    
    # RYY h = rXY     
    rYY=correlate(ytrain,ytrain)
    rXY=correlate(xtrain,ytrain)
    
    # Ah=b
    A=scl.toeplitz(rYY[:degre])
    b=rXY[: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 A')
    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('X1A.pdf')
    plt.close()
    
    if inspection:
        global X1Avars
        X1Avars=inspect.currentframe().f_locals

def test_toeplitz(N=50):
    print("Testing circulant linear algebra...")
    x = np.linspace(0, 10, N)
    y = np.vstack((np.sin(x), np.cos(x), x, x**2)).T
    c_row = np.exp(-0.5 * x ** 2)
    c_row[0] += 0.1
    cnum = circulant(c_row)
    cmat = CirculantMatrix(c_row)

    # Test dot products.
    assert np.allclose(np.dot(cnum, y[:, 0]), cmat.dot(y[:, 0]))
    assert np.allclose(np.dot(cnum, y), cmat.dot(y))

    # Test solves.
    assert np.allclose(np.linalg.solve(cnum, y[:, 0]), cmat.solve(y[:, 0]))
    assert np.allclose(np.linalg.solve(cnum, y), cmat.solve(y))

    # Test eigenvalues.
    ev = np.linalg.eigvals(cnum)
    ev = ev[np.argsort(np.abs(ev))[::-1]]
    assert np.allclose(np.abs(cmat.eigvals()), np.abs(ev))

    print("Testing Toeplitz linear algebra...")
    tnum = toeplitz(c_row)
    tmat = ToeplitzMatrix(c_row)

    # Test dot products.
    assert np.allclose(np.dot(tnum, y[:, 0]), tmat.dot(y[:, 0]))
    assert np.allclose(np.dot(tnum, y), tmat.dot(y))

    # Test solves.
    assert np.allclose(np.linalg.solve(tnum, y[:, 0]),
                       tmat.solve(y[:, 0], tol=1e-12, verbose=True))
    assert np.allclose(np.linalg.solve(tnum, y),
                       tmat.solve(y, tol=1e-12, verbose=True))

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 __init__(self, n=80, sig=0.05, err=2, **kwargs):
        self.n = n  # Number of grid points
        self.sig = float(sig)  # Gaussian kernel width
        self.err = float(err) / 100  # Percent error in data

        # Setup grid
        h = 1.0 / n
        z = np.arange(h / 2, 1 - h / 2 + h, h)

        # Compute nxn matrix K = convolution with Gaussian kernel
        gaussKernel = 1 / sqrt(np.pi) / self.sig * np.exp(-(z - h / 2) ** 2 / self.sig ** 2)
        self.K = h * np.matrix(toeplitz(gaussKernel))

        # Setup true solution, blurred and noisy data
        trueimg = 0.75 * np.where((0.1 < z) & (z < 0.25), 1, 0)
        trueimg += 0.25 * np.where((0.3 < z) & (z < 0.32), 1, 0)
        trueimg += np.where((0.5 < z) & (z < 1), 1, 0) * np.sin(2 * np.pi * z) ** 4
        blurred = self.K * np.asmatrix(trueimg).T
        blurred = np.asarray(blurred.T)[0]  # np.matrix messes up your data
        noise = self.err * np.linalg.norm(blurred) * np.random.random(n) / sqrt(n)
        self.data = blurred + noise
        self.z = z
        self.trueimg = trueimg
        self.blurred = blurred
        self.setsolver()