def sinc_interp1d(x, s, r):
    """Interpolates `x`, sampled at times `s`
    Output `y` is sampled at times `r`

    inspired from from Matlab:
    http://phaseportrait.blogspot.com/2008/06/sinc-interpolation-in-matlab.html

    :param ndarray x: input data time series
    :param ndarray s: input sampling time series (regular sample interval)
    :param ndarray r: output sampling time series
    :return ndarray: output data time series (regular sample interval)
    """

    # init
    s = sp.asarray(s)
    r = sp.asarray(r)
    x = sp.asarray(x)
    if x.ndim == 1:
        x = sp.atleast_2d(x)
    else:
        if x.shape[0] == len(s):
            x = x.T
        else:
            if x.shape[1] != s.shape[0]:
                raise ValueError('x and s must be same temporal extend')
    if sp.allclose(s, r):
        return x.T
    T = s[1] - s[0]

    # resample
    sincM = sp.tile(r, (len(s), 1)) - sp.tile(s[:, sp.newaxis], (1, len(r)))
    return sp.vstack([sp.dot(xx, sp.sinc(sincM / T)) for xx in x]).T

def MakePulseDataRepLPC(pulse,spec,N,rep1,numtype = sp.complex128):
    """ This will make data by assuming the data is an autoregressive process.
        Inputs
            spec - The properly weighted spectrum.
            N - The size of the ar process used to model the filter.
            pulse - The pulse shape.
            rep1 - The number of repeats of the process.
        Outputs
            outdata - A numpy Array with the shape of the """

    lp = len(pulse)
    lenspec = len(spec)
    r1 = scfft.ifft(scfft.ifftshift(spec))
    rp1 = r1[:N]
    rp2 = r1[1:N+1]
    # Use Levinson recursion to find the coefs for the data
    xr1 = sp.linalg.solve_toeplitz(rp1, rp2)
    lpc = sp.r_[sp.ones(1), -xr1]
    # The Gain  term.
    G = sp.sqrt(sp.sum(sp.conjugate(r1[:N+1])*lpc))
    Gvec = sp.r_[G, sp.zeros(N)]
    Npnt = (N+1)*3+lp
    # Create the noise vector and normalize
    xin = sp.random.randn(rep1, Npnt)+1j*sp.random.randn(rep1, Npnt)
    xinsum = sp.tile(sp.sqrt(sp.mean(xin.real**2+xin.imag**2, axis=1))[:, sp.newaxis],(1, Npnt))
    xin = xin/xinsum
    outdata = sp.signal.lfilter(Gvec, lpc, xin, axis=1)
    outpulse = sp.tile(pulse[sp.newaxis], (rep1, 1))
    outdata = outpulse*outdata[:, 2*N:2*N+lp]
    return outdata

def phenSpecificEffects(snps,pheno1,pheno2,K=None,covs=None,test='lrt'):
    """
    Univariate fixed effects interaction test for phenotype specific SNP effects

    Args:
        snps:   [N x S] SP.array of S SNPs for N individuals (test SNPs)
        pheno1: [N x 1] SP.array of 1 phenotype for N individuals
        pheno2: [N x 1] SP.array of 1 phenotype for N individuals
        K:      [N x N] SP.array of LMM-covariance/kinship koefficients (optional)
                        If not provided, then linear regression analysis is performed
        covs:   [N x D] SP.array of D covariates for N individuals
        test:    'lrt' for likelihood ratio test (default) or 'f' for F-test

    Returns:
        limix LMM object
    """
    N=snps.shape[0]
    if K is None:
        K=SP.eye(N)
    assert (pheno1.shape[1]==pheno2.shape[1]), "Only consider equal number of phenotype dimensions"
    if covs is None:
        covs = SP.ones(N,1)
    assert (pheno1.shape[1]==1 and pheno2.shape[1]==1 and pheno1.shape[0]==N and pheno2.shape[0]==N and K.shape[0]==N and K.shape[1]==N and covs.shape[0]==N), "shapes missmatch"
    Inter = SP.zeros((N*2,1))
    Inter[0:N,0]=1
    Inter0 = SP.ones((N*2,1))
    Yinter=SP.concatenate((pheno1,pheno2),0)
    Xinter = SP.tile(snps,(2,1))
    Covitner= SP.tile(covs(2,1))
    lm = simple_interaction(snps=Xinter,pheno=Yinter,covs=Covinter,Inter=Inter,Inter0=Inter0,test=test)
    return lm

def MakePulseDataRep(pulse_shape, filt_freq, delay=16,rep=1,numtype = sp.complex128):
    """ This function will create a repxLp numpy array, where rep is number of independent
        repeats and Lp is number of pulses, of noise shaped by the filter who's frequency
        response is passed as the parameter filt_freq. The pulse shape is delayed by the parameter
        delay into the data. The noise vector that will be multiplied by the filter's frequency
        response will be zero mean complex white Gaussian noise with a power of 1. The user
        then will need to multiply the filter by its size to get the desired power from using
        the function.
        Inputs:
            pulse_shape: A numpy array that holds the shape of the single pulse.
            filt_freq - a numpy array that holds the complex frequency response of the filter
            that will be used to shape the noise data.
            delay - The number of samples that the pulse will be delayed into the
            array of noise data to avoid any problems with filter overlap.
            rep - Number of indepent samples/pulses shaped by the filter.
            numtype - The type of numbers used for the output.
        Output
            data_out - A repxLp of data that has been shaped by the filter. Points along
            The first axis are independent of each other while samples along the second
            axis are colored using the filter and multiplied by the pulse shape.
    """
    npts = len(filt_freq)
    filt_tile = sp.tile(filt_freq[sp.newaxis,:],(rep,1))
    shaperep = sp.tile(pulse_shape[sp.newaxis,:],(rep,1))
    noisereal = sp.random.randn(rep,npts).astype(numtype)
    noiseimag = sp.random.randn(rep,npts).astype(numtype)
    noise_vec =(noisereal+1j*noiseimag)/sp.sqrt(2.0)
#    noise_vec = noisereal
    mult_freq = filt_tile.astype(numtype)*noise_vec
    data = scfft.ifft(mult_freq,axis=-1)
    data_out = shaperep*data[:,delay:(delay+len(pulse_shape))]
    return data_out

def trueFeatureStats(T, R, fMap, discountFactor, stateProp=1, MAT_LIMIT=1e8):
    """ Gather the statistics needed for LSTD,
    assuming infinite data (true probabilities).
    Option: if stateProp is  < 1, then only a proportion of all 
    states will be seen as starting state for transitions """
    dim = len(fMap)
    numStates = len(T)
    statMatrix = zeros((dim, dim))
    statResidual = zeros(dim)
    ss = range(numStates)
    repVersion = False
    
    if stateProp < 1:
        ss = random.sample(ss, int(numStates * stateProp))
    elif dim * numStates**2 < MAT_LIMIT:
        repVersion = True
    
    # two variants, depending on how large we can afford our matrices to become.        
    if repVersion:    
        tmp1 = tile(fMap, (numStates,1,1))
        tmp2 = transpose(tmp1, (2,1,0))
        tmp3 = tmp2 - discountFactor * tmp1            
        tmp4 = tile(T, (dim,1,1))
        tmp4 *= transpose(tmp1, (1,2,0))
        statMatrix = tensordot(tmp3, tmp4, axes=[[0,2], [1,2]]).T
        statResidual = dot(R, dot(fMap, T).T)
    else:
        for sto in ss:
            tmp = fMap - discountFactor * repmat(fMap[:, sto], numStates, 1).T
            tmp2 = fMap * repmat(T[:, sto], dim, 1)
            statMatrix += dot(tmp2, tmp.T)             
            statResidual += R[sto] * dot(fMap, T[:, sto])
    return statMatrix, statResidual

def massmatrix_rowcols(complex,k):
    """
    Compute the row and column arrays in the COO
    format of the Whitney form mass matrix
    """
    simplices = complex[-1].simplices
    num_simplices = simplices.shape[0]
    p = complex.complex_dimension()
    
    if k == p:
        #top dimension
        rows = arange(num_simplices,dtype=simplices.dtype)
        cols = arange(num_simplices,dtype=simplices.dtype)
        return rows,cols
    
    k_faces = [tuple(x) for x in combinations(range(p+1),k+1)]

    faces_per_simplex = len(k_faces)
    num_faces = num_simplices*faces_per_simplex
    faces     = empty((num_faces,k+1),dtype=simplices.dtype)
   
    for n,face in enumerate(k_faces):
        for m,i in enumerate(face):
            faces[n::faces_per_simplex,m] = simplices[:,i]

    #faces.sort() #we can't assume that the p-simplices are sorted

    indices = simplex_array_searchsorted(complex[k].simplices,faces)

    rows = tile(indices.reshape((-1,1)),(faces_per_simplex,)).flatten()
    cols = tile(indices.reshape((-1,faces_per_simplex)),(faces_per_simplex,)).flatten()

    return rows,cols

def plot_results_interval(twosample_interval_object, xlabel='Time/hr', ylabel='expression level', title="", legend=False, *args, **kwargs):
        """
        Plot results of resampling of a (subclass of) 
        :py:class:`gptwosample.twosample.interval_smooth.GPTwoSampleInterval`.
        This method will predict some data new, for plotting purpose.
        
        **Parameters:**
        
        twosample_interval_object: :py:class:`gptwosample.twosample.interval_smooth`
            The GPTwosample resample object, from which to take the results.
        """
        
        predicted_indicators = twosample_interval_object.get_predicted_indicators()
        model_dist,Xp = twosample_interval_object.get_predicted_model_distribution()
        
        IS = SP.tile(~predicted_indicators, twosample_interval_object._n_replicates_ind)
        IJ = SP.tile(predicted_indicators, twosample_interval_object._n_replicates_comm)

        # predict GPTwoSample object with indicators as interval_indices
        if(IS.any() and IJ.any()):
            twosample_interval_object._twosample_object.predict_model_likelihoods(\
                interval_indices={individual_id:IS, common_id:IJ}, messages=False)
            twosample_interval_object._twosample_object.predict_mean_variance(Xp,\
                interval_indices={individual_id:IS, common_id:IJ})
        else:
            twosample_interval_object._twosample_object.predict_model_likelihoods(messages=False)
            twosample_interval_object._twosample_object.predict_mean_variance(Xp)
        #now plot stuff
        ax1 = PL.axes([0.15, 0.1, 0.8, 0.7])

        plot_results(twosample_interval_object._twosample_object, 
                     alpha=model_dist, 
                     legend=legend,#interval_indices={individual_id:IS, common_id:IJ},
                     xlabel=xlabel,
                     ylabel=ylabel,
                     title="", *args, **kwargs)
        
        PL.suptitle(title,fontsize=20)
        
        PL.xlim([Xp.min(), Xp.max()])
        yticks = ax1.get_yticks()[0:-1]
        ax1.set_yticks(yticks)
        
        data = twosample_interval_object._twosample_object.get_data(common_id)
        Ymax = data[1].max()
        Ymin = data[1].min()
        
        DY = Ymax - Ymin
        PL.ylim([Ymin - 0.1 * DY, Ymax + 0.1 * DY])
        #2nd. plot prob. of diff
        ax2 = PL.axes([0.15, 0.8, 0.8, 0.10], sharex=ax1)
        PL.plot(Xp, model_dist, 'k-', linewidth=2)
        PL.ylabel('$P(z(t)=1)$')
#            PL.yticks([0.0,0.5,1.0])
        PL.yticks([0.5])           
        #horizontal bar
        PL.axhline(linewidth=0.5, color='#aaaaaa', y=0.5)
        PL.ylim([0, 1])
        PL.setp(ax2.get_xticklabels(), visible=False)

    def _LMLgrad_covar(self,hyperparams,debugging=False):
        """
        evaluates the gradient of the log marginal likelihood with respect to the
        hyperparameters of the covariance function
        """
        try:
            KV = self.get_covariances(hyperparams,debugging=debugging)
        except LA.LinAlgError:
            LG.error('linalg exception in _LMLgrad_covar')
            return {'covar_r':SP.zeros(len(hyperparams['covar_r'])),'covar_c':SP.zeros(len(hyperparams['covar_c'])),'covar_r':SP.zeros(len(hyperparams['covar_r']))}
        except ValueError:
            LG.error('value error in _LMLgrad_covar')
            return {'covar_r':SP.zeros(len(hyperparams['covar_r'])),'covar_c':SP.zeros(len(hyperparams['covar_c'])),'covar_r':SP.zeros(len(hyperparams['covar_r']))}
 
        RV = {}
        Si = unravel(1./KV['S'],self.n,self.t)

        if 'covar_r' in hyperparams:
            theta = SP.zeros(len(hyperparams['covar_r']))
            for i in range(len(theta)):
                Kgrad_r = self.covar_r.Kgrad_theta(hyperparams['covar_r'],i)
                d=(KV['U_r']*SP.dot(Kgrad_r,KV['U_r'])).sum(0)
                LMLgrad_det = SP.dot(d,SP.dot(Si,KV['S_c']))
                UdKU = SP.dot(KV['U_r'].T,SP.dot(Kgrad_r,KV['U_r']))
                SYUdKU = SP.dot(UdKU,(KV['Ytilde']*SP.tile(KV['S_c'][SP.newaxis,:],(self.n,1))))
                LMLgrad_quad = - (KV['Ytilde']*SYUdKU).sum()
                LMLgrad = 0.5*(LMLgrad_det + LMLgrad_quad)
                theta[i] = LMLgrad

                if debugging:
                    Kd = SP.kron(KV['K_c'], Kgrad_r)
                    _LMLgrad = 0.5 * (KV['W']*Kd).sum()
                    assert SP.allclose(LMLgrad,_LMLgrad), 'ouch, gradient is wrong for covar_r'
                    
            RV['covar_r'] = theta

        if 'covar_c' in hyperparams:
            theta = SP.zeros(len(hyperparams['covar_c']))
            for i in range(len(theta)):
                Kgrad_c = self.covar_c.Kgrad_theta(hyperparams['covar_c'],i)

                d=(KV['U_c']*SP.dot(Kgrad_c,KV['U_c'])).sum(0)
                LMLgrad_det = SP.dot(KV['S_r'],SP.dot(Si,d))

                UdKU = SP.dot(KV['U_c'].T,SP.dot(Kgrad_c,KV['U_c']))
                SYUdKU = SP.dot((KV['Ytilde']*SP.tile(KV['S_r'][:,SP.newaxis],(1,self.t))),UdKU.T)
                LMLgrad_quad = -SP.sum(KV['Ytilde']*SYUdKU)
                LMLgrad = 0.5*(LMLgrad_det + LMLgrad_quad)
                theta[i] = LMLgrad
            
                if debugging:
                    Kd = SP.kron(Kgrad_c, KV['K_r'])
                    _LMLgrad = 0.5 * (KV['W']*Kd).sum()
                    assert SP.allclose(LMLgrad,_LMLgrad), 'ouch, gradient is wrong for covar_c'
                    
                RV['covar_c'] = theta

        return RV

def regular_cube_innerproduct(rcc,k):      
    """
    For a given regular_cube_complex, compute a matrix
    representing the k-form innerproduct.

    These elements are similar to Whitney forms,
    except using standard linear (bilinear,trilinear,..)
    elements for 0-forms.
    """

    N = rcc.complex_dimension()

    #standard cube is [0,0,..,0] [0,1,...,N]   
    standard_cube  = atleast_2d(array([0]*N + range(N),dtype='i'))
    standard_k_faces = standard_cube
    for i in range(N,k,-1):        
        standard_k_faces = cube_array_boundary(standard_k_faces,i)[0]

        
    k_faces_per_cube = standard_k_faces.shape[0]


    K = zeros((k_faces_per_cube,k_faces_per_cube)) #local stiffness matrix
    h = 1
    V = h**N #cube volume
    scale = V * (1/h)**2 * (1/3.0)**(N-k)
    for i,row_i in enumerate(standard_k_faces):
        for j,row_j in enumerate(standard_k_faces):
            if all(row_i[N:] == row_j[N:]):
                differences = (row_i[:N] != row_j[:N])
                differences[row_i[N:]] = 0                
                K[i,j] = scale * (1.0/2.0)**sum(differences)
            else:
                K[i,j] = 0
        

    CA = rcc[-1].cube_array[:,:N]
    num_cubes = CA.shape[0]

    k_faces  = tile(hstack((CA,zeros((CA.shape[0],k),dtype=CA.dtype))),(1,k_faces_per_cube)).reshape((-1,N+k))
    k_faces += tile(standard_k_faces,(num_cubes,1))
    
    k_face_array = rcc[k].cube_array

    face_indices = cube_array_search(k_face_array,k_faces)

    rows = face_indices.repeat(k_faces_per_cube)
    cols = face_indices.reshape((-1,k_faces_per_cube)).repeat(k_faces_per_cube,axis=0).reshape((-1,))
    data = K.reshape((1,-1)).repeat(num_cubes,axis=0).reshape((-1,))
    
    # temporary memory cost solution - eliminate zeros from COO representation
    nz_mask = data != 0.0
    rows = rows[nz_mask]
    cols = cols[nz_mask]
    data = data[nz_mask]

    shape = (len(k_face_array),len(k_face_array))
    return coo_matrix( (data,(rows,cols)), shape).tocsr()

def plot_stoch_value():    
    #Compute Solution==========================================================
    sigma = .5
    mu = 4*sigma
    K = 7
    Gamma, eps = discretenorm.discretenorm(K,mu,sigma)
    
    N = 100
    W = sp.linspace(0,1,N)
    V = sp.zeros((N,K))
    
    u = lambda c: sp.sqrt(c)
    beta = 0.99
    
    X,Y= sp.meshgrid(W,W)
    Wdiff = Y-X
    index = Wdiff < 0
    Wdiff[index] = 0
    
    util_grid = u(Wdiff)
    
    util3 = sp.tile(util_grid[:,:,sp.newaxis],(1,1,K))
    eps_grid = eps[sp.newaxis,sp.newaxis,:]
    eps_util = eps_grid*util3
    
    Gamma_grid = Gamma[sp.newaxis,:]
    
    delta = 1
    Vprime = V
    z = 0
    while (delta > 10**-9):
        z= z+1
        V = Vprime
        gamV = Gamma_grid*V
        Expval = sp.sum(gamV,1)
        Exp_grid = sp.tile(Expval[sp.newaxis,:,sp.newaxis],(N,1,K))
        arg = eps_util+beta*Exp_grid
        arg[index] = -10^10
        Vprime = sp.amax(arg,1)
        psi_ind = sp.argmax(arg,1)
        psi = W[psi_ind]
        delta = sp.linalg.norm(Vprime - V)
    
    #============================================================    
    #Plot 3D    
    x=sp.arange(0,N)
    y=sp.arange(0,K)
    X,Y=sp.meshgrid(x,y)
    fig1 = plt.figure()
    ax1= Axes3D(fig1)
    ax1.set_xlabel(r'$W$')
    ax1.set_ylabel(r'$\varepsilon$')
    ax1.set_zlabel(r'$V$')
    ax1.plot_surface(W[X],Y,sp.transpose(Vprime), cmap=cm.coolwarm)
    plt.savefig('stoch_value.pdf')

def krondiag(v1,v2):
    """calcualte diagonal of kronecker(diag(v1),diag(v2)))
    note that this returns a non-flattened matrix
    """
    M1 = SP.tile(v1[:,SP.newaxis],[1,v2.shape[0]])
    M2 = SP.tile(v2[SP.newaxis,:],[v1.shape[0],1])
    M1 *= M2
    #RV  = (M1).ravel()
    #naive:
    #r=SP.kron(SP.diag(v1), SP.diag(v2)).diagonal()
    return M1

 def _gradQuadrForm(self, hyperparams,dK,columns =True ):
     """derivative of the quadtratic form w.r.t. kernel derivative matrix (dK)"""
     KV = self.get_covariances(hyperparams)
     Si = KV['Si']
     Ytilde = (KV['YSi'])
     if columns:
         UdKU = SP.dot(KV['Uc'].T,SP.dot(dK,KV['Uc']))
         SYUdKU = SP.dot((Ytilde*SP.tile(KV['Sr'][:,SP.newaxis],(1,Ytilde.shape[1]))),UdKU.T)
     else:
         UdKU = SP.dot(KV['Ur'].T,SP.dot(dK,KV['Ur']))
         SYUdKU = SP.dot(UdKU,(Ytilde*SP.tile(KV['Sc'][SP.newaxis,:],(Ytilde.shape[0],1))))
     return -SP.dot(Ytilde.ravel(),SYUdKU.ravel())

 def _gradQuadrFormX(self, hyperparams,dKx,columns =True ):
     """derivative of the quadtratic form with.r.t. covarianceparameters for row or column covariance"""
     KV = self.get_covariances(hyperparams)
     Ytilde = (KV['YSi'])
     if columns:
         UY=SP.dot(KV['Uc'],Ytilde.T)
         UYS = UY*SP.tile(KV['Sr'][SP.newaxis,:],(Ytilde.shape[1],1))
     else:
         UY=SP.dot(KV['Ur'],Ytilde)
         UYS = UY*SP.tile(KV['Sc'][SP.newaxis,:],(Ytilde.shape[0],1))
     UYSYU=SP.dot(UYS,UY.T)
     trUYSYUdK=(UYSYU*dKx.T).sum(0)
     return -2.0*trUYSYUdK

def MNEfit(stim,resp,order):
    # in order for dlogloss to work, we need to know -<g(yt(n),xt)>data
    # == calculate the constrained averages over the data set
    Nsamples = sp.size(stim,0)
    Ndim = sp.size(stim,1)
    psp = sp.mean(sp.mean(resp)) #spike probability (first constraint)
    avg = (1.0*stim.T*resp)/(Nsamples*1.0)
    avgs = sp.vstack((psp,avg))
    if(order > 1):
        avgsqrd = (stim.T*1.0)*(sp.array(sp.tile(resp,(1,Ndim)))*sp.array(stim))/(Nsamples*1.0)
        avgsqrd = sp.reshape(avgsqrd,(Ndim**2,1))
        avgs = sp.vstack((avgs,avgsqrd))
    
    #initialize params:
    pstart = sp.log(1/avgs[0,0] - 1)
    pstart = sp.hstack((pstart,(.001*(2*sp.random.rand(Ndim)-1))))
    if(order > 1):
        temp = .0005*(2*sp.random.rand(Ndim,Ndim)-1)
        pstart = sp.hstack((pstart,sp.reshape(temp+temp.T,(1,Ndim**2))[0]))
    
    #redefine functions with fixed vals:
    def logLoss(p):
        return LLF.log_loss(p, stim, resp, order)
    def dlogLoss(p):
        return LLF.d_log_loss(p, stim, avgs, order)
    #run the function:
    #pfinal = opt.fmin_tnc(logLoss,pstart,fprime=dlogLoss)
    # conjugate-gradient:
    pfinal = opt.fmin_cg(logLoss,pstart,fprime=dlogLoss)
    #pfinal = opt.fmin(logLoss,pstart,fprime=dlogLoss)
    return pfinal

def _non_dominated_front_arr(iterable, key=lambda x: x, allowequality=True):
    """Return a subset of items from iterable which are not dominated by any
    other item in iterable.

    Faster version, based on boolean matrix manipulations.
    """
    items = list(iterable)
    fits = map(key, items)
    l = len(items)
    x = array(fits)
    a = tile(x, (l, 1, 1))
    b = a.transpose((1, 0, 2))
    if allowequality:
        ndom = sum(a <= b, axis=2)
    else:
        ndom = sum(a < b, axis=2)
    ndom = array(ndom, dtype=bool)
    res = set()
    for ii in range(l):
        res.add(ii)
        for ij in list(res):
            if ii == ij:
                continue
            if not ndom[ij, ii]:
                res.remove(ii)
                break
            elif not ndom[ii, ij]:
                res.remove(ij)
    return set(map(lambda i: items[i], res))

def Problem1Real():
    beta = 0.9;
    T = 10;
    N = 100;
    u = lambda c: sp.sqrt(c);
    W = sp.linspace(0,1,N);
    X, Y = sp.meshgrid(W,W);
    Wdiff = Y-X
    index = Wdiff <0;
    Wdiff[index] = 0;
    util_grid = u(Wdiff);
    util_grid[index] = -10**10;
    V = sp.zeros((N,T+2));
    psi = sp.zeros((N,T+1));


    for k in xrange(T,-1,-1):
        val = util_grid + beta*sp.tile(sp.transpose(V[:,k+1]),(N,1));
        vt = sp.amax(val, axis = 1);
        psi_ind = sp.argmax(val,axis = 1)
        V[:,k]    = vt;
        psi[:,k]    = W[psi_ind];

    
    return V,psi

 def _sig_surface(self, siglevel):
     '''
     Significance surface for plotting.
     '''
     sig = wave_signif(self, siglevel, lag1(self.series))
     sig = sp.tile(sig, (len(self.series), 1)).T
     return sig

def center_on_cos(raw_quadratures, phi0=None, omega=None, snap_omega=False):
    mean = scipy.average(raw_quadratures, axis=1)
    no_angles, no_pulses = raw_quadratures.shape
    model = Model(cos_model)
    offset, amplitude, phi0, omega = guess_initial_parameters(mean, phi0, omega)
    model.set_param_hint("offset", value=offset)
    model.set_param_hint("amplitude", min=0., value=amplitude)
    model.set_param_hint("phi0", value=phi0)
    model.set_param_hint("omega", min=0., value=omega)
    model.make_params(verbose=False)
    steps = scipy.arange(no_angles)
    res = model.fit(mean, x=steps, verbose=False)
    omega_param = res.params["omega"]
    if snap_omega:
        appx_omega = float(omega_param)
        no_pi_intervals = int(round(pi/appx_omega))
        omega = pi/no_pi_intervals
        omega_param.set(omega, vary=False)
        res.fit(mean, x=steps, verbose=False)
    d_value, p_value_ks = kstest(res.residual, 'norm')
    mean_fit = res.eval(x=steps)
    offset = mean-mean_fit
    aligned_quadratures = raw_quadratures - scipy.tile(offset, (no_pulses, 1)).T
    centered_quadratures = aligned_quadratures - float(res.params["offset"])
    return (centered_quadratures,
            float(omega_param), float(res.params["phi0"]), p_value_ks)

def make_batches(data, labels=None, batch_size=100):
    if labels is not None:
        num_labels = labels.shape[1]
        cls_data = [data[find(labels[:,i] == 1)] for i in range(num_labels)]
        cls_sizes = [d.shape[0] for d in cls_data]
        cls_sels = [permutation(range(s)) for s in cls_sizes]
        n = min(cls_sizes) * len(cls_sizes)
        batch_size = min(n, batch_size)
        lpb = batch_size / num_labels
        new_dat = []
        for i in range(n/batch_size):
            for sel, cd in zip(cls_sels, cls_data):
                new_dat.append(cd[sel[i*lpb:(i+1)*lpb]])
        if sparse.issparse(data):
            data = sparse.vstack(new_dat).tocsr()
        else:
            data = np.vstack(new_dat)
        labels = np.tile(np.repeat(np.eye(num_labels),lpb,0), (n/batch_size,1))
        n = len(labels)
        perm = range(n)
    else:
        n = data.shape[0]
        perm = permutation(range(n))
    i = 0
    while i < n:
        batch = perm[i:i+batch_size]
        i += batch_size
        yield (data[batch], None) if labels is None else (data[batch], labels[batch])

def con2vert(A, b):
    """
    Convert sets of constraints to a list of vertices (of the feasible region).
    If the shape is open, con2vert returns False for the closed property.
    """
    # Python implementation of con2vert.m by Michael Kleder (July 2005),
    #  available: http://www.mathworks.com/matlabcentral/fileexchange/7894
    #  -con2vert-constraints-to-vertices
    # Author: Michael Kelder (Original)
    #         Andre Campher (Python implementation)
    c = linalg.lstsq(mat(A), mat(b))[0]
    btmp = mat(b)-mat(A)*c
    D = mat(A)/matlib.repmat(btmp, 1, A.shape[1])

    fmatv = qhull(D, "Ft") #vertices on facets

    G  = zeros((fmatv.shape[0], D.shape[1]))
    for ix in range(0, fmatv.shape[0]):
        F = D[fmatv[ix, :], :].squeeze()
        G[ix, :] = linalg.lstsq(F, ones((F.shape[0], 1)))[0].transpose()

    V = G + matlib.repmat(c.transpose(), G.shape[0], 1)
    ux = uniqm(V)

    eps = 1e-13
    Av = dot(A, ux.T)
    bv = tile(b, (1, ux.shape[0]))
    closed = sciall(Av - bv <= eps)

    return ux, closed