def m_D(a_x,a_y,a_z):
    '''calculates the dynamic geometrical tensor of the ellipsoid:
       to be computed once for every geometry

    Parameters
    ----------
    'a_x,a_y,a_z' = three axes of the ellipsoid (in nm)

    Returns
    -------
    'D' = dynamic geometrical tensor of the ellipsoid'''

    D_x_int,err = sp_i.nquad(lambda t,z: f_Dx(t,z,a_x,a_y,a_z),
                             [[0,2.0*np.pi],[0.0,1.0]])
    D_y_int,err = sp_i.nquad(lambda t,z: f_Dy(t,z,a_x,a_y,a_z),
                             [[0,2.0*np.pi],[0.0,1.0]])
    D_z_int,err = sp_i.nquad(lambda t,z: f_Dz(t,z,a_x,a_y,a_z),
                             [[0,2.0*np.pi],[0.0,1.0]])
    D_y = (0.75/np.pi)*D_y_int
    D_x = (0.75/np.pi)*D_x_int
    D_z = (0.75/np.pi)*D_z_int

    D = np.array([[D_x,0.0,0.0],[0.0,D_y,0.0],[0.0,0.0,D_z]])

    return D

		def Integrandlosv(r):   #It's using r, because I define below has limits in r, not x
                        xi = r/rs
                        def sigr(vt,vr):
                                return f(xi, vr, vt, param) * vt * (vr*vr) * 2*np.pi 

			def sigt(vt,vr):
                                return f(xi, vr, vt, param) * vt * (vt*vt) * 2*np.pi
			
			def rhofunc2(vt,vr):
                                return f(xi, vr, vt, param) * vt * 2*np.pi		


                        vmax = vesc(xi,[rlim, Vmax, rmax, alpha, beta, gamma])
                        def bounds_vr():
                                return [0, vmax]  #should be +/-vmax, but the function is symmetric, so.

                        def bounds_vt(vr):
                                lim = (vmax*vmax - vr*vr)**0.5
                                return [0, lim]

                        rhosigr2 = integrate.nquad(sigr, [bounds_vt, bounds_vr])[0]
			rhosigt2 = integrate.nquad(sigt, [bounds_vt, bounds_vr])[0]
			rhor     = integrate.nquad(rhofunc2, [bounds_vt, bounds_vr])[0]
                        zz = (xi*xi - Xi*Xi)
                        dz = rs*xi / np.sqrt(zz)
			result = ( zz*rhosigr2 + 0.5*Xi*Xi*rhosigt2 ) * dz / (xi*xi)  #!!! 0.5 or no???
                        return 2*result

 def expected_value( self, x):
     """
     Currently, this does the whole sum/integral over the cube support of Z.
     We may be able to improve this by taking into account how the joint
     and conditionals factorize, and/or finding a more efficient support.
     
     This should be reasonably fast for |Z| <= 2 or 3, and small enough discrete
     variable cardinalities.  It runs in O(n_1 n_2 ... n_k) in the cardinality of
     the discrete variables, |Z_1| = n_1, etc.  It likewise runs in O(V^n) for n
     continuous Z variables.  Factorizing the joint/conditional distributions in
     the sum could linearize the runtime.
     """
     causal_effect = 0.
     x = x[self.causes]
     if self.discrete_Z:
         discrete_variable_ranges = [ xrange(*(int(self.support[variable][0]), int(self.support[variable][1])+1)) for variable in self.discrete_Z]
         for z_vals in itertools.product(*discrete_variable_ranges):
             z_discrete = pd.DataFrame({k : [v] for k, v in zip(self.discrete_Z, z_vals)})
             if self.continuous_Z:
                 continuous_Z_ranges = [self.support[variable] for variable in self.continuous_Z]
                 args = z_discrete.join(x).values[0]
                 causal_effect += nquad(self.expectation_integration_function,continuous_Z_ranges,args=args)[0]
             else:
                 z_discrete = z_discrete[self.admissable_set]
                 exog_predictors = x.join(z_discrete)[self.conditional_density_vars]
                 causal_effect += self.conditional_expectation.fit(data_predict=exog_predictors.values)[0] * self.density.pdf(data_predict=z_discrete.values)
         return causal_effect
     elif self.continuous_Z:
         continuous_Z_ranges = [self.support[var] for var in self.continuous_Z]
         causal_effect, error = nquad(self.expectation_integration_function,continuous_Z_ranges,args=tuple(x.values[0]))
         return causal_effect
     else:
         return self.conditional_expectation.fit(data_predict=x[self.causes])[0]

def main():
	n=6
	liste=[q for q in xrange(1,n+1)]
	
	func= lambda x1,x2,x3: function([x1,x2,x3])
	print integrate.nquad(func,[[0,1],[0,1],[0,1]], opts=[{"epsabs": 1.49e-20},{"epsabs": 1.49e-20},{"epsabs": 1.49e-20}])
	
	
	
	
	func2= lambda y1,y2,y3: another_function(**{"eins": y1,"zwei": y2,"drei": y3})
	print integrate.nquad(func2,[[0,1],[0,1],[0,1]], opts=[{"epsabs": 1.49e-20},{"epsabs": 1.49e-20},{"epsabs": 1.49e-20}])

    def testZernikeAnnularEval(self):

        # Obscuration
        e = 0.61

        # Calculate the radius
        dd = np.sqrt(self.xx**2 + self.yy**2)

        # Define the invalid range
        idx = (dd > 1) | (dd < e)

        # Create the Zernike terms
        Z = np.zeros(22)

        # Generate the map of Z12
        Z[11] = 1

        # Calculate the map of Zernike polynomial
        Zmap = ZernikeAnnularEval(Z, self.xx, self.yy, e)
        Zmap[idx] = np.nan

        # Put the elements to be 0 in the invalid region
        Zmap[np.isnan(Zmap)] = 0

        # Check the normalization for Z1 - Z28
        e = 0.61
        ansValue = np.pi*(1-e**2)
        for ii in range(28):
            Z = np.zeros(28)
            Z[ii] = 1

            normalization = nquad(self._genNormalizedFunc,
                                  [[e, 1], [0, 2*np.pi]], args=(Z, e))[0]

            self.assertAlmostEqual(normalization, ansValue)

        # Check the orthogonality for Z1 - Z28
        for jj in range(28):
            Z1 = np.zeros(28)
            Z1[jj] = 1
            for ii in range(28):
                if (ii != jj):
                    Z2 = np.zeros(28)
                    Z2[ii] = 1

                    orthogonality = nquad(self._genOrthogonalFunc,
                                          [[e, 1], [0, 2*np.pi]],
                                          args=(Z1, Z2, e))[0]

                    self.assertAlmostEqual(orthogonality, 0)

 def _calc_spec_pgamma(self, Eeminus):
     Eeminus = Eeminus.to('eV').value
     x_range = [0, 1]
     eta_range = [0.66982, 31.3]
     spec_hi = self._mpc2 * nquad(self._H_integrand, [x_range, eta_range],
                                  args=[Eeminus], )[0]
     return spec_hi.value

def selectModel(npafIndepVar, npafData, funErrorPDF, lfunPriorParamPDF,
                lllfPriorParamBounds, lfunModels, lfPriorOdds=[]):
    '''
        Select between two models using a naive quadrature method for
        integration of probability distributions. Inputs are the data to model,
        the error pdf for the data, the parameter probability density, a list
        of 2 element lists of the [lower bound, upper bound] of the arguments
        of parameter distribution, the list of models to select from and a list
        of lists of 2 element lists of the model parameter bounds. For discrete
        parameters, select between functions defined with the different values
        of the discrete parameters. The outputs are the index of the best
        model and a list of the odds ratios for each model relative to the best
        model.
    '''
    from scipy.integrate import nquad
    import numpy as np
    npaOR = []
    if lfPriorOdds == []:
        lfPriorOdds = np.ones(len(lfunModels))
    nIdx = 0
    for funParamPDF, llfParamBounds, funModel in zip(lfunPriorParamPDF,
                                                     lllfPriorParamBounds,
                                                     lfunModels):
        funIntegrand = lambda *x: (funParamPDF(x) *
                                   funErrorPDF(npafData -
                                   funModel(npafIndepVar, x)))
        npaOR.append(nquad(funIntegrand, llfParamBounds)*lfPriorOdds[nIdx])
        nIdx += 1
    npaOR = np.array(npaOR)
    npaOR = npaOR/np.max(npaOR)
    nBestModel = np.argmin(npaOR)
    return nBestModel, npaOR

    def test_matching_tplquad(self):
        def func3d(x0, x1, x2, c0, c1):
            return x0 ** 2 + c0 * x1 ** 3 - x0 * x1 + 1 + c1 * np.sin(x2)

        res = tplquad(func3d, -1, 2, lambda x: -2, lambda x: 2, lambda x, y: -np.pi, lambda x, y: np.pi, args=(2, 3))
        res2 = nquad(func3d, [[-np.pi, np.pi], [-2, 2], (-1, 2)], args=(2, 3))
        assert_almost_equal(res, res2)

def simu_YC(Alpha,beta,Z,U,V1,V2,T,Sigma):
    N=len(T)
    D=len(U[0])
    Y=np.zeros(N,float)
    for i in range(N):
        Y[i]=Z[i].dot(Alpha)
        for j in range(D):
            def Fbeta(t):
                return U[i,j]+np.sum(V1[i,j,:]*np.sin( (np.arange(10)+1.)*2.*ma.pi*t/100.)+V2[i,j,:]*np.cos((np.arange(10)+1.)*2.*ma.pi*t/100.))*beta.val([t,T[i]]);
            Y[i]= Y[i] + integrate.nquad(Fbeta,[[0,T[i]]],opts=[{'epsabs':5e-04}])[0] / T[i]
    # Erreur sur le label
    if Sigma>0:
        Y = Y + np.random.normal(0,Sigma,N)
    # Labels censurés et temps de censures
    C=np.random.normal(np.mean(Y)+np.std(Y)/2,np.sqrt(np.std(Y)),N)
    Yc=np.zeros(N,float)
    Delta=np.zeros(N)
    for i in range(N):
        if C[i]==min(C[i],Y[i]):
            Yc[i]=C[i]
            Delta[i]=0
        else:
            Yc[i]=Y[i]
            Delta[i]=1
    return (Y,C,Yc,Delta);

 def integral(self, f):
     """
     Returns the integral of `f` with respect to the currwnt measure
     over the support.
     """
     from types import FunctionType
     m = 0
     if self.DomType == "set":
         if type(f) not in [dict, FunctionType]:
             raise Exception(
                 "The integrand must be a `function` or a `dict`")
         if type(f) == dict:
             for p in self.supp:
                 if p in f:
                     m += self.weight[p] * f[p]
         else:
             for p in self.supp:
                 try:
                     m += self.weight[p] * f(p)
                 except:
                     pass
     else:
         if type(f) != FunctionType:
             raise Exception("The integrand must be a `function`")
         from scipy import integrate
         fw = lambda *x: f(*x) * self.weight(*x)
         m = integrate.nquad(fw, self.supp)[0]
     return m

 def f_int():
     y = 0
     if (delta != 0):
         y, err = integrate.nquad(f,  [[-delta, delta], [-delta, delta], [-delta, delta], [-delta, delta]])
     else:
         y = f(0)
     return 1 / 16.0 / delta**4 * y

def SigR(param):
	a,d,e, Ec, rlim, b, q, Jb, Vmax, rmax, alpha, beta, gamma = param
        rs = rmax/2.16           # rmax=2.16*rs
	xlim = rlim/rs
	
	def rhoRI(theta, v, x, X):
		z   = (x*x - X*X +10**-10)**.5
		aux = 2 * ftheta(x, v, theta, param) * pow(v,2) * np.sin(theta) * 2 * np.pi  * rs
		return aux * x/z

	def bounds_theta(v,x):
                return [0, np.pi]
        def bounds_v(x):
                vmax = vesc( x, [rlim, Vmax, rmax, alpha, beta, gamma] )                        
		return [0, vmax]
	#def bounds_r(X):
        #        return [X, xlim]

	#def result(Xi):
	#	ans = integrate.nquad(rhoRI, [bounds_theta, bounds_v, bounds_r], args=(Xi,))[0]
	#	return ans

	Xarr = np.linspace(0, xlim, 10)
	projlist = []
	for Xi in Xarr:
		def rhoRI2(theta,v,x):
			return rhoRI(theta,v,x,Xi)
		ans = integrate.nquad(rhoRI2, [bounds_theta, bounds_v, [Xi,xlim]])[0]
		projlist.append(ans)
		print 'SigR Xi: ', Xi
	proj = np.array(projlist)

	return Xarr*rs, proj/max(proj)

    def test_variable_limits(self):
        scale = .1
        def func2(x0, x1, x2, x3, t0, t1):
            return x0*x1*x3**2 + np.sin(x2) + 1 + (1 if x0 + t1*x1 - t0 > 0 else 0)
        def lim0(x1, x2, x3, t0, t1):
            return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
                    scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
        def lim1(x2, x3, t0, t1):
            return [scale * (t0*x2 + t1*x3) - 1,
                    scale * (t0*x2 + t1*x3) + 1]
        def lim2(x3, t0, t1):
            return [scale * (x3 + t0**2*t1**3) - 1,
                    scale * (x3 + t0**2*t1**3) + 1]
        def lim3(t0, t1):
            return [scale * (t0 + t1) - 1, scale * (t0 + t1) + 1]
        def opts0(x1, x2, x3, t0, t1):
            return {'points':[t0 - t1*x1]}
        def opts1(x2, x3, t0, t1):
            return {}
        def opts2(x3, t0, t1):
            return {}
        def opts3(t0, t1):
            return {}

        res = nquad(func2, [lim0, lim1, lim2, lim3], args=(0,0),
                    opts=[opts0, opts1, opts2, opts3])
        assert_quad(res, 25.066666666666663)

    def test_square_aliased_ranges_and_opts(self):
        def f(y, x):
            return 1.0

        r = [-1, 1]
        opt = {}
        assert_quad(nquad(f, [r, r], opts=[opt, opt]), 4.0)

def psi_Kolonderski_full(lp,ls,li,wp,ws,wi,axp,axs,axi):
# Integrale come somma; non funziona
#    (ksx,dksx) = (ksy,dksy) = (kix,dkix) = (kiy,dkiy) = np.linspace(-1000000,1000000,50,retstep=True)
#    
#    f = func_Kolonderski(ksx[:,None,None,None],ksy[None,:,None,None],kix[None,None,:,None],kiy[None,None,None,:],lp,ls,li,wp,ws,wi,axp,axs,axi)
#    return np.sum(np.sum(np.sum(np.sum(f))))*dksx*dksy*dkix*dkiy
# Integrale usando le funzioni di scipy: eterno
    return integ.nquad(func_Kolonderski,[[-np.inf,np.inf],[-np.inf,np.inf],[-np.inf,np.inf],[-np.inf,np.inf]],args=(lp,ls,li,wp,ws,wi,axp,axs,axi))

    def test_matching_quad(self):
        def func(x):
            return x**2 + 1

        res, reserr = quad(func, 0, 4)
        res2, reserr2 = nquad(func, ranges=[[0, 4]])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2)

    def test_matching_dblquad(self):
        def func2d(x0, x1):
            return x0**2 + x1**3 - x0 * x1 + 1

        res, reserr = dblquad(func2d, -2, 2, lambda x:-3, lambda x:3)
        res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2)

def multi_integrate(func, options):
    r = c.parse_range(options['range'])

    try:
        result = i.nquad(func, r)
        return c.clean_integrate(result)
    except Exception as e:
        return str(e)

    def test_densities(self):
        causes = ['c']
        effects = ['d']
        admissable_set = ['a']
        variable_types={'a': 'c','b': 'c','c': 'c','d' : 'c'}
        effect = CausalEffect(self.X,causes,effects,admissable_set,variable_types)
        density =  lambda x: effect.density.pdf( data_predict=[x])
        integral = nquad( density, [effect.support[d_var] for d_var in admissable_set])[0]
        print integral
        assert(abs(integral - 1.) < TOL)

        x_vals = [np.mean(effect.support[var]) for var in causes]
        z_vals = [np.mean(effect.support[var]) for var in admissable_set]
        density = lambda x: effect.conditional_density.pdf(endog_predict=[x], exog_predict=x_vals + z_vals)
        integral = nquad(density, [effect.support[d_var] for d_var in effects])[0]
        print x_vals, z_vals,integral
        assert(abs(integral - 1.) < TOL)

def Posterior():
	func= lambda s,b,mumu: norm.pdf(s,mean_s,sigma_s)*norm.pdf(b,mean_b,sigma_b)*ExtendedPoisson(N,mumu*s+b)
	func2= lambda mumumu: integrate.nquad(func,[[mean_s - 5*sigma_s,mean_s + 5*sigma_s],[mean_b - 5*sigma_b,mean_b + 5*sigma_b]],[mumumu])[0]#/integrate.nquad(func,[[mean_s - 5*sigma_s,mean_s + 5*sigma_s],[mean_b - 5*sigma_b,mean_b + 5*sigma_b],[1e-15,10]])[0]

	x=arange(0,5.0,0.2)
	y=map(func2,x)
	plt.plot(x,y)
	plt.show()