def YICgen(X,Z,alpha, delta):
	"""
	Yt = (Kt**alpha)*((exp(Zt)*Lt)**(1-alpha))
	This function generates the output levels
	given the previously defined X = sp.array([[K],[L]])
	(X is a 2xT period matrix of capital on top and
	labor on the bottom row) and Z (previously generated
	technology shocks).
	
	It = Ktp1 - (1 - delta)*Kt
	This function generates the investment levels per period
	delta = depreciation rate of capital.
	
	Ct = Yt - It
	This function defines the consumption levels as a
	difference between output and investment.
	"""
	K = X[0,:]
	L = X[1,:]
	t = sp.shape(X)[1]
	Y = sp.zeros(t)
	I = sp.zeros(t)
	C = sp.zeros(t)
	#solve for Y in each period t
	for i in range(t):
		Y[i] = (K[i]**alpha)*((sp.exp(Z[i])*L[i])**(1.-alpha))
	#solve for I in each period t
	for i in range(t-1):
		I[i] = K[i+1] - (1. - delta)*K[i]
	#solve for C in each period t
	for i in range(t-1):
		C[i] = Y[i] - I[i]
	return Y, I, C

def cov_dvrpmllbb_to_vxyz_single(d,e_d,e_vr,pmll,pmbb,cov_pmllbb,l,b):
    """
    NAME:
       cov_dvrpmllbb_to_vxyz
    PURPOSE:
       propagate distance, radial velocity, and proper motion uncertainties to
       Galactic coordinates for scalar inputs
    INPUT:
       d - distance [kpc, as/mas for plx]
       e_d - distance uncertainty [kpc, [as/mas] for plx]
       e_vr  - low velocity uncertainty [km/s]
       pmll - proper motion in l (*cos(b)) [ [as/mas]/yr ]
       pmbb - proper motion in b [ [as/mas]/yr ]
       cov_pmllbb - uncertainty covariance for proper motion
       l - Galactic longitude [rad]
       b - Galactic lattitude [rad]
    OUTPUT:
       cov(vx,vy,vz) [3,3]
    HISTORY:
       2010-04-12 - Written - Bovy (NYU)
    """
    M= _K*sc.array([[pmll,d,0.],[pmbb,0.,d]])
    cov_dpmllbb= sc.zeros((3,3))
    cov_dpmllbb[0,0]= e_d**2.
    cov_dpmllbb[1:3,1:3]= cov_pmllbb
    cov_vlvb= sc.dot(M,sc.dot(cov_dpmllbb,M.T))
    cov_vrvlvb= sc.zeros((3,3))
    cov_vrvlvb[0,0]= e_vr**2.
    cov_vrvlvb[1:3,1:3]= cov_vlvb
    R= sc.array([[m.cos(l)*m.cos(b), m.sin(l)*m.cos(b), m.sin(b)],
                 [-m.sin(l),m.cos(l),0.],
                 [-m.cos(l)*m.sin(b),-m.sin(l)*m.sin(b), m.cos(b)]])
    return sc.dot(R.T,sc.dot(cov_vrvlvb,R))

    def deactivated_test_extreme_index(self):
        """Set of parameters know to have cased issues in the past with
        numerical stability."""

        nf = 40
        nt = 150
        n = nf * nt
        dt = 0.26214
        BW = 1. / dt / 2.
        time_stream = sp.zeros((nf, nt))
        time_stream = al.make_vect(time_stream, axis_names=("freq", "time"))
        time = dt * (sp.arange(nt) + 50)
        N = dirty_map.Noise(time_stream, time)
        # Thermal.
        thermal = sp.zeros(nf, dtype=float) + 0.0002 * BW * 2.
        thermal[22] = dirty_map.T_infinity**2
        N.add_thermal(thermal)
        # Time mean and slope.
        N.deweight_time_mean()
        N.deweight_time_slope()
        # Extreem index over_f bit.
        mode = -sp.ones(nf, dtype=float) / sp.sqrt(nf - 1)
        mode[22] = 0
        # Parameters measured from one of the data sets.  Known to screw things
        # up.
        #N.add_over_f_freq_mode(8.128e-7, -4.586, 1.0, 1.422e-7, mode, True)
        N.add_over_f_freq_mode(0.001729, -0.777, 1.0, 1e-8, mode, True)
        #N.orthogonalize_modes()
        N.finalize()
        # Check if the fast inverse works.
        N_mat = N.get_mat()
        N_mat.shape = (n, n)
        N_inv = N.get_inverse()
        N_inv.shape = (n, n)

    def __init__(self, imageData=None):
        MarkerWindowInteractor.__init__(self)

        print "PlaneWidgetsXYZ.__init__()"

        self.vtksurface = None

        self.interactButtons = (1,2,3)
        self.sharedPicker = vtk.vtkCellPicker()
        #self.sharedPicker.SetTolerance(0.005)
        self.SetPicker(self.sharedPicker)
        
        self.pwX = vtk.vtkImagePlaneWidget()
        self.pwY = vtk.vtkImagePlaneWidget()
        self.pwZ = vtk.vtkImagePlaneWidget()

        
        self.textActors = {}
        self.boxes = {}

        self.set_image_data(imageData)
        self.Render()

        self.vtk_translation = zeros(3, 'd')
        self.vtk_rotation = zeros(3, 'd')

 def fgmres(self,rhs,tol=1e-6,restrt=None,maxiter=None,callback=None):
     if maxiter == None:
         maxiter = len(rhs)
     if restrt == None:
         restrt = 2*maxiter
     # implemented as in [Saad, 1993]
     # start
     x = zeros(len(rhs))
     H = zeros((restrt+1, restrt))
     V = zeros((len(rhs),restrt))
     Z = zeros((len(rhs),restrt))
     # Arnoldi process (with modified Gramm-Schmidt)
     res = 1.
     j = 0
     r = rhs - self.point.matvec(x)
     beta = norm(r)
     V[:,0]=r/beta
     while j < maxiter and res > tol:
         Z[:,j] = self.point.psolve(V[:,j])
         w = self.point.matvec(Z[:,j])
         for i in range(j+1):
             H[i,j]=dot(w,V[:,i])
             w = w - H[i,j]*V[:,i]
         H[j+1,j] = norm(w)
         V[:,j+1]=w/H[j+1,j]
         e = zeros(j+2)
         e[0]=1.
         y, res, rank, sing_val = lstsq(H[:j+2,:j+1],beta*e)
         j += 1
         print "# GMRES| iteration :", j, "res: ", res/beta
         self.resid = r_[self.resid,res/beta]
         Zy = dot(Z[:,:j],y)
     x = x + Zy
     info = 1
     return (x,info)

    def _execute(self, x, *args, **kwargs):
        """run the clustering on a set of observations"""

        # init
        self._labels = sp.zeros((len(self.crange) * self.repeats,
                                 x.shape[0]), dtype=int) - 1
        self._gof = sp.zeros(len(self.crange) * self.repeats,
                             dtype=self.dtype)
        self._ll = sp.zeros(len(self.crange) * self.repeats,
                            dtype=self.dtype)
        self._parameters = [None] * len(self.crange) * self.repeats

        # clustering
        fit_func = {
            'kmeans': self._fit_kmeans,
            'gmm': self._fit_gmm,
            #'vbgmm': self._fit_vbgmm,
            'dpgmm': self._fit_dpgmm,
            'spectral': self._fit_spectral,
            'meanshift': self._fit_mean_shift,
            'dbscan': self._fit_dbscan
        }[self.clus_type](x)

        self._winner = sp.nanargmin(self._gof)
        self.parameters = self._parameters[self._winner]
        self.labels = self._labels[self._winner]

def extract_spikes(data, epochs):
    """extract spike waveforms according to :epochs: from :data:

    :type data: ndarray
    :param data: the signal to extract from [samples, channels]
    :type epochs: ndarray
    :param epochs: epochs to cut [[start,end]], should have common length!
    :type mc: bool
    :returns: ndarray, extracted spike waveforms from :data:
    """

    # inits and checks
    if not all(map(isinstance, [data, epochs], [sp.ndarray] * 2)):
        raise TypeError('pass sp.ndarrays!')
    ns, nc = epochs.shape[0], data.shape[1]
    if epochs.shape[0] == 0:
        return sp.zeros((0, 0))
    tf = epochs[0, 1] - epochs[0, 0]

    # extract
    rval = sp.zeros((ns, tf * nc), dtype=data.dtype)
    for s in xrange(ns):
        for c in xrange(nc):
            correct_beg = min(0, epochs[s, 0])
            correct_end = max(0, epochs[s, 1] - data.shape[0])
            rval[s, c * tf - correct_beg:(c + 1) * tf - correct_end] =\
            data[epochs[s, 0] - correct_beg:epochs[s, 1] - correct_end, c]
    return rval

def makesumrule(ptype,plen,ts,lagtype='centered'):
    """ This function will return the sum rule.
        Inputs
            ptype - The type of pulse.
            plen - Length of the pulse in seconds.
            ts - Sample time in seconds.
            lagtype -  Can be centered forward or backward.
        Output
            sumrule - A 2 x nlags numpy array that holds the summation rule.
    """
    nlags = sp.round_(plen/ts)
    if ptype.lower()=='long':
        if lagtype=='forward':
            arback=-sp.arange(nlags,dtype=int)
            arforward = sp.zeros(nlags,dtype=int)
        elif lagtype=='backward':
            arback = sp.zeros(nlags,dtype=int)
            arforward=sp.arange(nlags,dtype=int)
        else:
            arback = -sp.ceil(sp.arange(0,nlags/2.0,0.5)).astype(int)
            arforward = sp.floor(sp.arange(0,nlags/2.0,0.5)).astype(int)
        sumrule = sp.array([arback,arforward])
    elif ptype.lower()=='barker':
        sumrule = sp.array([[0],[0]])
    return sumrule

def makeinputh5(Iono,basedir):
    """This will make a h5 file for the IonoContainer that can be used as starting
    points for the fitter. The ionocontainer taken will be average over the x and y dimensions
    of space to make an average value of the parameters for each altitude.
    Inputs
    Iono - An instance of the Ionocontainer class that will be averaged over so it can
    be used for fitter starting points.
    basdir - A string that holds the directory that the file will be saved to.
    """
    # Get the parameters from the original data
    Param_List = Iono.Param_List
    dataloc = Iono.Cart_Coords
    times = Iono.Time_Vector
    velocity = Iono.Velocity
    zlist,idx = sp.unique(dataloc[:,2],return_inverse=True)
    siz = list(Param_List.shape[1:])
    vsiz = list(velocity.shape[1:])

    datalocsave = sp.column_stack((sp.zeros_like(zlist),sp.zeros_like(zlist),zlist))
    outdata = sp.zeros([len(zlist)]+siz)
    outvel = sp.zeros([len(zlist)]+vsiz)
    #  Do the averaging across space
    for izn,iz in enumerate(zlist):
        arr = sp.argwhere(idx==izn)
        outdata[izn] = sp.mean(Param_List[arr],axis=0)
        outvel[izn] = sp.mean(velocity[arr],axis=0)

    Ionoout = IonoContainer(datalocsave,outdata,times,Iono.Sensor_loc,ver=0,
                            paramnames=Iono.Param_Names, species=Iono.Species,velocity=outvel)
    Ionoout.saveh5(basedir/'startdata.h5')

    def split(self, sagi, meri):
        """ utilizes geometry.grid to change the rectangle into a generalized surface,
        it is specified with a single set of basis vectors to describe the meridonial,
        normal, and sagittal planes."""
        ins = float((sagi - 1))/sagi
        inm = float((meri - 1))/meri
        stemp = self.norm.s/sagi
        mtemp = self.meri.s/meri

        z,theta = scipy.meshgrid(scipy.linspace(-self.norm.s*ins,
                                                self.norm.s*ins,
                                                sagi),
                                 scipy.linspace(-self.meri.s*inm,
                                                self.meri.s*inm,
                                                meri))

        vecin =geometry.Vecr((self.sagi.s*scipy.ones(theta.shape),
                              theta+scipy.pi/2,
                              scipy.zeros(theta.shape))) #this produces an artificial
        # meri vector, which is in the 'y_hat' direction in the space of the cylinder
        # This is a definite patch over the larger problem, where norm is not normal
        # to the cylinder surface, but is instead the axis of rotation.  This was
        # done to match the Vecr input, which works better with norm in the z direction
               
        pt1 = geometry.Point(geometry.Vecr((scipy.zeros(theta.shape),
                                            theta,
                                            z)),
                             self)

        pt1.redefine(self._origin)

        vecin = vecin.split()

        x_hat = self + pt1 #creates a vector which includes all the centers of the subsurface

        out = []
        #this for loop makes me cringe super hard
        for i in xrange(meri):
            try:
                temp = []
                for j in xrange(sagi):
                    inp = self.rot(vecin[i][j])
                    temp += [Cyl(geometry.Vecx(x_hat.x()[:,i,j]),
                                 self._origin,
                                 [2*stemp,2*mtemp],
                                 self.sagi.s,
                                 vec=[inp, self.norm.copy()],
                                 flag=self.flag)]
                out += [temp]
            except IndexError:
                inp = self.rot(vecin[i])
                out += [Cyl(geometry.Vecx(x_hat.x()[:,i]),
                            self._origin,
                            [2*stemp,2*mtemp],
                            self.norm.s,
                            vec=[inp, self.norm.copy()],
                            flag=self.flag)]
                

        return out

 def crossOver(self, parents, nbChildren):
     """ generate a number of children by doing 1-point cross-over """
     """ change as the <choice> return quite often the same p1 and even
         several time p2 was return the same than p1 """
     xdim = self.numParameters
     shuffle(parents)
     children = []
     for i in range(len(parents)/2):
         p1 = parents[i]
         p2 = parents[i+(len(parents)/2)]
         if xdim < 2:
             children.append(p1)
             children.append(p2)
         else:
             point = choice(range(xdim-1))
             point += 1
             res = zeros(xdim)
             res[:point] = p1[:point]
             res[point:] = p2[point:]
             children.append(res)
             res = zeros(xdim)
             res[:point] = p2[:point]
             res[point:] = p1[point:]
             children.append(res)
     shuffle(children)
     if len(children) > nbChildren:
         children = children[:nbChildren]  
     elif len(children) < nbChildren:
         k = True
         while k:
            children +=sample(children,len(children)) 
            if len(children) >= nbChildren:
               children = children[:nbChildren]
               k = False
     return children

 def init(self, values):
     self.values = values.copy()
     self.prev_values = values.copy()
     self.more_prev_values = values.copy()
     self.previous_gradient = zeros(values.shape)
     self.step = zeros(values.shape)
     self.previous_error = float("-inf")

 def __init__(self,linear_solver,parameters=None):
     """ 
     input: 
     =====
         linear_solver (LinearSolver) 
             contains the linear solver that will be used in 
             each Newton iteration
         parameters (dict) 
             look at the docstring of getDefaultParameters() 
             to find out which fields there are
     behaviour:
     =========
         This class implements a Newton solver that stops when the 
         maximum number of iterations has been reached, OR the 
         relative OR absolute tolerance have been reached.
     """
     Solver.Solver.__init__(self,parameters)
     if isinstance(linear_solver,LinearSolver.LinearSolver):
         self.linsolv=linear_solver
     else:
         raise TypeError, "input argument " + linear_solver \
             + " should be a linear solver"
     self.nb_newt = 0 
     self.newton_residual = zeros((0,))
     self.newton_res_norm = zeros(0,)
    # self.newton_states = zeros((  param['max_iter'],len( self.point.getCurrentGuess()) ))
     self.newton_states = zeros((  0 ))

 def GetMat(self, s, sym=False):
     """Return the element transfer matrix for the RigidMass
     element.  If sym=True, 's' must be a symbolic string and a
     matrix of strings will be returned.  Otherwise, 's' is a
     numeric value (probably complex) and the matrix returned will
     be complex."""
     if sym:
         myparams=self.symparams
     else:
         myparams=self.params
     if self.maxsize==4 and self.usez:
         rigidmat1=rigidmatz(s,myparams)
     else:
         rigidmat1=rigidmaty(s,myparams)
     if self.maxsize==4:
         return rigidmat1
     elif self.maxsize>4:
         rigidmat2=rigidmatz(s,myparams)
         zmat=scipy.zeros(scipy.shape(rigidmat2))
     if self.maxsize==8:
         bigmat1=c_[rigidmat1,zmat]
         bigmat2=c_[zmat,rigidmat2]
         temp=r_[bigmat1,bigmat2]
         return Transform8by8(temp)
     elif self.maxsize==12:
         rigidmat0=rigidmatx(s,myparams)
         row1=c_[rigidmat0,zmat,zmat]
         t1=c_[rigidmat1,zmat]
         t2=c_[zmat,rigidmat2]
         temp=r_[t1,t2]
         temp=Transform8by8(temp)
         part2=c_[scipy.zeros((8,4)),temp]
         return r_[row1, part2]

def make_line(m=1.0, b=25.0, points=100, xstep=1.0, ysigma=5.0, data=1):
    
    #generate arrays
    x = sc.zeros(points)
    y = sc.zeros(points)
    y_center = sc.zeros(points)

    #initialize random seed
    nu.random.seed(10)
    
    #fill arrays
    for i in range(points):
        x[i] = i*xstep
    y_center = m*x + b
    y = (ysigma*nu.random.randn(points) ) + y_center
    print x,y
    
    #save array to file
    if data==1:
        data_out = raw_input('data file name, ending in .txt: ')
        if data_out == '':
            data_out = 'new_file.txt'
        for j in range(points):
            f = open(data_out, "w")
            f.write(str(x[j])+','+'\t'+str(y[j]))  
            #f.write(y[j])
        f.close()
    print 'file', data_out, 'successfully written and closed'

    #plot line
    plt.scatter(x, y)
    plt.show()

    print 'ending program...'

    def _LMLgrad_lik(self,hyperparams):
        """derivative of the likelihood parameters"""

	logtheta = hyperparams['covar']
        try:   
            KV = self.get_covariances(hyperparams)
        except linalg.LinAlgError:
            LG.error("exception caught (%s)" % (str(hyperparams)))
            return 1E6
	
        #loop through all dimensions
        #logdet term:
        Kd = 2*KV['Knoise']
        dldet = 0.5*(Kd*KV['Si']).sum(axis=0)
        #quadratic term
        y_roti = KV['y_roti']
        dlquad = -0.5 * (y_roti * Kd * y_roti).sum(axis=0)
        if VERBOSE:
            dldet_  = SP.zeros([self.d])
            dlquad_ = SP.zeros([self.d])
            for d in xrange(self.d):
                _K = KV['K'] + SP.diag(KV['Knoise'][:,d])
                _Ki = SP.linalg.inv(_K)
                dldet_[d] = 0.5* SP.dot(_Ki,SP.diag(Kd[:,d])).trace()
                dlquad_[d] = -0.5*SP.dot(self.y[:,d],SP.dot(_Ki,SP.dot(SP.diag(Kd[:,d]),SP.dot(_Ki,self.y[:,d]))))

            assert (SP.absolute(dldet-dldet_)<1E-3).all(), 'outch'
            assert (SP.absolute(dlquad-dlquad_)<1E-3).all(), 'outch'


        LMLgrad = dldet + dlquad
        RV = {'lik': LMLgrad}
    
        return RV

def GetSqAmpl(filename,Nsamp=1,channel=None,V=None,O=None,r=sc.zeros((1,2)),rp=sc.zeros((1,2)),q=None):
    """
    For the transverse channel:
    Calculates and returns Sq(sample,n,r)=<q,r|q,n><q,n|Sqp|GS>.
    Coordinates are in the order Sq(sample,n,r).
    For the longitudinal channel: input r has no effect.
    Calculates and return Sq(sample,n)=|<q,n|Sqz|GS>|^2.
    """
    if type(filename)==str:
        filename=[filename]
    attrs=GetAttr(filename[0])
    if channel==None:
        channel=attrs['channel']
    L=int(attrs['L'])
    if q==None:
        q=[float(attrs['qx']/L),float(attrs['qy'])/L]
    else:
        q=[float(q[0])/L,float(q[1])/L]
    shift=None
    if 'phasex' in attrs.keys():
        shift=[attrs['phasex']/2.0,attrs['phasey']/2.0]
    else:
        shift=[attrs['phase_shift_x']/2.0,attrs['phase_shift_y']/2.0]
    phi=attrs['phi']
    neel=attrs['neel']
    if O==None or V== None:
        H,O,E,V=GetEigSys(filename,Nsamp=Nsamp,channel=channel,q=q)
    if channel=='long':
        return sf.sqwlongamp(V,O,L,L,q,shift,phi,neel)
    elif channel=='trans':
        return sf.sqwtransamp(V,O,L,L,q,shift,phi,neel,r,rp)
    pass

def EBTransitPhase(tset, kid_x):
    ebp = atpy.Table('%s/eb_pars.txt' %dir, type = 'ascii')
    
    lc = tset.tables[1]   
    time = lc.TIME
    flux = lc.PDCSAP_FLUX
    nobs = len(time)
    lg = scipy.isfinite(time)
    pylab.figure(52)
    pylab.clf()
    pylab.plot(time[lg], flux[lg])
    npl = 2
    phase = scipy.zeros((npl, nobs))
    inTr = scipy.zeros((npl, nobs), 'int')
    period = ebp.P[ebp.KID == kid_x]
    for ipl in scipy.arange(npl):
        if ipl == 0: t0 = ebp.Ep1[ebp.KID == kid_x]
        if ipl == 1: t0 = ebp.Ep2[ebp.KID == kid_x]
        if ipl == 0: dur = ebp.Dur1[ebp.KID == kid_x]
        if ipl == 1: dur = ebp.Dur2[ebp.KID == kid_x]
        dur /= period
        counter = 0
        while (time[lg] - t0).min() < 0:
            t0 -= period
            counter += 1
            if counter > 1000: break
        ph = ((time - t0) % period) / period
        ph[ph < -0.5] += 1
        ph[ph > 0.5] -= 1
        phase[ipl,:] = ph
        inTr[ipl,:] = (abs(ph) <= dur/1.5)
    return phase, inTr

 def __init__(self, npts, k=None):
     self.npts = npts
     from ..gaussian_process import GaussianProcess
     if k is None:
         from .squared_exponential import SquaredExponentialKernel
         k = SquaredExponentialKernel(fixed_params=[True, False])
     self.gp = GaussianProcess(k, X=scipy.zeros(npts), y=scipy.zeros(npts))

def Xgen(X0,Z,PP,QQ,Xbar):
    """
    This function generates a history of X given a history
    technology shocks (Z), a P matrix, a Q matrix, and an
    intial X (X0).
    Note Xt(tilde) = PXt-1(tilde) +QZt(tilde)
    Xt=Xbar*e^Xt(tilde)
    """
    num_endog=sp.shape(PP)[1]
    T=len(Z)#sp.shape(Z)[0]
    #display(T)
    X=sp.zeros((num_endog,T))
    X[:,0]=X0
    for i in range(1,T):
        Zt=Z[i]
        Xt_1=sp.zeros((num_endog,1))
        for j in range(num_endog):
            Xt_1[j,0]=X[j,i-1]
        Xt=sp.dot(PP,Xt_1)+sp.dot(QQ,Zt)
        for k in range(num_endog):
            X[k,i]=Xt[k,0]
    exponents=sp.exp(X)
    for p in range(T):
        for q in range(num_endog):
            X[q,p]=Xbar[0,q]*exponents[q,p]
    return X