def _initialize_defaults(self):
        if self._w0 != -1.0 or self._wa != 0.0:
            a_array = numpy.logspace(
                -16, 0,
                defaults.default_precision["cosmo_npoints"])
            self._de_presure_array = self._de_presure(1/a_array - 1.0)
            self._de_presure_spline = InterpolatedUnivariateSpline(
                numpy.log(a_array), self._de_presure_array)

        self._chi = integrate.romberg(
            self.E, 0.0, self._redshift, vec_func=True,
            tol=defaults.default_precision["cosmo_precision"])

        self.growth_norm = integrate.romberg(
            self._growth_integrand, 1e-16, 1.0, vec_func=True,
            tol=defaults.default_precision["cosmo_precision"])
        self.growth_norm *= 2.5*self._omega_m0*numpy.sqrt(self.E0(0.0))

        a = 1.0/(1.0 + self._redshift)
        growth = integrate.romberg(
            self._growth_integrand, 1e-16, a, vec_func=True,
            tol=defaults.default_precision["cosmo_precision"])
        growth *= 2.5*self._omega_m0*numpy.sqrt(self.E0(self._redshift))
        self._growth = growth/self.growth_norm

        self._sigma_norm = 1.0
        self._sigma_norm = self._sigma_8*self._growth/self.sigma_r(8.0)

    def raw_window_function(self, chi):
        a = 1.0/(1.0 + self.cosmo.redshift(chi))

        try:
            g_chi = numpy.empty(len(chi))
            for idx,value in enumerate(chi):
                chi_bound = value
                if chi_bound < self._g_chi_min: chi_bound = self._g_chi_min

                g_chi[idx] = integrate.romberg(
                    self._lensing_integrand, chi_bound,
                    self.chi_max, args=(value,), vec_func=True,
                    tol=defaults.default_precision["window_precision"])
        except TypeError:
            chi_bound = chi
            if chi_bound < self._g_chi_min: chi_bound = self._g_chi_min

            g_chi = integrate.romberg(
                self._lensing_integrand, chi_bound,
                self.chi_max, args=(chi,), vec_func=True,
                tol=defaults.default_precision["window_precision"])

        g_chi *= self.cosmo.H0*self.cosmo.H0*chi

        return 3.0/2.0*self.cosmo._omega_m0*g_chi/a

    def _normalize(self):
        self.f_norm = 1.0
        norm = integrate.romberg(
            self.f_nu, self.nu_min, self.nu_max, vec_func=True,
            tol=defaults.default_precision["global_precision"],
            rtol=defaults.default_precision["mass_precision"],
            divmax=defaults.default_precision["divmax"])
        self.f_norm = 1.0/norm

        self.bias_norm = 1.0
        norm = integrate.romberg(
            lambda x: self.f_nu(x)*self.bias_nu(x),
            self.nu_min, self.nu_max, vec_func=True,
            tol=defaults.default_precision["global_precision"],
            rtol=defaults.default_precision["mass_precision"],
            divmax=defaults.default_precision["divmax"])
        self.bias_norm = 1.0/norm
        
        self.bias_2_norm = 0.0
        norm = integrate.romberg(
            lambda x: self.f_nu(x)*self.bias_2_nu(x),
            self.nu_min, self.nu_max, vec_func=True,
            tol=defaults.default_precision["global_precision"],
            rtol=defaults.default_precision["mass_precision"],
            divmax=defaults.default_precision["divmax"])
        self.bias_2_norm = -norm

    def _initialize_pp_gm(self):
        pp_gm_cent_array = numpy.zeros_like(self._ln_k_array)

        for idx in xrange(self._ln_k_array.size):
            pp_gm_cent = integrate.romberg(
                self._pp_gm_integrand, numpy.log(self.mass.nu_min),
                numpy.log(self.mass.nu_max), vec_func=True,
                tol=defaults.default_precision["halo_precision"],
                args=(self._ln_k_array[idx], "central_first_moment"))
            pp_gm_cent_array[idx] = pp_gm_cent/self.n_bar

        self._pp_gm_cent_spline = InterpolatedUnivariateSpline(
            self._ln_k_array, pp_gm_cent_array)

        pp_gm_sat_array = numpy.zeros_like(self._ln_k_array)

        for idx in xrange(self._ln_k_array.size):
            pp_gm_sat = integrate.romberg(
                self._pp_gm_integrand, numpy.log(self.mass.nu_min),
                numpy.log(self.mass.nu_max), vec_func=True,
                tol=defaults.default_precision["halo_precision"],
                args=(self._ln_k_array[idx], "satellite_first_moment"))
            pp_gm_sat_array[idx] = pp_gm_sat/self.n_bar

        self._pp_gm_sat_spline = InterpolatedUnivariateSpline(
            self._ln_k_array, pp_gm_sat_array)

        self._initialized_pp_gm = True

def complex_integral(func,a,b,args,intype='stupid'):
    real_func = lambda z: np.real(func(z,*args))
    imag_func = lambda z: np.imag(func(z,*args))
    if intype == 'quad':
        real_int = integ.quad(real_func,a,b)
        imag_int = integ.quad(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int[0] + 1j * imag_int[0]
    elif intype == 'quadrature':
        real_int = integ.quadrature(real_func,a,b)
        imag_int = integ.quadrature(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int[0] + 1j * imag_int[0]
    elif intype == 'romberg':
        real_int = integ.romberg(real_func,a,b)
        imag_int = integ.romberg(imag_func,a,b)
#        print(real_int)
#        print(imag_int)
        return real_int + 1j * imag_int
    elif intype == 'stupid':
        Npoints = 500
        z,dz = np.linspace(a,b,Npoints,retstep=True)
        real_int = np.sum(real_func(z))*dz
        imag_int = np.sum(imag_func(z))*dz
#        print(real_int)
#        print(imag_int)
        return real_int + 1j*imag_int

    def raw_correlation(self, r):
        """
        Compute the value of the correlation at array values r

        Args:
            r: float array of position values in Mpc/h
        """
        try:
            xi_out = numpy.empty(len(r))
            for idx, value in enumerate(r):
                xi_out[idx] = integrate.romberg(
                    self._correlation_integrand,
                    self._ln_k_min, self._ln_k_max, args=(value,),
                    vec_func=True,
                    tol=defaults.default_precision["global_precision"],
                    rtol=defaults.default_precision["corr_precision"],
                    divmax=defaults.default_precision["divmax"])
        except TypeError:
            xi_out = integrate.romberg(
                    self._correlation_integrand,
                    self._ln_k_min, self._ln_k_max, args=(r,), vec_func=True,
                    tol=defaults.default_precision["global_precision"],
                    rtol=defaults.default_precision["corr_precision"],
                    divmax=defaults.default_precision["divmax"])
        return xi_out

def dirichlet_integrate(alpha):
    normalizer = exp(sum(gammaln(alpha)) - gammaln(sum(alpha)))

    def f_recur(x, idx, upper, vals):
        if idx == 1:
            # base case.
            # set values for last two components
            vals[1] = x
            vals[0] = 1.0 - sum(vals[1:])
            # compute Dirichlet value
            print(vals.T, prod(vals ** (alpha - 1)) , normalizer, alpha)
            return prod(vals.T ** (alpha - 1)) / normalizer
        else:
            vals[idx] = x
            split = alpha[idx-1] / sum(alpha)
            if (split < upper - x):
                return romberg(f_recur, 0, split, args=(idx - 1, upper - x, vals), vec_func=False) + \
                    romberg(f_recur, split, upper - x, args=(idx - 1, upper - x, vals), vec_func=False)
            else:
                return romberg(f_recur, 0, upper - x, args=(idx - 1, upper - x, vals), vec_func=False)

    split = alpha[-1] / sum(alpha)
    print(alpha / sum(alpha))
    return romberg(f_recur, 0, split, args=(len(alpha) - 1, 1.0, zeros((len(alpha), 1), float64)), vec_func=False) + \
        romberg(f_recur, split, 1, args=(len(alpha) - 1, 1.0, zeros((len(alpha), 1), float64)), vec_func=False) 

	def __call__(self,Phi_Mz,p_Mz,par):
		#
		integrand = lambda M,z: Phi_Mz(M,z,par) * p_Mz(M,z) * self.dVdzdO(z)
		inner = lambda z: romberg(integrand,*self.Mrange,args=(z,),
		                          **self.int_kwargs)
		outer = romberg(inner,*self.zrange,**self.int_kwargs)
		return outer

    def correlation(self, l):
        """
        Compute the value of the correlation at array values theta

        Args:
            theta: float array of angular values in radians to compute the
                correlation
        """
        try:
            power = numpy.empty(len(l))
            for idx, value in enumerate(l):
                power[idx] = integrate.romberg(
                    self._correlation_integrand, 
                    self.kernel.chi_min, self.kernel.chi_max, args=(value,), 
                    vec_func=True,
                    tol=defaults.default_precision["global_precision"],
                    rtol=defaults.default_precision["corr_precision"],
                    divmax=defaults.default_precision["divmax"])
        except TypeError:
            power = integrate.romberg(
                self._correlation_integrand, 
                self.kernel.chi_min, self.kernel.chi_max, args=(l,), vec_func=True,
                tol=defaults.default_precision["global_precision"],
                rtol=defaults.default_precision["corr_precision"],
                divmax=defaults.default_precision["divmax"])
        return power

    def synphot (self, wlen, flam):
        """`wlen` and `flam` give a tabulated model spectrum in wavelength and f_λ
        units. We interpolate linearly over both the model and the bandpass
        since they're both discretely sampled.

        Note that quadratic interpolation is both much slower and can blow up
        fatally in some cases. The latter issue might have to do with really large
        X values that aren't zero-centered, maybe?

        I used to use the quadrature integrator, but Romberg doesn't issue
        complaints the way quadrature did. I should probably acquire some idea
        about what's going on under the hood.

        """
        from scipy.interpolate import interp1d
        from scipy.integrate import romberg

        d = self._ensure_data ()

        mflam = interp1d (wlen, flam,
                          kind='linear',
                          bounds_error=False, fill_value=0)

        mresp = interp1d (d.wlen, d.resp,
                          kind='linear',
                          bounds_error=False, fill_value=0)

        bmin = d.wlen.min ()
        bmax = d.wlen.max ()

        numer = romberg (lambda x: mresp (x) * mflam (x),
                         bmin, bmax, divmax=20)
        denom = romberg (lambda x: mresp (x),
                         bmin, bmax, divmax=20)
        return numer / denom

 def ri(u):
     ''' Epälineaarisen osan integrointi Gaussin kaavalla '''
     J = h/2 # Jakobiaani
     I = lambda k: N(k)*np.cos(np.sum(N(k)*u))*J # Integrandi
     I0 = lambda k: I(k)[0]
     I1 = lambda k: I(k)[1]
     #return 2.0*I(0.0)
     #return quad(I,-1,1)[0]
     return np.array([romberg(I0,-1,1),romberg(I1,-1,1)])

def romberg_slit_1d(q, width, height, form, pars):
    """
    Romberg integration for slit resolution.

    This is an adaptive integration technique.  It is called with settings
    that make it slow to evaluate but give it good accuracy.
    """
    from scipy.integrate import romberg, dblquad

    if any(k not in form.info['defaults'] for k in pars.keys()):
        keys = set(form.info['defaults'].keys())
        extra = set(pars.keys()) - keys
        raise ValueError("bad parameters: [%s] not in [%s]"%
                         (", ".join(sorted(extra)), ", ".join(sorted(keys))))

    if np.isscalar(width):
        width = [width]*len(q)
    if np.isscalar(height):
        height = [height]*len(q)
    _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars)
    _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars)
    # If both width and height are defined, then it is too slow to use dblquad.
    # Instead use trapz on a fixed grid, interpolated into the I(Q) for
    # the extended Q range.
    #_int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars)
    q_calc = slit_extend_q(q, np.asarray(width), np.asarray(height))
    Iq = eval_form(q_calc, form, pars)
    result = np.empty(len(q))
    for i, (qi, w, h) in enumerate(zip(q, width, height)):
        if h == 0.:
            r = romberg(_int_w, 0, w, args=(qi,),
                        divmax=100, vec_func=True, tol=0, rtol=1e-8)
            result[i] = r/w
        elif w == 0.:
            r = romberg(_int_h, -h, h, args=(qi,),
                        divmax=100, vec_func=True, tol=0, rtol=1e-8)
            result[i] = r/(2*h)
        else:
            w_grid = np.linspace(0, w, 21)[None,:]
            h_grid = np.linspace(-h, h, 23)[:,None]
            u = sqrt((qi+h_grid)**2 + w_grid**2)
            Iu = np.interp(u, q_calc, Iq)
            #print np.trapz(Iu, w_grid, axis=1)
            Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:,0])
            result[i] = Is / (2*h*w)
            """
            r, err = dblquad(_int_wh, -h, h, lambda h: 0., lambda h: w,
                             args=(qi,))
            result[i] = r/(w*2*h)
            """

    # r should be [float, ...], but it is [array([float]), array([float]),...]
    return result

    def __massEjected(self, m_min):
        """
        Return the mass integration of the mass ejected by the collapse of the
        star
        """
        if(self.imfType == "S"):
            return self.__massEjectedSalpeter(m_min)

        else:
            mEject = (romberg(self.__mPhi, m_min, self.__amsup1, tol=1.48e-04)
                 - romberg(self.__mrPhi, m_min, self.__amsup1, tol=1.48e-04)
                      )
        return mEject

def dL2(z,h):
    c=3.e5
    H0=100.*h
    if type(z) == ndarray:
	n=len(z)
	DL=zeros(n,'d')
	for i in range(n):
	    s=romberg(func,0.,z[i])#,tol=1.e-6)
	    DL[i]=c/H0*(1+z[i])*s #(Mpc/h)
    else:
	s=romberg(func,0.,z)#,tol=1.e-6)
	DL=c/H0*(1+z)*s #(Mpc/h)
    return DL

def dL(z,h):
    c=3.e5
    H0=100.*h

    try:#multiple objects
	n=len(z)
	DL=zeros(n,'d')
	for i in range(n):
	    s=romberg(func,0.,z[i])#,tol=1.e-6)
	    DL[i]=c/H0*(1+z[i])*s #(Mpc/h)
    except TypeError:
	s=romberg(func,0.,z)#,tol=1.e-6)
	DL=c/H0*(1+z)*s #(Mpc/h)
    return DL

    def _de_presure(self, redshift):
        dpresuredz = lambda z: (1 + self.w(z))/(1 + z)

        try:
            presure = numpy.empty(len(redshift))
            for idx, z in enumerate(redshift):
                presure[idx] = 3.0*integrate.romberg(
                    dpresuredz, 0, z, vec_func=True,
                    tol=defaults.default_precision["cosmo_precision"])
        except TypeError:
            presure = 3.0*integrate.romberg(
                dpresuredz, 0, redshift, vec_func=True,
                tol=defaults.default_precision["cosmo_precision"])
        return presure

def _test_one(name, p, w, tol):
    import scipy.integrate as sum
    # Check that the pdf approximately matchs the numerical integral
    # Check that integral(-inf, 0) of pdf sums to 0.5
    err = abs(sum.romberg(p.pdf, -20*w, 0, tol=1e-15) - 0.5)
    assert err < tol, "%s cdf(0) == 0.5 yields %g"%(name, err)

    # Check that integral(-inf, x) of pdf sums to cdf when x != 0
    err = abs(sum.romberg(p.pdf, -20*w, w/6, tol=1e-15) - p.cdf(w/6))
    assert err < tol, "%s cdf(w/6) == w/6 yields %g"%(name, err)

    # Check that P = cdf(ppf(P))
    P = 0.002
    err = abs(p.cdf(p.ppf(P)) - P)
    assert err < tol, "%s p(lo) = P yields %g"%(name, err)

def RhoCoul(r1,r2,theta,tau,xkappa,z,kmax,lmax,nmax,D):
    uklim = sqrt(float(kmax)/(tau*xkappa))
    rhocoul = 0.
    if D == 2:
        rhocoul = integrate.romberg(CGrand2D, 0., 2.*uklim, args=(r1,r2,theta,tau,xkappa,z,lmax))
    elif D == 3:
        rhocoul = integrate.romberg(CGrand3D, 0., 2.*uklim, args=(r1,r2,theta,tau,xkappa,z,lmax))
    if (z > 0.):
        return rhocoul
    else:
        if D == 2:
            rhocoul += Bound2D(r1,r2,theta,tau,xkappa,z,nmax)
        elif D == 3:
            rhocoul += Bound3D(r1,r2,theta,tau,xkappa,z,nmax)
    return rhocoul

def yieldcode(tmin, tmax, time, RBsol, fracinterp, polynomial, c, IMF):
    """
    The yieldcode module will use the step solution from the rombergquad
    module to calculate the yield in any time period, by calculating and
    excess outside the step ranges.
    """

    if tmin < time[-1]:
        tmin = time[-1]
        
    if tmax > time[0]:
        tmax = time[0]

    deriv = polynomial.deriv()    

    lowerindex = np.where(time<=tmax)[0][0]
    upperindex = np.where(time>=tmin)[0][-1]
   

    if lowerindex >= upperindex: # this is to avoid counting romberg intervals when the times don't actually cover any interval.
        m = lambda t: np.exp(polynomial(np.log(t)))
        y = lambda t: IMF(m(t))*c*deriv(np.log(t))*fracinterp(t)/t*m(t)    
        intsum = romberg(y, tmax, tmin)
        
    else:
        intervals = RBsol[lowerindex:upperindex] 
    
        intsum = np.sum(intervals)
    
        
        # Finding the excess integrals
        m = lambda t: np.exp(polynomial(np.log(t)))
        y = lambda t: IMF(m(t))*c*deriv(np.log(t))*fracinterp(t)/t*m(t) 
       
        t1 = time[lowerindex]   # corrected time is a decreasing vector
        t2 = tmax

        sol = romberg(y, t2, t1)
        intsum += sol
        
        
        t1 = tmin   
        t2 = time[upperindex]

        sol = romberg(y, t2, t1)
        intsum += sol
        
    return intsum

 def test_romberg_rtol(self):
     # Typical function with two extra arguments:
     def myfunc(x, n, z):       # Bessel function integrand
         return 1e19*cos(n*x-z*sin(x))/pi
     val = romberg(myfunc, 0, pi, args=(2, 1.8), rtol=1e-10)
     table_val = 1e19*0.30614353532540296487
     assert_allclose(val, table_val, rtol=1e-10)