algorithm - Python Numerical Integration for Volume of Region -

for program, need algorithm compute volume of solid. shape specified function that, given point p(x,y,z), returns 1 if p point of solid , 0 if p not point of solid.

i have tried using numpy using following test:

import numpy scipy.integrate import * def integrand(x,y,z):     if x**2. + y**2. + z**2. <=1.:         return 1.     else:         return 0. g=lambda x: -2. f=lambda x: 2. q=lambda x,y: -2. r=lambda x,y: 2. i=tplquad(integrand,-2.,2.,g,f,q,r) print 

but fails giving me following errors:

warning (from warnings module): file "c:\python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, integrationwarning) integrationwarning: maximum number of subdivisions (50) has been achieved. if increasing limit yields no improvement advised analyze integrand in order determine difficulties. if position of local difficulty can determined (singularity, discontinuity) 1 gain splitting interval , calling integrator on subranges. perhaps special-purpose integrator should used.

warning (from warnings module): file "c:\python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, integrationwarning) integrationwarning: algorithm not converge. roundoff error detected in extrapolation table. assumed requested tolerance cannot achieved, , returned result (if full_output = 1) best can obtained.

warning (from warnings module): file "c:\python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, integrationwarning) integrationwarning: occurrence of roundoff error detected, prevents requested tolerance being achieved. error may underestimated.

warning (from warnings module): file "c:\python27\lib\site-packages\scipy\integrate\quadpack.py", line 321 warnings.warn(msg, integrationwarning) integrationwarning: integral divergent, or convergent.

so, naturally, looked "special-purpose integrators", not find needed.

then, tried writing own integration using monte carlo method , tested same shape:

import random  # monte carlo method def get_volume(f,(x0,x1),(y0,y1),(z0,z1),prec=0.001,init_sample=5000):     xr=(x0,x1)     yr=(y0,y1)     zr=(z0,z1)     vdomain=(x1-x0)*(y1-y0)*(z1-z0)     def rand((p0,p1)):         return p0+random.random()*(p1-p0)     vol=0.     points=0.     s=0.  # sum part of variance of f     err=0.     percent=0     while err>prec or points<init_sample:         p=(rand(xr),rand(yr),rand(zr))         rpoint=f(p)         vol+=rpoint         points+=1         s+=(rpoint-vol/points)**2         if points>1:             err=vdomain*(((1./(points-1.))*s)**0.5)/(points**0.5)         if err>0:             if int(100.*prec/err)>=percent+1:                 percent=int(100.*prec/err)                 print percent,'% complete\n  error:',err     print int(points),'points used.'     return vdomain*vol/points f=lambda (x,y,z): ((x**2)+(y**2)<=4.) , ((z**2)<=9.) , ((x**2)+(y**2)>=0.25) print get_volume(f,(-2.,2.),(-2.,2.),(-2.,2.)) 

but works slowly. program using numerical integration 100 times or so, , doing on larger shapes, take minutes if not hour or 2 @ rate goes now, not mention want better precision 2 decimal places.

i have tried implementing miser monte carlo method, having difficulties , i'm still unsure how faster be.

so, asking if there libraries can asking, or if there better algorithms work several times faster (for same accuracy). suggestions welcome, i've been working on quite while now.

edit:

if cannot working in python, open switching other language both compilable , has relatively easy gui functionality. suggestions welcome.

your function not continuous function, think it's difficult integration.

how about:

import numpy np  def sphere(x,y,z):     return x**2 + y**2 + z**2 <= 1  x, y, z = np.random.uniform(-2, 2, (3, 2000000))  sphere(x, y, z).mean() * (4**3), 4/3.0*np.pi 

output:

(4.1930560000000003, 4.1887902047863905) 

or vtk:

from tvtk.api import tvtk  n = 151 r = 2.0 x0, x1 = -r, r y0, y1 = -r, r z0, z1 = -r, r x,y,z = np.mgrid[x0:x1:n*1j, y0:y1:n*1j, z0:z1:n*1j] s = sphere(x, y, z)  img = tvtk.imagedata(spacing=((x1-x0)/(n-1), (y1-y0)/(n-1), (z1-z0)/(n-1)),                       origin=(x0, y0, z0), dimensions=(n, n, n))  img.point_data.scalars = s.astype(float).ravel()  blur = tvtk.imagegaussiansmooth(input=img)  blur.set_standard_deviation(1) contours = tvtk.contourfilter(input = blur.output)   contours.set_value(0, 0.5) mp = tvtk.massproperties(input = contours.output) mp.volume, mp.surface_area 

output:

4.186006622559839, 12.621690438955586