"""Multi-dimensional integration.
Integrates a function over a hypercube in multiple dimensions.
"""
__all__ = ['dcuhre']
import warnings
try:
import _dcuhre
except ImportError, err:
warnings.warn("Could not import _dcuhre. Skipping...")
_dcuhre = None
import numpy
if _dcuhre:
[docs] def dcuhre(f,a,b,epsabs=1.49e-08,epsrel=1.49e-08,
full_output=False,minpts=0,maxpts=100000,
key=0,nproc=1,nw=None):
"""
Description:
Return the integral of f over the hypercube a <= x <= b.
Inputs:
f -- f(x) should return a one dimensional vector (ndarray for
performance) of the function(s) evaluated at the point x.
a -- Lower limits of integration (2 <= ndim = len(a) <= 15)
b -- Upper limits of integration (2 <= ndim = len(b) <= 15)
epsabs -- Desired absolute tolerance
epsabs -- Desired relative tolerance
minpts -- Minimum number of points to sample function at
maxpts -- Maximum number of function evaluations
key -- Key to selected local integration rule.
key = 0 is the default value.
For ndim = 2 the degree 13 rule is selected.
For ndim = 3 the degree 11 rule is selected.
For ndim > 3 the degree 9 rule is selected.
key = 1 gives the user the 2 dimensional degree 13
integration rule that uses 65 evaluation
points.
key = 2 gives the user the 3 dimensional degree 11
integration rule that uses 127 evaluation
points.
key = 3 gives the user the degree 9 integration rule.
key = 4 gives the user the degree 7 integration rule.
This is the recommended rule for problems that
require great adaptivity.
nproc -- Number of processors (if DCUHRE compiled with
multiple processor support)
nw -- Defines the length of the working array. Should be
automatically determined (see below).
full_output -- non-zero to return a dictionary of integration
information. If non-zero, warning messages are
also suppressed and the message is appended to
the output tuple.
The result hopefully satisfies for each component k of f the
following claim for accuracy:
abs(ans - exact) <= max(epsabs,epsrel*abs(ans))
Based on the FORTRAN code DCUHRE which is a driver for the
integration routine DADHRE, which repeatedly subdivides the region
of integration and estimates the integrals and the errors over the
subregions with greatest estimated errors until the error request
is met or maxpts function evaluations have been used.
Let ndim = len(a) == len(b) be the dimension of the domain
hypercube:
- For ndim = 2 the default integration rule is of degree 13 and
uses 65 evaluation points.
- For ndim = 3 the default integration rule is of degree 11 and
uses 127 evaluation points.
- For ndim greater then 3 the default integration rule
is of degree 9 and uses num evaluation points where
num = 1 + 4*2*ndim + 2*ndim*(ndim-1) + 4*ndim*(ndim-1) +
4*ndim*(ndim-1)*(ndim-2)/3 + 2**ndim
The degree 9 rule may also be applied for ndim = 2 and ndim = 3.
A rule of degree 7 is available in all dimensions. The number of
evaluation points used by the degree 7 rule is
num = 1 + 3*2*ndim + 2*ndim*(ndim-1) + 2**ndim
The number of function evaluations over each subregion
is num:
if (key == 0 or key == 1) and (ndim == 2):
num = 65
elif (key == 0 or key == 2) and (ndim == 3):
num = 127
elif (key == 0 and ndim > 3) or (key == 3):
num = 1 + 4*2*ndim + 2*ndim*(ndim-1) + 4*ndim*(ndim-1) +\
4*ndim*(ndim-1)*(ndim-2)/3 + 2**ndim
elif (key == 4):
num = 1 + 3*2*ndim + 2*ndim*(ndim-1) + 2**ndim
maxpts >= 3*num and maxpts >= minpts
for 3 < ndim < 13 the minimum values for maxpts are:
ndim = 4 5 6 7 8 9 10 11 12
key = 3: 459 819 1359 2151 3315 5067 7815 12351 20235
key = 4: 195 309 483 765 1251 2133 3795 7005 13299
Let maxsub denote the maximum allowed number of subregions for the
given values of maxpts, key and ndim.
maxsub = (maxpts-num)/(2*num) + 1
nw >= maxsub*(2*ndim + 2*numfun + 2) + 17*numfun + 1
For efficient execution on parallel computers nw should be chosen
greater or equal to maxsub*(2*ndim+2*numfun+2)+17*numfun*mdiv+1
where mdiv is the number of subregions that are divided in each
subdivision step. mdiv is default set internally in DCUHRE equal
to 1. For efficient execution on parallel computers with nproc
processors mdiv should be set equal to the smallest integer such
that (2*mdiv)%nproc == 0.
When DCUHRE computes estimates to a vector of integrals, all
components of the vector are given the same treatment. That is,
I(f[j]) and I(f[k]) for j != k, are estimated with the same
subdivision of the region of integration. For integrals with
enough similarity, we may save time by applying DCUHRE to all
integrands in one call. For integrals that vary continuously as
functions of some parameter, the estimates produced by DCUHRE will
also vary continuously when the same subdivision is applied to all
components. This will generally not be the case when the different
components are given separate treatment.
On the other hand this feature should be used with caution when
the different components of the integrals require clearly
different subdivisions.
Credit:
Jarle Berntsen, The Computing Centre,
University of Bergen, Thormohlens gt. 55,
N-5008 Bergen, Norway
Phone.. 47-5-544055
Email.. jarle@eik.ii.uib.no
Terje O. Espelid, Department of Informatics,
University of Bergen, Thormohlens gt. 55,
N-5008 Bergen, Norway
Phone.. 47-5-544180
Email.. terje@eik.ii.uib.no
Alan Genz, Computer Science Department, Washington State
University, Pullman, WA 99163-1210, USA
Phone.. 509-335-2131
Email.. acg@cs2.cs.wsu.edu
>>> def f(x):
... ndim = len(x)
... nf = ndim + 1
... s = 0
... for n in xrange(ndim):
... s = s + (n+1)*x[n]**2
... fv = numpy.zeros(nf)
... fv[0] = numpy.exp(-s/2.0)
... for n in xrange(ndim):
... fv[n+1] = x[n]*fv[0]
... return fv
>>> ndim = 5
>>> a = numpy.zeros(ndim,float)
>>> b = numpy.ones(ndim,float)
>>> abseps = 0
>>> releps = 0.001
>>> ans = dcuhre(f=f,a=a,b=b,epsabs=abseps,epsrel=releps)
Success!
>>> def f(x):
... return sum(x*x)
>>> Ndim = 4.0
>>> abs_tol = 0.0001
>>> rel_tol = 0.001
>>> ans = dcuhre(f,a=numpy.zeros(Ndim),b=numpy.ones(Ndim),\
epsabs=abs_tol,epsrel=rel_tol)
Success!
>>> abs(ans[0] - Ndim/3.0) < abs_tol
True
>>> abs(ans[0] - Ndim/3.0)/Ndim*3.0 < rel_tol
True
"""
ndim = len(a)
assert len(b) == ndim
assert a.shape == b.shape
fa = f(a)
if numpy.isscalar(fa):
numfun = 1
else:
numfun = len(fa)
if (key == 0 or key == 1) and (ndim == 2):
num = 65
elif (key == 0 or key == 2) and (ndim == 3):
num = 127
elif (key == 0 and ndim > 3) or (key == 3):
num = 1 + 4*2*ndim + 2*ndim*(ndim-1) + 4*ndim*(ndim-1) +\
4*ndim*(ndim-1)*(ndim-2)/3 + 2**ndim
elif (key == 4):
num = 1 + 3*2*ndim + 2*ndim*(ndim-1) + 2**ndim
else:
raise ValueError("key must be in {0,1,2,3,4}.")
assert maxpts >= 3*num and maxpts >= minpts
maxsub = (maxpts-num)/(2*num) + 1
assert 0 < nproc, isinstance(nproc,int)
mdiv = 1
while (2*mdiv)%nproc is not 0:
mdiv = mdiv + 1
nw = maxsub*(2*ndim+2*numfun+2)+17*numfun*mdiv+1
(res,abs_err,neval,ifail) = \
_dcuhre.dcuhre(numfun=numfun,a=a,b=b,
minpts=minpts,maxpts=maxpts,
funsub=f,
epsabs=epsabs,epsrel=epsrel,
nw=nw,key=key)
rel_err = abs(abs_err/res)
Message = {0:"Success!",
1:"Premature termination: try increasing maxpts.\n"+\
"Absolute Error: "+`abs_err`+" "+\
"(Request: epsabs = "+`epsabs`+")\n"+\
"Relative Error: "+`rel_err`+" "+\
"(Request: epsrel = "+`epsrel`+")",
2:"Invalid key = "+`key`+": (0,1,2,3,4 only)",
3:"Invalid ndim = "+`ndim`+": (2 to 15 dimensions only)",
4:"key = 1 only valid for ndim = 2 dimensions",
5:"key = 2 only valid for ndim = 3 dimensions",
6:"Invalid numfun (Internal error in python "+\
"code. Please inform maintainer...)",
7:"Invalid integration region (volume == 0)",
8:"Invalid maxpts < 3*num (please increase "+\
"maxpts >= "+`3*num`+")",
9:"Invalid maxpts < minpts (please increase "+\
"maxpts >= "+`minpts`+")",
10:"Invalid tolrances (must both be greater than 0).\n"+\
"epsabs = "+`epsabs`+", epsrel = "+`epsrel`,
11:"Workspace nw = "+`nw`+" too small",
12:"Invalid restar (Internal error in python "+\
"code. Please inform maintainer...)"}
result = {"Message":Message[ifail],
"Result":res,
"Estimated Error":abs_err,
"Function Evaluations":neval}
if full_output:
return result
else:
print result['Message']
return res
def subregions(a,b):
"""Iterator over each subregion of the region a<=x<=b.
>>> for (a,b) in subregions([0,0,0],[2,2,2]):
... print a,b
[ 0. 0. 0.] [ 1. 1. 1.]
[ 0. 0. 1.] [ 1. 1. 2.]
[ 0. 1. 0.] [ 1. 2. 1.]
[ 0. 1. 1.] [ 1. 2. 2.]
[ 1. 0. 0.] [ 2. 1. 1.]
[ 1. 0. 1.] [ 2. 1. 2.]
[ 1. 1. 0.] [ 2. 2. 1.]
[ 1. 1. 1.] [ 2. 2. 2.]
"""
dx2 = (numpy.array(b)-numpy.array(a))/2.0
Ndim = len(a)
n = 0
for n in xrange(2**Ndim):
inds = numpy.array(list(numpy.binary_repr(n).rjust(Ndim,'0'))
).astype(float)
yield (a+dx2*inds,a+dx2+dx2*inds)
return
def integrate(f,a,b,abstol=0.01):
"""Integrate the function f(x) over the hypercube a<=x<=b.
The function f(x) should broadcast well over an array of points.
>>> def f(x):
... return sum(x*x)
>>> Ndim = 3
>>> tol = 0.1
>>> ans = integrate(f=f,a=numpy.zeros(Ndim),b=numpy.ones(Ndim),abstol=tol)
>>> abs(ans - Ndim/3.0) < tol
True
"""
Ndim = len(a)
xmid = (a+b)/2.0
dx = b-a
fmid = f(xmid)
F = fmid*numpy.prod(dx)
if abs(F) < abstol:
# Return value if integral is small enough.
return F
else:
F = 0.0
for (an,bn) in subregions(a,b):
F = F + integrate(f,an,bn,abstol)
return F