def beamfrackernel(kernelx, kernely, length, angle):
    """
    The beam fraction intercepted by a sample, used for calculating footprints.

    Parameters
    ----------
    kernelx: array-like
        x axis for the probability kernel
    kernely: array-like
        probability kernel describing the intensity distribution of the beam
    length: float
        length of the sample
    angle: float
        angle of incidence (degrees)

    Returns
    -------
    fraction: float
        The fraction of the beam intercepted by a sample.
    """
    height_of_sample = length * np.sin(np.radians(angle))
    total = integrate.simps(kernely, kernelx)
    lowlimit = np.where(-height_of_sample / 2. >= kernelx)[0][-1]
    hilimit = np.where(height_of_sample / 2. <= kernelx)[0][0]

    area = integrate.simps(kernely[lowlimit: hilimit + 1],
                           kernelx[lowlimit: hilimit + 1])
    return area / total

def exp(x,slope,xlo,xhi,efficiency=None,num_int_points=100,subranges=None):

    xnorm = np.linspace(xlo,xhi,num_int_points)
    ynorm = np.exp(-slope*xnorm)

    if efficiency!=None:
        ynorm *= efficiency(xnorm)

    normalization = integrate.simps(ynorm,x=xnorm)

    # Subranges of the normalization.
    if subranges!=None:
        normalization = 0.0
        for sr in subranges:
            xnorm = np.linspace(sr[0],sr[1],num_int_points)
            ynorm = np.exp(-slope*xnorm)

            if efficiency!=None:
                ynorm *= efficiency(xnorm)

            normalization += integrate.simps(ynorm,x=xnorm)

    y = np.exp(-slope*x)/normalization

    #'''
    if efficiency!=None:
        y *= efficiency(x)
    #'''

    return y

def print_model_results(model_name, means, sigmas, header=True):

	med_mean = np.median(means)
	n_data = float(np.size(means))
	n_params = likelihood.ndim

	model_grids = model_grids(model_name)
	likelihood = likelihood(means, sigmas, med_mean, *model_grids[0:6])
	posterior, evidence = posterior_evidence(likelihood, *model_grids[6:8])
	AIC, AICc, BIC = info_criteria(likelihood, n_data)

	marg_post = 1.0*posterior
	while marg_post.ndim != 1:
		marg_post = integrate.simps(marg_post)

	cdf = np.zeros_like(marg_post)
	for j in range(1,cdf.size):
		cdf[j] = integrate.simps(p_ld[0:j+1],ld[0:j+1])
	interp_from_cdf = interp1d(cdf,ld)
	cdf_prctls = (int_cdf(0.5), int_cdf(0.84)-int_cdf(0.5), int_cdf(0.5)-int_cdf(0.16))

	if header==True:
		print 'model (M_k) \t N_free \t P(mu0 | D, M_k) \t log evidence \t AIC \t AICc \t BIC'
	print model_name+'\t %1i \t %.2f^{+%.2f}_{-%.2f} \t %.2f \t %.2f \t %.2f \t %.2f' % \
			(n_params, interp_from_cdf(0.5), interp_from_cdf(0.84)-interp_from_cdf(0.5), 
				interp_from_cdf(0.5)-interp_from_cdf(0.16), np.log10(evidence), AIC, AICc, BIC)

	return

    def __init__(self, filepath, maxbins=1000):
        """

        :type filepath: string
        :type maxbins: int
        """
        super(FileFilter, self).__init__()
        self.path = filepath

        data = np.loadtxt(self.path)
        wf = data[:, 0]
        tp = data[:, 1]
        if len(data[:, 0]) < maxbins:  # Re-sample large filters for performance
            wfx = np.linspace(wf[0], wf[-1], maxbins)
            tpx = griddata(wf, tp, wfx)

            wf = wfx
            tp = tpx

        self.wave = wf * u.angstrom
        self.response = tp

        self.freq = (c.c / self.wave).to(u.Hz)

        nmax = np.argmax(self.response)
        halfmax_low = self.wave[:nmax][np.argmin(np.abs(self.response[nmax] - 2 * self.response[:nmax]))]
        halfmax_hi = self.wave[nmax:][np.argmin(np.abs(self.response[nmax] - 2 * self.response[nmax:]))]

        self.fwhm = halfmax_hi - halfmax_low

        self.lambda_c = (simps(self.wave * self.response, self.wave) /
                         simps(self.response, self.wave)) * u.angstrom

        self.nu_c = (simps(self.freq * self.response, self.freq) /
                     simps(self.response, self.freq)) * u.Hz

def calculateViolation(predictLine, allocateline, startTime ,endTime):
    stepSize = 1
    violateArea = 0
    violateTime = 0
    area_violations = []
    time_violations  = []
    for i in  drange(startTime + stepSize, endTime, stepSize):
        predicted_i0 = getValue(predictLine,i - stepSize)
        predicted_i1 = getValue(predictLine,i)
        ##print("Predicty=  i0 :%s i1:%s" %(predicted_i0,predicted_i1))
        allocated_i0 = getValue(allocateline, i - stepSize)
        allocated_i1 = getValue(allocateline,i)
        area_under_predicted  = simps(y = [predicted_i0, predicted_i1] , dx = stepSize)
        area_under_allocated  = simps(y = [allocated_i0, allocated_i1] , dx = stepSize)
        ##print("Areas : %s, %s"  %(area_under_predicted,area_under_allocated))
        if area_under_allocated < area_under_predicted:
            violateArea += (area_under_predicted -area_under_allocated)
            violateTime += stepSize
        area_violations.append(violateArea)
        time_violations.append(violateTime)

    area = np.array(area_violations)
    time = np.array(time_violations)
    xvalue =  np.arange(startTime+stepSize, endTime, stepSize)
    f , (plt1, plt2) = plt.subplots(1,2, sharex= True)
    plt1.plot(xvalue,area)
    plt2.plot(xvalue,time)
    #plt.show()
    #print("ViolateArea : %s ViolateTime : %s" %(violateArea, violateTime))
    return  violateArea, violateTime

def sfq0(rdfX,rdfY,ndens,Lmax=20.0,qbins=1024,damped=None):
    minq,maxq,dq=0,Lmax,Lmax/qbins
    qs=[i*dq+minq for i in range(qbins)]
    qs[0]=1E-10

    rdfY=np.array(rdfY)
    rdfX=np.array(rdfX)
    dx = rdfX[1]-rdfX[0]

    #Extend the h(r) to get a better estimate near q0
    cr,grExtx,grExty = rdfExtend(rdfX,rdfY,ndens,rmax=50.0,Niter=25,T=1000.0,rm=2.5,eps=-1,damped=0.1)

    sf1=list()
    for i,q in enumerate(qs):
        sf1 += [1 + 4*pi*ndens * integrate.simps((rdfY-1.0)*np.sin(q*rdfX)*rdfX/q ,dx=dx)]

    import pylab as pl
    
    sf2=list()
    for i,q in enumerate(qs):
        sf2 += [1 + 4*pi*ndens * integrate.simps((grExty-1.0)*np.sin(q*grExtx)*grExtx/q ,dx=dx)]
    #        R=1.0/q
    #        alpha = np.array([(1-rij/2.0/R)**2 * (1+rij/4.0/R) if rij<2*R else 0 for rij in rdfX])
    #        sf2 += [1 + 4*pi*ndens * integrate.simps((rdfY-1.0)*np.sin(q*rdfX)*rdfX/q*alpha,dx=dx)]

    f=open("/home/acadien/Dropbox/sfq0.dat","w")
    a=map(lambda x: "%f %f\n"%(x[0],x[1]),zip(qs,sf2))
    f.writelines(a)
#    exit(0)
    pl.plot(qs,sf1)
    pl.plot(qs,sf2)
    pl.show()
    exit(0)

    def fit_single_line(self, x, y, zero_lev, err_continuum, fitting_parameters, bootstrap_iterations = 1000):

        #Simple fit
        if self.fit_dict['MC_iterations'] == 1:
            fit_output = lmfit_minimize(residual_gauss, fitting_parameters, args=(x, y, zero_lev, err_continuum))
            self.fit_dict['area_intg'] = simps(y, x) - simps(zero_lev, x)
            self.fit_dict['area_intg_err'] = 0.0
             
        #Bootstrap
        else:
            mini_posterior  = Minimizer(lnprob_gaussCurve, fitting_parameters, fcn_args = ([x, y, zero_lev, err_continuum]))
            fit_output      = mini_posterior.emcee(steps=200, params = fitting_parameters)
            
            #Bootstrap for the area of the lines
            area_array = empty(bootstrap_iterations) 
            len_x_array = len(x)
            for i in range(bootstrap_iterations):
                y_new =  y + np_normal_dist(0.0, err_continuum, len_x_array)
                area_array[i] = simps(y_new, x) - simps(zero_lev, x)
            self.fit_dict['area_intg'] = mean(area_array)
            self.fit_dict['area_intg_err'] = std(area_array)           
        
        #Store the fitting parameters
        output_params = fit_output.params
        for key in self.fit_dict['parameters_list']:
            self.fit_dict[key + '_norm'] = output_params[key].value
            self.fit_dict[key + '_norm_er'] = output_params[key].stderr
            
        return

    def __init__(self,filepath,minbins=200):
        self.path = filepath
        
#        try:
        data = numpy.loadtxt(self.path)
        wf = data[:,0]
        tp = data[:,1]
        if len(data[:,0]) < minbins: #Re-sample large filters for performance
            wfx = numpy.linspace(wf[0],wf[-1],minbins)
            tpx = griddata(wf,tp,wfx)

            wf = wfx
            tp = tpx
            
        self.wave = wf * U.angstrom
        self.response = tp
        
        self.freq = (C.c/self.wave).to(U.Hz)
        
        nmax = numpy.argmax(self.response)
        halfmax_low = self.wave[:nmax][numpy.argmin(numpy.abs(self.response[nmax] - 2*self.response[:nmax]))]
        halfmax_hi = self.wave[nmax:][numpy.argmin(numpy.abs(self.response[nmax] - 2*self.response[nmax:]))]
        print self.wave[nmax],halfmax_low, halfmax_hi
        self.fwhm = halfmax_hi-halfmax_low
        
        self.lambda_c = (simps(self.wave*self.response,self.wave) / 
                         simps(self.response,self.wave))

        self.nu_c = (simps(self.freq*self.response,self.freq) / 
                     simps(self.response,self.freq))

def SusceptibilityHF(U,GF_A,X_A):
	''' susceptibility calculated from the full spectral self-energy derivative '''
	Int1_A = FD_A*sp.imag(GF_A**2*(1.0-U*X_A))
	Int2_A = FD_A*sp.imag(GF_A**2*X_A)
	I1 = simps(Int1_A,En_A)/sp.pi
	I2 = simps(Int2_A,En_A)/sp.pi
	return 2.0*I1/(1.0+U**2*I2)

def synthesize_photometric_point(filt, wl, fl):
    import photometry
    '''Takes in a spectrum in [ang], and [erg/(s cm^2 ang)], and then returns Jy'''
    if filt in {"u","g","r","i","z"}:
        wl_key = "lam"
        p_key = "air1.0"
    else:
        wl_key = "WL"
        p_key = "RES"
    p = interp1d(photometry.responses[filt][wl_key],photometry.responses[filt][p_key])
    filt_re = photometry.responses[filt][wl_key]
    p_min = min(filt_re)
    p_max = max(filt_re)
    if p_min < wl[0] or p_max > wl[-1]:
        print("ERROR, spectrum is out of range of filter: %s" % filt)
    def px(x):
        if (x < p_min) or (x > p_max):
            return 0
        else:
            return p(x)
    pxs = map(px, wl)
    f_num = wl * pxs * fl
    f_denom = wl * pxs
    #print f_num, f_denom
    num = simps(f_num, x=wl)
    denom = simps(f_denom, x=wl)
    flux = num/denom
    #This should be in spectral flux [erg/(s cm^2 A)]
    central = denom/simps(pxs, wl)
    #Returned product is in Jy  
    return photometry.Flamb_to_Jy(flux,central)

def Compute_Delta(niom, T, mu, Sig):
  # Matsubara Frequencies
  iom = pi*T*(2*arange(niom)+1)
  # Load DOS
  DOSfile = loadtxt('2D_SL_DOS')
  # 1st column as energies
  ommesh = DOSfile[:,0]
  # 2nd column as DOS
  DOS = DOSfile[:,1]
  # Normalize
  DOS = DOS / integrate.simps(DOS, ommesh)
  # Local Green function
  Gloc = zeros(niom, dtype=complex)
  for i in range(niom):
    Re = mu - ommesh - Sig[i].real
    Im = iom[i] - Sig[i].imag
    denom = 1/(Re**2 + Im**2)
    ReInt = DOS*Re*denom
    ImInt = DOS*Im*denom
    Gloc[i] = integrate.simps(ReInt, ommesh) - 1j*integrate.simps(ImInt, ommesh)

  Delta = 1j*iom+mu-Sig-1./Gloc
  with open('Delta.inp', 'w') as f:
    for i in range(niom):
      f.write('%.8f\t%.8f\t%.8f\n'%(iom[i], Delta[i].real, Delta[i].imag))

def cosine_content(pca_space, i):
    """Measure the cosine content of the PCA projection.

    The cosine content of pca projections can be used as an indicator if a
    simulation is converged. Values close to 1 are an indicator that the
    simulation isn't converged. For values below 0.7 no statement can be made.
    If you use this function please cite [BerkHess1]_.


    Parameters
    ----------
    pca_space: array, shape (number of frames, number of components)
        The PCA space to be analyzed.
    i: int
        The index of the pca_component projectection to be analyzed.

    Returns
    -------
    A float reflecting the cosine content of the ith projection in the PCA
    space. The output is bounded by 0 and 1, with 1 reflecting an agreement
    with cosine while 0 reflects complete disagreement.

    References
    ----------
    .. [BerkHess1] Berk Hess. Convergence of sampling in protein simulations.
                   Phys. Rev. E 65, 031910 (2002).
    """
    from scipy.integrate import simps
    t = np.arange(len(pca_space))
    T = len(pca_space)
    cos = np.cos(np.pi * t * (i + 1) / T)
    return ((2.0 / T) * (simps(cos*pca_space[:, i])) ** 2 /
            simps(pca_space[:, i] ** 2))

def process_fid(fid, absint=False):
    ft = cut_fft(do_fft(cut_fid(fid)))
    if absint:
        res = integrate.simps(np.abs(ft['fft']), ft.index)
    else:
        res = integrate.simps(np.real(ft['fft']), ft.index)
    return res

def test_simpson():
    ncalls = 10
    func = lambda x: np.sin(x - 0.2414)*x + 2.

    x = np.linspace(0, 10, 250001)
    y = func(x)

    t0 = time.time()
    for i in range(ncalls):
        s1 = simpson(y, dx=x[1]-x[0])
    print("cython (odd): {0} sec for {1} calls".format(time.time() - t0,ncalls))

    t0 = time.time()
    for i in range(ncalls):
        s2 = simps(y, x=x)
    print("python (odd): {0} sec for {1} calls".format(time.time() - t0,ncalls))
    np.testing.assert_allclose(s1, s2)

    # -----------------------------------------------------
    print()
    x = np.linspace(0, 10, 250000)
    y = func(x)
    t0 = time.time()
    for i in range(ncalls):
        s1 = simpson(y, dx=x[1]-x[0])
    print("cython (even): {0} sec for {1} calls".format(time.time() - t0,ncalls))

    t0 = time.time()
    for i in range(ncalls):
        s2 = simps(y, x=x)
    print("python (even): {0} sec for {1} calls".format(time.time() - t0,ncalls))

    np.testing.assert_allclose(s1, s2)

def int_func( xin, fin=None, simple=None):
    if fin is None :
        f = copy.deepcopy(xin)
        x = numpy.arange(numpy.size(f)).astype(float)
    else:
        f = copy.deepcopy(fin)
        x = copy.deepcopy(xin)
    
    n = numpy.size(f)

    g = numpy.zeros(n)

    if simple is not None :
     # Just use trapezium rule
     
        g[0] = 0.0
        for i in range (1, n) :
            g[i] = g[i-1] + 0.5*(x[i] - x[i-1])*(f[i] + f[i-1])
         
    else:
     
        n2 = numpy.int(old_div(n,2))
     
        g[0] = 0.0
        for i in range (n2, n) :
            g[i] = simps( f[0:i+1], x[0:i+1])


            
        for i in range (1, n2) :
            g[i] = g[n-1] - simps( f[i::], x[i::])
             
    return g

def _prob_I(v, x_t):
    if v == 1:
        if prob_I_cache.has_key(tuple(x_t)):
            return prob_I_cache.get(tuple(x_t))
        else:
            sigma = 4 # todo change this

            delta = 2
            span = sigma*4
            x_range = range(1, span, delta)
            k = len(x_range)
            matx = np.zeros((k, k))
            for i in x_range:
                for j in x_range:
                    i_f = x_t[0] - (span/2) + i
                    j_f = x_t[1] - (span/2) + j

                    if (np.array([i_f, j_f]) < 0).any() or (np.array([i_f, j_f]) >= 512).any():
                        matx[i/delta, j/delta] = 0
                    else:
                        matx[i/delta, j/delta] = f([i_f, j_f], x_t, sigma)

            first = scint.simps(matx, x_range, axis=0)
            _sum = scint.simps(first, x_range)
            val = 1-_sum
            prob_I_cache[tuple(x_t)] = val
            return val
    else:
        return 0

def int_func(**kwargs):
    if iter(list(kwargs.items())) == 1 :
        f = kwargs.xin
        x = numpy.arrange(numpy.size(f))*1.
    else:
        f = kwargs.fin
        x = kwargs.xin
    
    n = numpy.size(f)
  
    g = numpy.zeros(n)
    
    if kwargs.simple != None :
     # Just use trapezium rule
     
        g[0] = 0.0
        for i in range (1, n):
            g[i] = g[i-1] + 0.5*(x[i] - x[i-1])*(f[i] + f[i-1])
         
    else:
     
        n2 = numpy.int(old_div(n,2))
     
        g[0] = 0.0
        for i in range (n2, n) :
            g[i] = integrate.simps( f[0:i], x[0:i] )
      
      
        for i in range (1, n2) :
            g[i] = g[n-1] - integrate.simps( f[i:], x[i:] )
      
   
    
    return g 

def calcQ(lamin0, specin0, mstar=1.0, helium=False, f_nu=False):
    '''
    Claculate the number of lyman ionizing photons for given spectrum
    Input spectrum must be in ergs/s/A!!
    Q = int(Lnu/hnu dnu, nu_0, inf)
    '''
    lamin = np.asarray(lamin0)
    specin = np.asarray(specin0)
    c = 2.9979e18 #ang/s
    h = 6.626e-27 #erg/s
    if helium:
        lam_0 = 304.0
    else:
        lam_0 = 911.6
    if f_nu:
        nu_0 = c/lam_0
        inds, = np.where(c/lamin >= nu_0)
        hlam, hflu = c/lamin[inds], specin[inds]
        nu = hlam[::-1]
        f_nu = hflu[::-1]
        integrand = f_nu/(h*nu)
        Q = simps(integrand, x=nu)
    else:
        inds, = np.nonzero(lamin <= lam_0)
        lam = lamin[inds]
        spec = specin[inds]
        integrand = lam*spec/(h*c)
        Q = simps(integrand, x=lam)*mstar
    return Q

def KramersKronig(f, re, im, usezero=False):
    """Return real/imaginary parts retrieved by Kramers-Kronig relations.

    formulas including singularity removal according to Boukamp (1993)
    """
    from scipy.integrate import simps

    x = f * 2. * pi
    im2 = np.zeros(im.shape)
    re2 = np.zeros(im.shape)
    re3 = np.zeros(im.shape)
    drdx = np.diff(re) / np.diff(x)
    dredx = np.hstack((drdx[0], (drdx[:-1] + drdx[1:]) / 2, drdx[-1]))
    didx = np.diff(im) / np.diff(x)
    dimdx = np.hstack((didx[0], (didx[:-1] + didx[1:]) / 2, didx[-1]))
    for num, w in enumerate(x):
        x2w2 = x**2 - w**2
        x2w2[num] = 1e-12
        fun1 = (re - re[num]) / x2w2
        fun1[num] = dredx[num] / 2 / w
        im2[num] = -simps(fun1, x) * 2. * w / pi
        fun2 = (im * w / x - im[num]) / x2w2
        re2[num] = simps(fun2, x) * 2. * w / pi + re[0]
        fun3 = (im * x - im[num] * w) / x2w2
        fun3[num] = (im[num] / w + dimdx[num]) / 2
        re3[num] = simps(fun3, x) * 2. / pi + re[-1]

    if usezero:
        re3 = re2

    return re3, im2

    def Matsubara(self, iom):
        """ evaluate self-energy at matsubara axis
        """
        gm = 1/self.a2
        (x0, dh0) = swing_make_mesh(500, 1e-5*gm, 300*gm, gm)
        F0 = array([self.Fun(x) for x in x0])
        F00 = self.Fun(0.0)
        weigh0 = abs(self.expan_i[0]) / (pi*F00)

        datai=zeros(len(x0),dtype=float)
        datar=zeros(len(x0),dtype=float)
        F0i=[]
        F0r=[]
        for n in range(len(iom)):
            omn = iom[n]
            if (omn<0.3): subtract=1
            else: subtract=0
            for i in range(len(F0)):
                datai[i] = (F0[i]-F00*subtract)/(omn**2+x0[i]**2)
                datar[i] = F0[i]*x0[i]/(omn**2+x0[i]**2)
            wi = -(omn*integrate.simps(datai, x0) + F00*pi*subtract)
            wr = -integrate.simps(datar, x0)
            F0i.append(wi)
            F0r.append(wr)
            
        return (array(F0r), array(F0i), weigh0)