"""Integration made easy.
Potentially useful formula (from the method of exhaustion):
integral(: f(x).dx ←x; a ≤x&leb :)
= sum(: 2**(1+n) * sum(: (-1)**m f(a +(b-a)*(m+1)/2**(1+n)) ←m :2**(1+n) -2)
←n :{naturals})*(b-a)
see: http://en.wikipedia.org/wiki/Method_of_exhaustion
$Id: integrate.py,v 1.5 2007/03/24 22:42:21 eddy Exp $
"""
class Integrator:
"""Base class for integrators.
Provides crude trapezium rule integrator, which derived classes might wish
to over-ride. More importantly, defines an API for integrators:
measure(func) -- integrator scaling self.integrand by func
between(start, stop) -- integrates over an interval
before(stop) -- as between, but with minus infinity as start
beyond(start) -- as between, but with plus infinity as stop
Both before() and beyond() presume that they're integrating a `tail' of a
distribution: to compute an integral between one side of a distribution's
mode and infinity on the other side, it is probably wisest to call between()
on the given bound and a more-or-less arbitrary position on the same side of
the mode as the infinity in question, then to use before() or beyond() to
compute the remainder, using between()'s output as offset (see below).
All three integrator methods take an optional argument, test, following the
given required arguments; and an optional keyword argument, offset, after
that. These are used in deciding when the integral has been determined
accurately enough: after each iteration, the integrator calls test with the
change in its estimate of the integral as first argument and, as second
argument, the new estimate, optionally with offset added to it. If the
change is small enough that further refinement is a waste of time, test
should return true. The integrator will then return the given estimate
added to an error bar whose width is the last change.
By default, the code is geared up to deal with values of class Quantity (see
basEddy.quantity) and, indeed, the integrators return Quantity()s. If the
first estimate at the integral (plus optional offset) has a positive .width
attribute, it is presumed to be a Quantity() and the default test simply
compares the change in estimate with the .width of the new estimate; if the
change is within this existing error bar, it is taken to be small enough.
In the absence of .width, or if the .width is not positive, the default test
simply looks to see whether the change in estimate is smaller than 1e-6
times the estimate (plus optional offset); in the absence of an offset, this
last test will work poorly if the integral should yield zero. """
def __init__(self, func, width=None):
"""Initialises an integrator.
Required first argument is the function to be integrated.
This will be stored as self.integrand.
Optional second argument, width - only needed if .beyond() or .before()
is liable to be called with bound zero - is a `unit' change in input.
It should ideally be approximately the difference between highest and
lowest inputs for which the integrand differs significantly from 0;
e.g. if func describes a Gaussian distribution with variance 1, a width
of about 5 would be prudent. """
self.integrand = func
if width is not None: self.__unit = abs(width)
def measure(self, func, width=None):
"""New integrator scaling self.integrand pointwise by a given function.
Required argument, func, is a function accepting the same inputs as
self.integrand(). Optional argument, width, is as for the Integrator
constructor; if omitted, the value used when self was constructed (if
any) is used.
Returns an Integrator whose integrand is (: func(x) * self.integrand(x)
←x :) which can be construed as integrating func using self as
measure; or as integrating self using func as measure. If self is a
probability distribution, its moments can be computed using func()s of
form (: x**i ←x :) for i = 1, 2, 3 ..."""
if width is None:
try: width = self.__unit
except AttributeError: pass
return Integrator(lambda x, f=func, i=self.integrand: f(x) * i(x), width)
def between(self, start, stop, test=None, offset=None):
"""Integral over a range.
Two required arguments, start and stop, give the bounds of the range.
Also accepts the usual optional tolerance specifiers, test and offset:
see class doc for details. """
return self.__between(start,
(stop - start) * 1., # avoid integral type gotchas ...
(self.integrand(start) + self.integrand(stop)) * .5,
test, offset)
def before(self, stop, test=None, offset=None):
"""Integration from minus infinity.
Equivalent to between() with start set to a value less than any at which
self's distribution is distinguishable from zero. Required argument is
the upper bound of the integral; also accepts the usual optional
tolerance specifiers, test and offset: see class doc for details. """
return -self.__outwards(stop, -self.__step(stop), test, offset)
def beyond(self, start, test=None, offset=None):
"""Integration to plus infinity.
As before() but required argument is the lower bound of integration, the
upper bound being greater than any value at which self's distribution is
distinguishable from zero. """
return self.__outwards(start, self.__step(start), test, offset)
# hairy implementation follows: no further exports.
from study.value.quantity import tophat
def __blur(mid, spread, H=tophat): return mid + H * spread
del tophat
# __blur and __gettest will be del'd shortly ... they're *not* methods
def __gettest(eg,
microclose = lambda d, n: abs(d) <= 1e-6 * abs(n),
bywidth = lambda d, n: abs(d.best) <= n.width):
"""Returns a sensible tolerance test.
Where between, below and beyond are not given a tolerance test function,
this determines one given an example value kindred to the final output
(plus optional offset) relative to which the tolerance is computed. """
try:
w = eg.width
try:
count = 10
while count > 0: eg, count = eg.best, count-1
while eg.width > 0 * eg.best: eg = eg.best
except AttributeError: pass
if w > 0 * eg: return bywidth
except AttributeError: pass
return microclose
# __sum is also del'd imminently
def __sum(row, plus=lambda a, b: a+b): return reduce(plus, row)
def __between(self, start, gap, edge, test, offset,
blur=__blur, gettest=__gettest, sum=__sum):
now = edge * gap # first crude estimate
# get advertised default for offset:
if offset is None:
try: offset = now - now.best
except AttributeError: offset = now * 0
# get advertised default test:
if test is None: test = gettest(offset + now)
n = 1
while 1:
n, was = 3 * n, now
h = gap / n
now = (sum(map(lambda i, b=start, s=h, f=self.integrand: f(b+i*s),
range(1, n))) + edge) * h
dif = now - was
if test(dif, now + offset): return blur(now, dif)
del __sum
def __outwards(self, bound, step, test, offset,
blur=__blur, gettest=__gettest):
small = abs(self.integrand(bound)) / 1e3
while self.__probe(bound, step) < small: step = step / 7.
while self.__probe(bound, step) > small: step = step * 7
if test is None: test = gettest(small * step)
next = bound + step
total, bound = self.between(bound, next, test, offset=offset), next
if offset is None:
try: offset = total - total.best
except AttributeError: offset = 0 * total
while 1:
step = step * 3
next = bound + step
more = self.between(bound, next, test, offset=total + offset)
total, bound = total + more, next
if test(more, total + offset): return blur(total, more)
del __blur, __gettest
# need better initial step than abs(bound) ... bound could be zero.
def __step(self, bound):
try: return self.__unit
except AttributeError: pass
try: ans = abs(bound.copy(lambda x: 1)) # for Quantity()s
except (AttributeError, TypeError): ans = abs(bound)
if ans: self.__unit = ans
else: raise ValueError, \
'Integrator needs a width parameter for .before(0) or .beyond(0)'
return ans
import math
def __probe(self, base, scale,
samples=[1, math.exp(1/math.pi), math.sqrt(2.0), 2, math.e, 3, math.pi]):
"""Scale of integrand's values for inputs base + of order scale. """
return max(map(lambda x, f=self.integrand, s=scale, b=base: abs(f(b + x*s)), samples))
del math