diff options
Diffstat (limited to 'vquad/core.py')
| -rw-r--r-- | vquad/core.py | 240 |
1 files changed, 240 insertions, 0 deletions
diff --git a/vquad/core.py b/vquad/core.py new file mode 100644 index 0000000..1866276 --- /dev/null +++ b/vquad/core.py @@ -0,0 +1,240 @@ +# Copyright 2017 Christoph Groth (CEA). +# +# This file is part of Vquad. It is subject to the license terms in the file +# LICENSE.rst found in the top-level directory of this distribution. + +import numpy as np +from scipy.linalg import norm + +from . import tables as tbls + +eps = np.spacing(1) + +# If the relative difference between two consecutive approximations is +# lower than this value, the error estimate is considered reliable. +# See section 6.2 of Pedro Gonnet's thesis. +hint = 0.1 + +# Smallest acceptable relative difference of points in a rule. This was chosen +# such that no artifacts are apparent in plots of (i, log(a_i)), where a_i is +# the sequence of estimates of the integral value of an interval and all its +# ancestors.. +min_sep = 16 * eps + +min_level = 1 +max_level = 4 + +ndiv_max = 20 +max_ivals = 200 + +_sqrt_one_half = np.sqrt(0.5) + + +def _calc_coeffs(vals, level): + nans = np.flatnonzero(~np.isfinite(vals)) + if nans.size: + # Replace vals by a copy and zero-out non-finite elements. + vals = vals.copy() + vals[nans] = 0 + # Prepare things for the loop further down. + b = tbls.newton_coeffs[level].copy() + m = len(b) - 2 # = len(tbls.nodes[level]) - 1 + coeffs = tbls.inv_Vs[level] @ vals + + # This is a variant of Algorithm 7 from the thesis of Pedro Gonnet where no + # linear system has to be solved explicitly. Instead, Algorithm 5 is used. + for i in nans: + b[m + 1] /= tbls.alpha[m] + x = tbls.nodes[level][i] + b[m] = (b[m] + x * b[m + 1]) / tbls.alpha[m - 1] + for j in range(m - 1, 0, -1): + b[j] = ((b[j] + x * b[j + 1] - tbls.gamma[j + 1] * b[j + 2]) + / tbls.alpha[j - 1]) + b = b[1:] + + coeffs[:m] -= coeffs[m] / b[m] * b[:m] + coeffs[m] = 0 + m -= 1 + + return coeffs + + +class DivergentIntegralError(ValueError): + def __init__(self, msg, igral, err, nr_points): + self.igral = igral + self.err = err + self.nr_points = nr_points + super().__init__(msg) + + +class _Interval: + __slots__ = ['a', 'b', 'coeffs', 'vals', 'igral', 'err', 'level', 'depth', + 'ndiv', 'c00', 'unreliable_err'] + + def __init__(self, a, b, level, depth): + self.a = a + self.b = b + self.level = level + self.depth = depth + + def points(self): + a = self.a + b = self.b + return (a + b) / 2 + (b - a) * tbls.nodes[self.level] / 2 + + def interpolate(self, vals, coeffs_old=None): + self.vals = vals + self.coeffs = coeffs = _calc_coeffs(self.vals, self.level) + if self.level == min_level: + self.c00 = coeffs[0] + if coeffs_old is None: + coeffs_diff = norm(coeffs) + else: + coeffs_diff = np.zeros(max(len(coeffs_old), len(coeffs))) + coeffs_diff[:len(coeffs_old)] = coeffs_old + coeffs_diff[:len(coeffs)] -= coeffs + coeffs_diff = norm(coeffs_diff) + w = self.b - self.a + self.igral = w * coeffs[0] * _sqrt_one_half + self.err = w * coeffs_diff + self.unreliable_err = coeffs_diff > hint * norm(coeffs) + + +class Vquad: + """Evaluate an integral using adaptive quadrature. + + The algorithm uses Clenshaw-Curtis quadrature rules of increasing + degree in each interval. The error estimate is + sqrt(integrate((f0(x) - f1(x))**2)), where f0 and f1 are two + successive interpolations of the integrand. To fall below the + desired total error, intervals are worked on ranked by their own + absolute error: either the degree of the rule is increased or the + interval is split if either the function does not appear to be + smooth or a rule of maximum degree has been reached. + + Reference: "Increasing the Reliability of Adaptive Quadrature + Using Explicit Interpolants", P. Gonnet, ACM Transactions on + Mathematical Software, 37 (3), art. no. 26, 2008. + """ + + def __init__(self, f, a, b, level=max_level - 1): + ival = _Interval(a, b, level, 1) + vals = f(ival.points()) + ival.interpolate(vals) + ival.c00 = 0.0 # Will go away. + ival.ndiv = 0 + + self.ivals = [ival] + self.f = f + self.nr_points = len(vals) + self.igral_excess = 0 + self.err_excess = 0 + self.i_max = 0 + + def split(self, ival): + m = (ival.a + ival.b) / 2 + f_center = ival.vals[(len(ival.vals) - 1) // 2] + + depth = ival.depth + 1 + children = [_Interval(ival.a, m, min_level, depth), + _Interval(m, ival.b, min_level, depth)] + points = np.concatenate([child.points()[1:-1] for child in children]) + self.nr_points += len(points) + valss = np.empty((2, tbls.sizes[min_level])) + valss[:, 0] = ival.vals[0], f_center + valss[:, -1] = f_center, ival.vals[-1] + valss[:, 1:-1] = self.f(points).reshape((2, -1)) + + for child, vals, T in zip(children, valss, tbls.Ts): + child.interpolate(vals, T[:, :ival.coeffs.shape[0]] @ ival.coeffs) + child.ndiv = (ival.ndiv + + (ival.c00 and child.c00 / ival.c00 > 2)) + if child.ndiv > ndiv_max and 2*child.ndiv > child.depth: + msg = ('Possibly divergent integral in the interval ' + '[{}, {}]! (h={})') + raise DivergentIntegralError( + msg.format(child.a, child.b, child.b - child.a), + child.igral * np.inf, None, self.nr_points) + return children + + def refine(self, ival): + """Increase degree of interval.""" + ival.level += 1 + points = ival.points() + vals = np.empty(points.shape) + vals[0::2] = ival.vals + vals[1::2] = self.f(points[1::2]) + self.nr_points += (len(vals) - 1) // 2 + ival.interpolate(vals, ival.coeffs) + return points + + def improve(self): + i_max = self.i_max + + if self.ivals[i_max].level == max_level: + split = True + else: + points = self.refine(self.ivals[i_max]) + split = self.ivals[i_max].unreliable_err + + if (points[1] - points[0] < points[0] * min_sep + or points[-1] - points[-2] < points[-2] * min_sep + or (self.ivals[i_max].err + < (abs(self.ivals[i_max].igral) * eps + * tbls.V_cond_nums[self.ivals[i_max].level]))): + # Remove the interval (while remembering the excess integral + # and error), since it is either too narrow, or the estimated + # relative error is already at the limit of numerical accuracy + # and cannot be reduced further. + self.err_excess += self.ivals[i_max].err + self.igral_excess += self.ivals[i_max].igral + self.ivals[i_max] = self.ivals[-1] + self.ivals.pop() + return + + if split: + self.ivals.extend(self.split(self.ivals[i_max])) + self.ivals[i_max] = self.ivals.pop() + + def totals(self): + # Compute the total error and new max. + i_max = 0 + i_min = 0 + err = self.err_excess + igral = self.igral_excess + for i in range(len(self.ivals)): + if self.ivals[i].err > self.ivals[i_max].err: + i_max = i + elif self.ivals[i].err < self.ivals[i_min].err: + i_min = i + err += self.ivals[i].err + igral += self.ivals[i].igral + + # If there are too many intervals, remove the one with smallest + # contribution to the error. + if len(self.ivals) > max_ivals: + self.err_excess += self.ivals[i_min].err + self.igral_excess += self.ivals[i_min].igral + self.ivals[i_min] = self.ivals[-1] + self.ivals.pop() + if i_max == len(self.ivals): + i_max = i_min + + self.i_max = i_max + return igral, err + + def improve_until(self, tol): + while True: + self.improve() + igral, err = self.totals() + + if (err == 0 + or err < abs(igral) * tol + or err - self.err_excess < abs(igral) * tol < self.err_excess + or not self.ivals): + return igral, err, self.nr_points + + +def vquad(f, a, b, tol): + igrator = Vquad(f, a, b) + return igrator.improve_until(tol) |