def clarray(aps, lmax, zarray, zromb=3, zwidth=None):
    """Calculate an array of C_l(z, z').

    Parameters
    ----------
    aps : function
        The angular power spectrum to calculate.
    lmax : integer
        Maximum l to calculate up to.
    zarray : array_like
        Array of z's to calculate at.
    zromb : integer
        The Romberg order for integrating over frequency samples.
    zwidth : scalar, optional
        Width of frequency channel to integrate over. If None (default),
        calculate from the separation of the first two bins.

    Returns
    -------
    aps : np.ndarray[lmax+1, len(zarray), len(zarray)]
        Array of the C_l(z,z') values.
    """

    if zromb == 0:
        return aps(np.arange(lmax+1)[:, np.newaxis, np.newaxis],
                   zarray[np.newaxis, :, np.newaxis], zarray[np.newaxis, np.newaxis, :])

    else:
        zsort = np.sort(zarray)
        zhalf = np.abs(zsort[1] - zsort[0]) / 2.0 if zwidth is None else zwidth / 2.0
        zlen = zarray.size
        zint = 2**zromb + 1
        zspace = 2.0*zhalf / 2**zromb

        za = (zarray[:, np.newaxis] + np.linspace(-zhalf, zhalf, zint)[np.newaxis, :]).flatten()

        lsections = np.array_split(np.arange(lmax+1), lmax / 50)

        cla = np.zeros((lmax+1, zlen, zlen), dtype=np.float64)

        for lsec in lsections:
            clt = aps(lsec[:, np.newaxis, np.newaxis],
                      za[np.newaxis, :, np.newaxis], za[np.newaxis, np.newaxis, :])

            clt = clt.reshape(-1, zlen, zint, zlen, zint)

            clt = si.romb(clt, dx=zspace, axis=4)
            clt = si.romb(clt, dx=zspace, axis=2)

            cla[lsec] = clt / (2*zhalf)**2 # Normalise

        return cla

def ricker_int_1d(p, u):
    """Equação a diferenças do tipo Ricker com difusão espacial:
    v[t+1] (x) = a * \int K(y-x) * v[t](y) * exp(-g * v[t](y)) dy

    Com K(y-x) := exp(-(y-x)^2 / s^2)

    v representa uma população de larvas.
    """ 
    n = len(u) - 1
    uu = np.zeros(u.size)
    for i in range(1, u.size):
        uu[i] = p['a'] * romb(u * np.exp(- p['g'] * u) * p['K'][n-i:-i], p['dx'])
    uu[0] =  p['a'] * romb(u * np.exp(- p['g'] * u) * p['K'][n:], p['dx'])
    return [uu]

def calc_ip_numint(f, g, rs, method=0):
    """
    f : continuum wave function
    g : L2 function which is not normalized
    """
    fgs = np.array([f(r) * g(r) for r in rs])
    ggs = np.array([g(r) * g(r) for r in rs])
    if method == 0:
	int_fgs = integrate.simps(fgs, rs)
	int_ggs = integrate.simps(ggs, rs)
    else:
	int_fgs = integrate.romb(fgs, rs[1]-rs[0])
	int_ggs = integrate.romb(ggs, rs[1]-rs[0])
    return int_fgs / np.sqrt(int_ggs) 

def integrate_f_jnjn(f, n, mean, delta, x_max):
    """Doesn't work for n=0, because we set f=0 for the first point.
    """

    if n == 0:
        raise NotImplementedError()

    x_max = float(x_max)
    mean = float(mean)
    delta = float(delta)

    # Reduce x_max if envelope is significant.
    if delta == 0:
        envelop_width = 1000 * x_max
    else:
        envelop_width = 5 * (2 * np.pi / delta)
    x_max = min(x_max, 5 * envelop_width)

    x, ntuple, delta_tuple = sample_jnjn(n, mean, delta, x_max)

    # Envelope of a Gaussian with width of several oscillations. This
    # controls the oscillations out to high k.
    envelope = np.exp(-0.5*(x - n / mean)**2 / envelop_width**2)

    f_eval = np.empty_like(x)
    f_eval[0] = 0
    f_eval[1:] = f(x[1:])

    lower = np.s_[:ntuple[0] + 1]
    jnjn_lower = np.empty_like(x[lower])
    jnjn_lower[0] = 0
    jnjn_lower[1:] =  jnjn(n, mean, delta, x[lower][1:])
    integral = integrate.romb(f_eval[lower] * jnjn_lower, dx=delta_tuple[0])
    
    if ntuple[1]:
        upper = np.s_[ntuple[0]:]
        jnjn_upper = approx_jnjn(n, mean, delta, x[upper])
        integral += integrate.romb(f_eval[upper] * jnjn_upper * envelope[upper],
                                   dx=delta_tuple[1])

    # XXX
    #print "n points:", ntuple
    #plt.plot(x[lower], jnjn_lower * f_eval[lower])
    #plt.plot(x[upper], jnjn_upper * f_eval[upper])
    #plt.plot(x[lower], jnjn_lower * f_eval[lower] * envelope[lower])
    #plt.plot(x[upper], jnjn_upper * f_eval[upper] * envelope[upper])
    #plt.show()

    return integral

    def spectralWeight(self, limits=None):
        """Calculates the spectral weight of the model. If an energy
         window is give, a partial spectral weight is returned.

         Parameter:
         limits -- A tuple indicating beginning and end where to calculate
                   the partial spectral weight of the model.

         Returns:
            The calculated dielectric function.
         """

        _sw = 0.0

        if not limits:
            for oscillator in self.__oscillators:
                _sw += oscillator.spectralWeight

        else:
            # Using Romberg: See http://young.physics.ucsc.edu/242/romberg.pdf
            interval = limits[1]-limits[0]

            # Finding minimal k for a smaller than 0.02 integration step
            k = ceil(log(interval/0.02-1)/log(2))

            # Determining final step size
            dx = interval/(2.0**k)

            # Create a 2**k+1 equally spaced sample
            x = np.linspace(limits[0], limits[1], 2**k+1)

            _sw = romb(np.imag(self.dielectricFunction(x)), dx)

        return _sw

    def dlnsdlnM(self,M):
        """
        Uses a top-hat window function to calculate |dlnsigma/dlnM| at a given radius.
        
        Input: R: the radius of the top-hat function. 
        Output: integral: the derivatiave of log sigma with respect to log Mass.
        """
        

        R = self.Radius(M)
        # define the vector g = kR
        g = np.exp(self.k_extended)*R
            
        # Find the mass variance.
        sigma = np.sqrt(self.MassVariance(M))
            
        # Define the derivative of the window function squared.
        window = (np.sin(g)-g*np.cos(g))*(np.sin(g)*(1-3.0/(g**2))+3.0*np.cos(g)/g)
            
        # Compile into the integrand.
        integrand = (np.exp(self.power_spectrum)/np.exp(self.k_extended))*window
            
        # Integrate using romberg integration and multiply by pre-factors.
        integral =(3.0/(2.0*sigma**2*np.pi**2*R**4))*intg.romb(integrand,dx=self.steps)
        
        return integral

def integrator_linspace(y, dx=1, method="simps", axis=-1):
    """
    Integrate y using samples along the given axis with spacing dx.
    method is in ['simps', 'trapz', 'romb']. Note that romb requires len(y) = 2**k + 1.

    Example:
    >>> nx = 2**8+1
    >>> x = np.linspace(0, np.pi, nx)
    >>> dx = np.pi/(nx-1)
    >>> y = np.sin(x)
    >>> for method in ["simps", "trapz", "romb"]: print(integrator_linspace(y, dx=dx, method=method))
    2.00000000025
    1.99997490024
    2.0
    """
    from scipy.integrate import simps, trapz, romb

    if method == "simps":
        return simps(y, x=None, dx=dx, axis=axis, even='avg')
    elif method == "trapz":
        return trapz(y, x=None, dx=dx, axis=axis)
    elif method == "romb":
        return romb(y, dx=dx, axis=axis, show=False)
    else:
        raise ValueError("Wrong method: %s" % method)

def diff_off_diag(p, u, v):
    """Equação a diferenças tipo logística com difusão espacial:
    u[t+1] (x) = a1 * v[t] * (1 - g1 * v[t])
    v[t+1] (x) = a2 * \int K(y-x) * u[t](y) dy

    Com K(y-x) := exp(-(y-x)^2 / s^2)

    u representa uma população de moscas e v a quantidade de larvas.
    """ 
    n = len(u) - 1
    uu = p['a1'] * v * (1 - p['g1'] * v)
    vv = np.zeros(v.size)
    for i in range(1, v.size):
        vv[i] = p['a2'] * romb(u * p['K'][n-i:-i], p['dx'])
    vv[0] =  p['a2'] * romb(u * p['K'][n:], p['dx'])
    return [uu, vv]

    def _compute_std_romberg(self, psi, p_sep, xmin, xmax):
        '''
        Compute the variance of the distribution function psi from xmin to xmax
        along the contours p_sep using numerical integration methods.
        '''

        x_arr = np.linspace(xmin, xmax, 257)
        dx = x_arr[1] - x_arr[0]

        Q, V = 0, 0
        for x in x_arr:
            y = np.linspace(0, p_sep(x), 257)
            dy = y[1] - y[0]
            z = psi(x, y)
            Q += romb(z, dy)
            z = x**2 * psi(x, y)
            V += romb(z, dy)
        Q *= dx
        V *= dx

        return np.sqrt(V/Q)

    def test_romb_gh_3731(self):
        # Check that romb makes maximal use of data points
        x = np.arange(2 ** 4 + 1)
        y = np.cos(0.2 * x)
        val = romb(y)
        val2, err = quad(lambda x: np.cos(0.2 * x), x.min(), x.max())
        assert_allclose(val, val2, rtol=1e-8, atol=0)

        # should be equal to romb with 2**k+1 samples
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=AccuracyWarning)
            val3 = romberg(lambda x: np.cos(0.2 * x), x.min(), x.max(), divmax=4)
        assert_allclose(val, val3, rtol=1e-12, atol=0)

    def test_romb_gh_3731(self):
        # Check that romb makes maximal use of data points
        x = np.arange(2**4+1)
        y = np.cos(0.2*x)
        val = romb(y)
        val2, err = quad(lambda x: np.cos(0.2*x), x.min(), x.max())
        assert_allclose(val, val2, rtol=1e-8, atol=0)

        # should be equal to romb with 2**k+1 samples
        with suppress_warnings() as sup:
            sup.filter(AccuracyWarning, "divmax .4. exceeded")
            val3 = romberg(lambda x: np.cos(0.2*x), x.min(), x.max(), divmax=4)
        assert_allclose(val, val3, rtol=1e-12, atol=0)

def Cmp_XC_Energy(Nc,MM,LL,mumesh,Ximesh,in12,index,R0,FX=1.0,FC=1.0):
    "Computes XC energy within LDA"
    code1="""
            #line 117 "LDA.py"
            double rho=0.0;
            for (int i=0; i<an12.extent(0); i++){
                int i1 = andex(an12(i,0));
                int i2 = andex(an12(i,1));
                rho += MM(i1,ix)*LL(i1,ir)*Nc(i)*MM(i2,ix)*LL(i2,ir);
            }
            return_val = rho/(2*M_PI*R0*R0*R0/8);
    """
    exc = ExchangeCorrelation()
    Exc = zeros((len(mumesh),len(Ximesh)),dtype=float)

    an12 = array(in12)
    andex = array(index)

    for ix,eta in enumerate(mumesh):
        for ir,xi in enumerate(Ximesh):
            rho = weave.inline(code1, ['ir', 'ix', 'Nc','an12','andex','MM','LL','R0'], type_converters=weave.converters.blitz, compiler = 'gcc')
            if rho<0.0:
                print 'WARNING: rho<0 for ', eta, xi, 'rho=', rho
                #print 'Nc=', Nc
                #dsum=0.0
                #for i,(i1,i2) in enumerate(in12):
                #    print i1, i2, MM[i1,ix]*LL[i1,ir]*Nc[i]*MM[i2,ix]*LL[i2,ir], dsum
                #    dsum += MM[i1,ix]*LL[i1,ir]*Nc[i]*MM[i2,ix]*LL[i2,ir]
                rho=1e-300
                
            rs = pow(3/(4*pi*rho),1/3.)
            Exc[ix,ir] = (2*exc.Ex(rs)*rho*FX + 2*exc.Ec(rs)*rho*FC)*(xi**2-eta**2)*(2*pi*R0**3/8)
            
    Ene = zeros(len(mumesh),dtype=float)
    for ix,eta in enumerate(mumesh):
        Ene[ix] = integrate.romb(Exc[ix,:]*(Ximesh-1.),dhXi)

    return integrate.romb(Ene,mumesh[1]-mumesh[0])

def kpc2pix(H, Om, z0):
    ''' Converts physical distances in kiloparsec to Chandra pixels at redshift z0
        H = Hubble constant in km/sec/Mpc
        Om = total matter density
    '''
    N = 64 # Number of evaluation points
    c = 3e5 # speed of light in km/s
    z = linspace(0, z0, N+1)
    E = sqrt( Om*(1+z)**3 + 1-Om )
    dC = 1e3 * c/H * romb(1/E, 1.*z0/N) # Comoving distance in Mpc
    dA = dC / (1+z0) # Angular diameter distance
    rad2arcsec = 180*3600/pi
    px2arcsec = 0.492
    return 1/dA * rad2arcsec / px2arcsec #kpc2px

def cmp_MPM(Ny, maxm,maxl, ist, ak, Modd):
    Xx = linspace(-1,1,Ny)
    
    MPM =  zeros((len(ist),maxl+1),dtype=float)
    MPxM = zeros((len(ist),maxl+1),dtype=float) 
    for i,(i1,i2) in enumerate(ist):
        (R1,Ene1,m1,p1,A1) = Sol[i1]
        (R2,Ene2,m2,p2,A2) = Sol[i2]
        
        print 'i1=', i1, 'i2=', i2, 'm1=', m1, 'm2=', m2 #, MPM,MPxM
        m = abs(m1-m2)
        
        Plm_all = array([special.lpmn(maxm,maxl,x)[0] for x in Xx])
        MM = array([Mm(x,m1,p1,ak[i1]) * Mm(x,m2,p2,ak[i2]) for x in Xx])
        Plm = transpose(Plm_all[:,m])
        for il in range(m,maxl+1):
            odd = (Modd[i1]+Modd[i2]+il+m)%2
            if not odd:  # If the function is odd, we do not calculate!
                MMP = MM*Plm[il,:]
                MMPx = MMP * Xx**2
                MPM [i,il] = integrate.romb(MMP, Xx[1]-Xx[0])
                MPxM[i,il] = integrate.romb(MMPx,Xx[1]-Xx[0])
    
    return (MPM,MPxM)

def test_marginalization():
    # Integrating out one of the variables of a 2D Gaussian should
    # yield a 1D Gaussian
    mean = np.array([2.5, 3.5])
    cov = np.array([[.5, 0.2], [0.2, .6]])
    n = 2 ** 8 + 1  # Number of samples
    delta = 6 / (n - 1)  # Grid spacing

    v = np.linspace(0, 6, n)
    xv, yv = np.meshgrid(v, v)
    pos = np.empty((n, n, 2))
    pos[:, :, 0] = xv
    pos[:, :, 1] = yv
    pdf = multivariate_normal.pdf(pos, mean, cov)

    # Marginalize over x and y axis
    margin_x = romb(pdf, delta, axis=0)
    margin_y = romb(pdf, delta, axis=1)

    # Compare with standard normal distribution
    gauss_x = norm.pdf(v, loc=mean[0], scale=cov[0, 0] ** 0.5)
    gauss_y = norm.pdf(v, loc=mean[1], scale=cov[1, 1] ** 0.5)
    assert_allclose(margin_x, gauss_x, rtol=1e-2, atol=1e-2)
    assert_allclose(margin_y, gauss_y, rtol=1e-2, atol=1e-2)

def local_term(ells, chi, delta, flat=False):
    
    out = []
    chi = float(chi)
    delta = float(delta)

    for ell in ells:
        ell = int(ell)
        p_k_interp = matter_power.p_k_interp(chi)

        if flat:
            if (ell/chi) > matter_power.K_MAX:
                out.append(0.)
                continue
            # Maximum k_par.
            k_max = np.sqrt(max(matter_power.K_MAX**2 - (ell/chi)**2, 0))
            # Reduce x_max if envelope is significant.
            if delta == 0:
                envelop_width = 1000 * k_max
            else:
                envelop_width = 5 * (2 * np.pi / delta)
            k_max = min(k_max, 5 * envelop_width)

            nk = sph_bessel.pow_2_gt(k_max * delta * 5) + 1
            k = np.linspace(0, k_max, nk, endpoint=True)

            delta_k = k_max / (nk - 1)
            
            # Envelope of a Gaussian with width of several oscillations. This
            # controls the oscillations out to high k.
            envelope = np.exp(-0.5*k**2 / envelop_width**2)

            p_k_local = p_k_interp(np.sqrt((ell/chi)**2 + k**2))
            p_k_local *= envelope
            # Fourier factor.
            p_k_local *= np.cos(k * delta)

            I1 = integrate.romb(p_k_local, dx=delta_k)
            I1 /= np.pi * chi**2
        else:
            p_k = lambda k: k**2 * p_k_interp(k)
            I1 = sph_bessel.integrate_f_jnjn(p_k, ell, chi, delta,
                    matter_power.K_MAX) * (2 / np.pi)
        out.append(I1)
    return out

def normalize(E, V, interval):
    """ Use the Romberg method of integration from the SciPy module to
        normalize the wave function.
    """
    min,max,steps = interval
    k = np.ceil(np.log2(steps - 1))
    int_interval = (min,max,np.exp2(k)+1)
    xpoints,ypoints = ode.solve(f_eq = schrod_eq,
                        dep_i = (0,1),
                        interval = int_interval,
                        order = 4,
                        return_points = True,
                        V = V,
                        E = E )
    ypoints = np.power(ypoints[0,:],2) # strip out phi values and square
    dx = (max-min)/(np.exp2(k))
    A2 = integrate.romb(ypoints,dx)
    return 1/np.sqrt(A2)

def invMFP_fast(eps, xs, gamma, field):
    """
    Calculate the inverse mean free path for given tabulated
    cross sections against an isotropic photon background:
    Fully vectorized version using Romberg integration.

    The tabulated cross sections need to be of length n = 2^i + 1
    and the tabulation points log-linearly spaced
    eps   : photon energies [J] in nucleus rest frame
    xs    : cross sections [m^2]
    gamma : (array of) nucleus Lorentz factors
    field : photon background, see photonField.py

    Returns : inverse mean free path [1/Mpc]
    """
    F = integrate.cumtrapz(x=eps, y=eps*xs, initial=0) / eps
    n = field.getDensity( np.outer(1 / (2 * gamma), eps) )
    dx = np.mean(np.diff(np.log(eps)))  # average log-spacing
    return integrate.romb(n * F, dx=dx) * Mpc / gamma

def d_plus(z, omegam, omegak, omegav):
    """
    Finds the factor D+(a), from Lukic et. al. 2007, eq. 8.
    
    Uses romberg integration with a suitable step-size. 
    
    Input: z: redshift.
    
    Output: dplus: the factor.
    """

    stepsize = step_size(0.0000001, a_from_z(z))
    a_vector = np.arange(0.0000001, a_from_z(z), stepsize)
    integrand = 1.0 / (a_vector * hubble_z(z_from_a(a_vector), omegam, omegak, omegav)) ** 3

    integral = intg.romb(integrand, dx=stepsize)
    dplus = 5.0 * omegam * hubble_z(z, omegam, omegak, omegav) * integral / 2.0

    return dplus

def rate(s_kin, xs, E,field):
    """
    Calculate interaction rate against an isotropic photon background
    for given tabulated cross sections.
    1/lambda = 1/(2E) * \int n((smax-m^2)/(4E)) / (smax-m^2) F(smax-m^2) dln(smax-m^2)
    F(smax-m^2) = \int_{smin-m^2}^{smax-m^2} sigma(s) (s - m^2) d(s-m^2)
    s = 2E eps (1 - cos(theta)) + m^2
    smax = 4E*eps + m^2
    smin = m^2
    field : n(eps) photon background, see photonField.py
    s_kin     : tabulated (s - m**2) for cross sections [J^2], size n=2^i+1, log-spaced
    xs    : tabulated cross sections [m^2]
    E     : (array of) cosmic ray energies [J]
    Returns : interaction rate [1/Mpc]
    """
    F = integrate.cumtrapz(x=s_kin, y=(s_kin)*xs, initial=0)
    n = field.getDensity(np.outer(1./(4*E), s_kin))
    ds = np.mean(np.diff(np.log(s_kin)))  # value of log-spacing
    return integrate.romb(n * F / (s_kin), dx=ds) / 2 / E * Mpc