def VegaFilterMagnitude(filter,spectrum,redshift):
    """
    Determines the Vega magnitude (up to a constant) given an input filter,
        SED, and redshift.
    """
    from scipy.interpolate import splev,splint,splrep
    from scipy.integrate import simps
    from math import log10

    wave = spectrum[0].copy()
    data = spectrum[1].copy()

    # Redshift the spectrum and determine the valid range of wavelengths
    wave *= (1.+redshift)
    data /= (1.+redshift)
    wmin,wmax = filter[0][0],filter[0][-1]
    cond = (wave>=wmin)&(wave<=wmax)

    # Evaluate the filter at the wavelengths of the spectrum
    response = splev(wave[cond],filter)

    # Determine the total observed flux (without the bandpass correction)
    observed = splrep(wave[cond],(response*data[cond]),s=0,k=1)
    flux = splint(wmin,wmax,observed)

    # Determine the magnitude of Vega through the filter
    vwave,vdata = getSED('Vega')
    cond = (vwave>=wmin)&(vwave<=wmax)
    response = splev(vwave[cond],filter)
    vega = splrep(vwave[cond],response*vdata[cond],s=0,k=1)
    vegacorr = splint(wmin,wmax,vega)

    return -2.5*log10(flux/vegacorr)#+2.5*log10(1.+redshift)

 def setGrid(self):
     from scipy import interpolate
     x0 = numpy.logspace(-3,1,81)
     etas = numpy.linspace(0.,2.,21)
     qs = numpy.linspace(0.2,1.,17)
     grid1 = numpy.empty((x0.size,x0.size,etas.size,qs.size))
     grid2 = numpy.empty(grid1.shape)
     for i in range(qs.size):
         q = qs[i]
         q2 = q**2
         b = 1-q2
         for j in range(etas.size):
             eta = etas[j]
             g = 0.5*eta-1. # g = -1*gamma
             for k in range(x0.size):
                 x = x0[k]
                 for l in range(x0.size):
                     y = x0[l]
                     qb = ((2*x*y)/b)**g  # q_bar
                     qt = q*(x*y)**0.5/b  # q_tilde
                     sb = 0.5*(x/y - y/x) + s**2*b/(2*x*y)
                     nu1 = s**2*b/(2*x*y)
                     nu2 = nu1+ 0.5*b*(x/y + y/(x*q2))
                     nu = numpy.logspace(nu1,nu2,1001)
                     mu = nu-sb
                     t = (1+mu**2)**0.5
                     f1 = (t-mu)**0.5/t
                     f2 = (t+mu)**0.5/t
                     ng = nu**g
                     I1 = interpolate.splrep(nu,f1*ng)
                     I2 = interpolate.splrep(nu,f2*ng)
                     grid1[k,l,i,j] = qt*interpolate.splint(nu1,nu2,I1)
                     grid2[k,l,i,j] = qt*interpolate.splint(nu1,nu2,I2)
             pylab.imshow(grid1[:,:,i,j])
             pylab.show()

def rho(z0, rhopar, pop, gp):
    vec = rhopar
    rho_at_rhalf = vec[0]
    vec = vec[1:]
    # get spline representation on gp.xepol, where rhopar are defined on
    spline_n = nr(gp.xepol, vec, pop, gp)

    # and apply it to these radii, which may be anything in between
    zs =  np.log(z0/gp.Xscale[pop]) # have to integrate in d log(r)
    logrright = []; logrleft = []
    if np.rank(zs) == 0:
        if zs>0:
            logrright.append(zs)
        else:
            logrleft.append(zs)
    else:
        logrright = zs[(zs>=0.)]
        logrleft  = zs[(zs<0.)]
        logrleft  = logrleft[::-1] # inverse order

    # integrate to left and right of halflight radius
    logrhoright = []
    for i in np.arange(0, len(logrright)):
        logrhoright.append(np.log(rho_at_rhalf) + \
                           splint(0., logrright[i], spline_n))
                           # integration along dlog(r) instead of dr

    logrholeft = []
    for i in np.arange(0, len(logrleft)):
        logrholeft.append(np.log(rho_at_rhalf) + \
                          splint(0., logrleft[i], spline_n))

    tmp = np.exp(np.hstack([logrholeft[::-1], logrhoright])) # still defined on log(r)
    gh.checkpositive(tmp, 'rho()')
    return tmp

    def test_x2_surrounds_x1_sine_spline(self):
        """
        x2 range is completely above x1 range
        using a random vector to build spline
        """
        # old size
        m = 5

        # new size
        n = 6

        # bin edges
        x_old = np.linspace(0., 1., m + 1)
        x_new = np.array([-.3, -.09, 0.11, 0.14, 0.2, 0.28, 0.73])

        subbins = np.array([-.3, -.09, 0., 0.11, 0.14, 0.2, 0.28, 0.4, 0.6,
                            0.73])

        y_old = 1. + np.sin(x_old[:-1] * np.pi)

        # compute spline ----------------------------------
        x_mids = x_old[:-1] + 0.5 * np.ediff1d(x_old)
        xx = np.hstack([x_old[0], x_mids, x_old[-1]])
        yy = np.hstack([y_old[0], y_old, y_old[-1]])

        # build spline
        spl = splrep(xx, yy)

        area_old = np.array([splint(x_old[i], x_old[i + 1], spl)
                             for i in range(m)])

        # computing subbin areas
        area_subbins = np.zeros((subbins.size - 1,))
        for i in range(area_subbins.size):
            a, b = subbins[i: i + 2]
            a = max([a, x_old[0]])
            b = min([b, x_old[-1]])
            if b > a:
                area_subbins[i] = splint(a, b, spl)

        # summing subbin contributions in y_new_ref
        y_new_ref = np.zeros((x_new.size - 1,))
        y_new_ref[1] = y_old[0] * area_subbins[2] / area_old[0]
        y_new_ref[2] = y_old[0] * area_subbins[3] / area_old[0]
        y_new_ref[3] = y_old[0] * area_subbins[4] / area_old[0]
        y_new_ref[4] = y_old[1] * area_subbins[5] / area_old[1]

        y_new_ref[5]  = y_old[1] * area_subbins[6] / area_old[1]
        y_new_ref[5] += y_old[2] * area_subbins[7] / area_old[2]
        y_new_ref[5] += y_old[3] * area_subbins[8] / area_old[3]

        # call rebin function
        y_new = rebin.rebin(x_old, y_old, x_new, interp_kind=3)

        assert_allclose(y_new, y_new_ref)

def kappa(r0fine, Mrfine, nufine, sigr2nu, intbetasfine, gp):
    # for the following: enabled calculation of kappa

    # kappa_r^4
    kapr4nu = np.ones(len(r0fine)-gp.nexp)
    xint  = r0fine                  # [pc]
    yint  = gu.G1__pcMsun_1km2s_2 * Mrfine/r0fine**2  # [1/pc (km/s)^2]
    yint *= nufine                  # [Munit/pc^4 (km/s)^2]
    yint *= sigr2nu               # [Munit/pc^4 (km/s)^4
    yint *= np.exp(intbetasfine)          # [Munit/pc^4 (km/s)^4]
    gh.checkpositive(yint, 'yint in kappa_r^4')
    yscale = 10.**(1.-min(np.log10(yint[1:])))
    yint *= yscale
    # power-law extrapolation to infinity
    C = max(0., gh.quadinflog(xint[-3:], yint[-3:], r0fine[-1], gp.rinfty*r0fine[-1]))

    splpar_nu = splrep(xint, yint, k=3) # interpolation in real space
    for k in range(len(r0fine)-gp.nexp):
        # integrate from minimal radius to infinity
        kapr4nu[k] = 3.*(np.exp(-intbetasfine[k])/nufine[k]) * \
            (splint(r0fine[k], r0fine[-1], splpar_nu) + C) # [(km/s)^4]

    kapr4nu /= yscale
    gh.checkpositive(kapr4nu, 'kapr4nu in kappa_r^4')

    splpar_kap = splrep(r0fine[:-gp.nexp], np.log(kapr4nu), k=3)
    kapr4ext = np.exp(splev(r0ext, splpar_kap))
    kapr4nu = np.hstack([kapr4nu, kapr4ext])
    gh.checkpositive(kapr4nu, 'kapr4nu in extended kappa_r^4')

    dbetafinedr = splev(r0fine, splrep(r0fine, betafine), der=1)
    gh.checknan(dbetafinedr, 'dbetafinedr in kappa_r^4')

    # kappa^4_los*surfdensity
    kapl4s = np.zeros(len(r0fine)-gp.nexp)
    for k in range(len(r0fine)-gp.nexp):
        xnew = np.sqrt(r0fine[k:]**2-r0fine[k]**2)      # [pc]
        ynew = g(r0fine[k:], r0fine[k], betafine[k:], dbetafinedr[k:]) # [1]
        ynew *= nufine[k:] * kapr4nu[k:]
        C = max(0., gh.quadinflog(xnew[-3:], ynew[-3:], xnew[-1], gp.rinfty*xnew[-1]))
        splpar_nu = splrep(xnew,ynew) # not s=0.1, this sometimes gives negative entries after int
        kapl4s[k] = 2. * (splint(0., xnew[-1], splpar_nu) + C)
        #kapl4s[k] /= yscale
        # LOG('ynew = ',ynew,', kapl4s =', kapl4s[k])

    gh.checkpositive(kapl4s, 'kapl4s in kappa_r^4')

    # project kappa4_los as well
    # only use middle values to approximate, without errors in center and far
    kapl4s_out = np.exp(splev(r0, splrep(r0fine[4:-gp.nexp], kapl4s[4:], k=3))) # s=0.
    gh.checkpositive(kapl4s_out, 'kapl4s_out in kappa_r^4')
    return sigl2s_out, kapl4s_out

 def update(self, x, dt):
     self.curr_time += dt
     self.samples.append(x)
     self.time.append(self.curr_time)
     if len(self.samples) > 4:
         self.samples = self.samples[1:]
         self.time = self.time[1:]
         tck = interpolate.splrep(self.time,self.samples)
         if self.warmup:
             self.integral += interpolate.splint(self.time[-2],self.time[-1],tck)
         else:
             self.integral += interpolate.splint(self.time[0],self.time[-1],tck)
             self.warmup = True

def cumulativeCurr(file_name):
  
  input_file = open(file_name,'r')
  index = -1
  all_time = []
  all_current = []
  Q = 0.0
  input_lines = input_file.readlines()
  input_file.close()
  for input_line in input_lines:
    index += 1
    tmp = input_line.split()
    t = float(tmp[0])
    I = float(tmp[1])
    all_time.append(t)
    all_current.append(I)
  
  
  sall_time = asarray(all_time)
  sall_current = asarray(all_current)
  
  
  
  
  
  splrepint = interpolate.splrep(sall_time, sall_current, s=0)
  
  for time in all_time:
    charge = interpolate.splint(sall_time[0], time, splrepint)
    print time,charge

def calc_tissue_tac(input_tac, mtt, bv, t, lag=0):
    """
    Calculate Time/Attenuation Curve (TAC) of tissue from input TAC smoothed with spline


    Args:
      input_tac (tuple): is argument to scipy.interpolate.splint(..., tck, ...)
      mtt (float): mean transit time of tissue in seconds
      bv (float): tissue blood volume. Should be between 0 and 1
      t (np.array): time steps of output TAC
      lag (float): time which input TAC needed to get to the tissue

    Returns:
      (np.array): tissue TAC with in defined time steps
    """
    if not 0 <= bv <= 1:
        raise ValueError('bv should be in interval from 0 to 1')
    if mtt == 0:
        mtt += 0.01

    t2 = t - lag
    t2[t2 < t[0]] = t[0]
    from_t = t2 - mtt
    from_t[from_t < t2[0]] = t2[0]
    final_arr = np.array([interpolate.splint(ft, tt, input_tac) for ft, tt in zip(from_t, t2)])
    return (final_arr * bv) / mtt

def cont_kldivergence(data, prior):
    """
    Calculates the Kullback–Leibler divergence for a continuous distribution.

    data: samples from the posterior distribution
    prior: pdf of the prior distribution
    """
    x, p = np.histogram(data, bins='auto')
    x = np.array(x)/np.sum(x)
    h = len(x)
    item = 0
    zeroitem = False
    x1 = []
    p1 = []
    for i in range(h):
        if x[i] > 0 and not zeroitem:
            x1.append(x[i] / (p[1] - p[0]))
            p1.append(np.mean((p[i], p[i + 1])))
            item += 1
        if x[i] == 0 and not zeroitem:
            zeroitem = True
            initial = i
        if x[i] > 0 and zeroitem:
            zeroitem = False
            x1.append(x[i] / (p[i + 1] - p[initial]))
            p1.append(np.mean(p[initial:i + 1]))
            item += 1
    x1 = np.array(x1)
    q1 = prior(p1)
    x2 = x1*np.log(x1/q1)
    tck = interpolate.splrep(p1, x2)
    return interpolate.splint(min(p1), max(p1), tck)

def ObjectiveFunction(sigma, knots, ts, ois_rate, lib_rate):
	"""
	constructs the B-spline for the given knots, times( in years), 
	OIS rates, and LIBOR rates for the corresponding times. Calculates
	the objective function for the given parameters.
	
	
	sigma = float (lambda from equation (40) on page 23 of lecture 1 notes)
	knots = numpy array
	ts = numpy array
	ois_rate = numpy array
	lib_rate = numpy array
	"""
		
	tck1 = interp.splrep(ts, ois_rate, t = knots)
	tck2 = interp.splrep(ts, lib_rate, t = knots)
	
	approx1 = interp.splev(ts,tck1)
	approx2 = interp.splev(ts,tck2)
	
	partial_sum = ((approx1 - ois_rate)**2).sum()
	partial_sum += (sum((approx2 - lib_rate)**2)).sum()
	partial_sum *= 0.5
		
	new_tck = []
	new_tck.append(knots)
	new_tck.append((interp.splev(ts,tck1,der=2))**2 + (interp.splev(ts, tck2, der=2)**2))
	new_tck.append(3)
	
	temp = interp.splint(ts[0], ts[-1], new_tck)
		
	partial_sum += 0.5*sigma*temp
	
	return partial_sum

 def GetSourceSize(self,z):
     self.z=z
     self.Da = astCalc.da(self.z)
     self.scale = self.Da*np.pi/180./3600.
     if len(self.srcs) == 1:
         self.Re = self.Ddic['Source 1 re']*0.05
         self.Re_lower = self.Ldic['Source 1 re']*0.05
         self.Re_upper = self.Udic['Source 1 re']*0.05
         self.Re_kpc = self.Re*self.scale
         return self.Re
     elif len(self.srcs) == 2:
         print 'test this out...!'
         Xgrid = np.logspace(-4,5,1501)
         Ygrid = np.logspace(-4,5,1501)
         Res = []
         for i in range(len(self.imgs)):
             source = self.fits[i][-3]*self.srcs[0].eval(Xgrid) + self.fits[i][-2]*self.srcs[1].eval(Xgrid)
             R = Xgrid.copy()
             light = source*2.*np.pi*R
             mod = splrep(R,light,t=np.logspace(-3.8,4.8,1301))
             intlight = np.zeros(len(R))
             for i in range(len(R)):
                 intlight[i] = splint(0,R[i],mod)
             model = splrep(intlight[:-300],R[:-300])
             reff = splev(0.5*intlight[-1],model)
             Res.append(reff*0.05)
         self.Re_v,self.Re_i = Res
         return self.Re_v, self.Re_i

 def pixeval(self,x,y):
     from numpy import cosh
     from math import pi
     from itertools import product
     from scipy.interpolate import splrep, splev, splint
     cos = numpy.cos(self.pa*pi/180.)
     sin = numpy.sin(self.pa*pi/180.)
     xp = (x-self.x)*cos+(y-self.y)*sin
     yp = (y-self.y)*cos-(x-self.x)*sin
     zp = numpy.logspace(-2,3,200)
     zp = numpy.concatenate((-zp[::-1],zp))
     array = numpy.zeros(xp.shape)
     print len(xp[0]), len(xp[1])
     for ii,j in product(range(len(xp[0])),range(len(xp[1]))):
         #print ii,j
         X = xp[ii,j]
         Y = yp[ii,j]*numpy.cos(self.i) + zp*numpy.sin(self.i)
         Z = -yp[ii,j]*numpy.sin(self.i) + zp*numpy.cos(self.i)
         rho = numpy.exp(-(X**2. + Y**2.)**0.5 / self.x0) /cosh(Z/self.y0)**2.
         mod = splrep(zp,rho)
         array[ii,j] = splint(zp[0],zp[-1],mod)
     import pylab as pl
     pl.figure()
     pl.plot(Z,rho)
     return array

 def GetSourceSize(self,kpc=False):
     self.z=source_redshifts[self.name]
     self.Da = astCalc.da(self.z)
     self.scale = self.Da*1e3*np.pi/180./3600.
     if len(self.srcs) == 1 or self.name == 'J0837':
         self.Re_v = self.Ddic['Source 1 re']*0.05
         self.Re_i = self.Re_v.copy()
         self.Re_lower = self.Ldic['Source 1 re']*0.05
         self.Re_upper = self.Udic['Source 1 re']*0.05
     elif len(self.srcs) == 2 and self.name != 'J0837':
         print 'test this out...!'
         Xgrid = np.logspace(-4,5,1501)
         Res = []
         for i in range(len(self.imgs)):
             #if self.name == 'J1605':
             #    source = 
             source = self.fits[i][-3]*self.srcs[0].eval(Xgrid) + self.fits[i][-2]*self.srcs[1].eval(Xgrid)
             R = Xgrid.copy()
             light = source*2.*np.pi*R
             mod = splrep(R,light,t=np.logspace(-3.8,4.8,1301))
             intlight = np.zeros(len(R))
             for i in range(len(R)):
                 intlight[i] = splint(0,R[i],mod)
             model = splrep(intlight[:-300],R[:-300])
             
             if len(model[1][np.where(np.isnan(model[1])==True)]>0):
                 print "arrays need to be increasing monotonically! But don't worry about it"
                 model = splrep(intlight[:-450],R[:-450])
             reff = splev(0.5*intlight[-1],model)
             Res.append(reff*0.05)
         self.Re_v,self.Re_i = Res
     if kpc:
         return [self.Re_v*self.scale, self.Re_i*self.scale]
     return [self.Re_v, self.Re_i]

 def integ(x, tck, constant=0):
     x = np.atleast_1d(x)
     out = np.zeros(x.shape[0], dtype=x.dtype)
     for n in xrange(len(out)):
         out[n] = interpolate.splint(0, x[n], tck)
 #    out += constant
     return out

def test_splint():
    """
    Evaluate the definite integral of a B-spline.
    Given the knots and coefficients of a B-spline, evaluate the definite
    integral of the smoothing polynomial between two given points.
    Parameters
    ----------
    a, b : float
        The end-points of the integration interval.
    tck : tuple
        A tuple (t,c,k) containing the vector of knots, the B-spline
        coefficients, and the degree of the spline (see `splev`).
    full_output : int, optional
        Non-zero to return optional output.
    Returns
    -------
    integral : float
        The resulting integral.
    wrk : ndarray
        An array containing the integrals of the normalized B-splines
        defined on the set of knots.
    """
    x = linspace(0,10,10)
    y = sin(x)
    tck = splrep(x, y)
    y2 = splint(x[0],x[-1],tck)
    print y2

def E_W(x,y): #ekvivalentna sirina, metodom cubic spline
    EkW=0
    tck=interpolate.splrep(x,y)
    i=0
    while(i<(N_tacaka-1)):
        EkW+=interpolate.splint(x[i],x[i+1],tck)
        i+=1 
    return EkW

def E_W(x,y):
    EkW=0
    tck=interpolate.splrep(x,y)
    i=0
    while(i<(N_tacaka-1)):
        EkW+=interpolate.splint(x[i],x[i+1],tck)
        i+=1 
    return EkW

def E_W(x,y):
    EW=0
    tck=interpolate.splrep(x,y)
    i=0
    while(i<99):
        EW+=interpolate.splint(x[i],x[i+1],tck)
        i+=1 
    return EW

def integ(x,tck,constant=-1):
    import numpy as np
    x = np.atleast_1d(x)
    out = np.zeros(x.shape, dtype=x.dtype)
    for n in xrange(len(out)):
        out[n] = interpolate.splint(0,x[n],tck)
    out += constant
    return out

def N(mu, E, T, g, n=None):
  """
  Find number of electrons by integrating g(E) * f(E,T,mu) * E**n over all E

  Parameters:
    mu: chemical potential
    E: energy grid
    T: temperature
    g: DOS
    n: momentum
  """
  if n is None:
   tck = splrep(E, g*fermi(E,T,mu))
   return splint(E[0], E[-1], tck)
  else:
   tck = splrep(E, g*fermi(E,T,mu)*E**n)
   return splint(E[0], E[-1], tck)