From 2c5446877794b7a900bb97f084f837253d78417a Mon Sep 17 00:00:00 2001
From: Christoph Groth <christoph.groth@cea.fr>
Date: Tue, 5 Jun 2018 11:55:07 +0200
Subject: Lyness-Kaganove benchmark

---
 README.rst         |  43 +++++++-
 vquad/__main__.py  |  19 ++++
 vquad/benchmark.py | 289 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 349 insertions(+), 2 deletions(-)
 create mode 100644 vquad/__main__.py
 create mode 100644 vquad/benchmark.py

diff --git a/README.rst b/README.rst
index 12eba6d..838add1 100644
--- a/README.rst
+++ b/README.rst
@@ -1,9 +1,48 @@
 Vquad
 =====
 
-Vquad is a work-in-progress towards a general-purpose, robust, efficient, and parallel library for numerical integration.
+Vquad is a work-in-progress towards a general-purpose, robust, efficient, and
+parallel library for numerical integration.
 
 History
 -------
 
-Vquad is based on algorithm 4 from the article "Increasing the Reliability of Adaptive Quadrature Using Explicit Interpolants", Pedro Gonnet, ACM TOMS 37, 26 (2010).
+Vquad is based on algorithm 4 from the article "Increasing the Reliability of
+Adaptive Quadrature Using Explicit Interpolants", Pedro Gonnet, ACM TOMS 37,
+26 (2010).
+
+Usage example
+-------------
+
+```
+import numpy
+import vquad
+
+# Simple usage.
+igral, err = vquad.vquad(numpy.cos, -1, 1, rtol=1e-10)
+
+# Object interface.
+it = vquad.Vquad(numpy.cos, -1, 1)
+igral, err = it.improve_until(rtol=1e-10)
+
+# Evaluate interpolant.
+xs = numpy.linspace(-1, 1, 101)
+ys = it(xs)
+```
+
+Benachmarks and tests
+---------------------
+
+Vquad includes extensive tests and benchmarks.  The tests focus on verifying
+correctness and can be run with
+```
+python3 -m pytest -s
+```
+The `-s` shows some statistics that are gathered during the testing.
+
+The benchmarks take longer to run and allow to quantitatively compare the
+performance of the algorithm.  The benchmark is run by executing the vquad
+module.  To get help, run:
+```
+python3 -m vquad -h
+```
diff --git a/vquad/__main__.py b/vquad/__main__.py
new file mode 100644
index 0000000..280de87
--- /dev/null
+++ b/vquad/__main__.py
@@ -0,0 +1,19 @@
+# Copyright 2018 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 sys
+import traceback
+from . import benchmark
+
+# Make sure that whenever there is an exception it is printed and an
+# appropriate exit code is set.
+rc = 1
+try:
+    benchmark.main()
+    rc = 0
+except (Exception, KeyboardInterrupt):
+    traceback.print_exc()
+finally:
+    sys.exit(rc)
diff --git a/vquad/benchmark.py b/vquad/benchmark.py
new file mode 100644
index 0000000..1feb3d8
--- /dev/null
+++ b/vquad/benchmark.py
@@ -0,0 +1,289 @@
+# Copyright 2018 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 math
+import itertools
+import argparse
+import pickle
+from concurrent.futures import ProcessPoolExecutor
+
+import numpy as np
+import scipy.integrate
+
+import matplotlib, matplotlib.cm
+from matplotlib import pyplot
+
+from . import core
+
+
+def f63(x, α, λ):
+    return np.abs(x - λ)**α
+
+def F63_0_1(α, λ):
+    def F(x):
+        return (x - λ) * abs(x - λ)**α / (α + 1)
+    return F(1) - F(0)
+
+
+def f64(x, α, λ):
+    return (x > λ) * np.exp(α * x)
+
+def F64_0_1(α, λ):
+    assert 0 <= λ <= 1
+    if α == 0:
+        return 1 - λ
+    return (math.exp(α) - math.exp(α * λ)) / α
+
+
+def f65(x, α, λ):
+    return np.exp(-α * np.abs(x - λ))
+
+def F65_0_1(α, λ):
+    assert 0 <= λ <= 1
+    if α == 0:
+        return 1
+    return (math.expm1(α*λ - α) + math.expm1(-α*λ)) / -α
+
+
+# As shown in Pedro's thesis (Eq. (6.6) and Eq. (6.7)), the term 10^α is not
+# squared in the denominator.  However, this is the correct version of the
+# integrand that is consistent with the benchmark results that are shown in the
+# thesis.
+def f67(x, α, *λs):
+    t = 10**α
+    return t * sum(1 / ((x - λ)**2 + t**2) for λ in λs)
+
+def F67_1_2(α, *λs):
+    t = 10**α
+    return sum(math.atan((2 - λ) / t) - math.atan((1 - λ) / t) for λ in λs)
+
+
+# Equation (6.8) is the only one where the benchmark results do not agree with
+# the ones listed in Pedro Gonnet's thesis.  So Eq. (6.8) seems to contain a
+# typo.  I have modified the function by adding a small constant so that the
+# definite integral is never zero.  Otherwise any required relative error
+# cannot be obtained which leads to extremely long running times.
+def f68(x, α, λ):
+    β = 10**α / max(λ**2, (1 - λ)**2)
+    return 2 + 2 * β * (x - λ) * np.cos(β * (x - λ)**2)
+
+def F68_0_1(α, λ):
+    β = 10**α / max(λ**2, (1 - λ)**2)
+    return 2 + math.sin(β * ((1 - λ)**2)) - math.sin(β * λ**2)
+
+
+families = [
+    (r"|x - \lambda|^\alpha",
+     (f63, 0, 1), F63_0_1, (-0.5, 0), (0, 1)),
+
+    (r"(x - \lambda)\mathrm{e}^{\alpha x}",
+     (f64, 0, 1), F64_0_1, (0, 1), (0, 1)),
+
+    (r"\exp(-\alpha|x - \lambda|)",
+     (f65, 0, 1), F65_0_1, (0, 4), (0, 1)),
+
+    (r"10^\alpha/((x - \lambda)^2 + 10^{2\alpha})",
+     (f67, 1, 2), F67_1_2, (-6, -3), (1, 2)),
+
+    (r"\sum_i 10^\alpha/((x - \lambda_i)^2 + 10^{2\alpha})",
+     (f67, 1, 2), F67_1_2, (-5, -3)) + ((1, 2),) * 4,
+
+    (r"2 + 2\beta(x - \lambda) \cos(\beta (x - \lambda)^2),\quad"
+     r"\beta = 10^\alpha / max(\lambda^2, (1 - \lambda)^2)",
+     (f68, 0, 1), F68_0_1, (1.8, 2), (0, 1)),
+]
+
+
+def lk_job_vquad(seed, family, rtols):
+    """Run a single integration for a sequence of requested relative
+       tolerances.  It's as if a len(rtols) integration were started,
+       but actually a single integrator is launched and the requested
+       tolerance is improved iteratively.
+    """
+    def ff(x):
+        nonlocal neval
+        neval += len(x)
+        return f(x, *args)
+
+    _, (f, a, b), F, *bounds = family
+    rng = np.random.RandomState(seed)
+    args = [rng.uniform(*b) for b in bounds]
+
+    relerrs = np.empty(len(rtols), float)
+    nevals = np.empty(len(rtols), int)
+
+    neval = 0
+    exact = F(*args)
+    vquad = core.Vquad(ff, a, b)
+    for j, rtol in enumerate(rtols):
+        igral, _ = vquad.improve_until(rtol)
+        relerrs[j] = abs((igral - exact) / exact)
+        nevals[j] = neval
+
+    return nevals, relerrs
+
+
+def lk_job_quadpack(seed, family, rtols):
+    """Run a single integration for a sequence of requested relative
+       tolerances.  It's as if a len(rtols) integration were started,
+       but actually a single integrator is launched and the requested
+       tolerance is improved iteratively.
+    """
+    def ff(x):
+        nonlocal neval
+        neval += 1
+        return f(x, *args)
+
+    _, (f, a, b), F, *bounds = family
+    rng = np.random.RandomState(seed)
+    args = [rng.uniform(*b) for b in bounds]
+
+    relerrs = np.empty(len(rtols), float)
+    nevals = np.empty(len(rtols), int)
+
+    exact = F(*args)
+    for j, rtol in enumerate(rtols):
+        neval = 0
+        igral, _ = scipy.integrate.quad(
+                    ff, a, b, epsabs=0, epsrel=rtol)
+        relerrs[j] = abs((igral - exact) / exact)
+        nevals[j] = neval
+
+    return nevals, relerrs
+
+
+def lk_runner(map, n, job, *args):
+    """Run job n times with different seed but otherwise same args.
+
+    return two arrays of shape (n, -1) where the second dimension is
+    determined by the args.
+    """
+    # map returns a sequence of (nevals, relerrs) pairs that is unzipped by the
+    # zip(*...) construct.
+    nevals, relerrs = zip(*map(job, range(n),
+                               *(itertools.repeat(a) for a in args)))
+    return np.stack(nevals), np.stack(relerrs)
+
+
+def run(runner, job, n, rtols, thresholds, ):
+    old_settings = np.seterr(all='ignore')
+
+    executor = ProcessPoolExecutor()
+    for family in families:
+        nevals, relerrs = runner(executor.map, n, job, family, rtols)
+        nevals = np.mean(nevals, axis=0)
+        relerrs.sort(axis=0)
+        successratios = [np.searchsorted(*args) / n
+                         for args in zip(relerrs.T, rtols)]
+        quantiles = {t: relerrs[int(t * n)] for t in thresholds}
+        yield nevals, quantiles, successratios
+
+    np.seterr(**old_settings)
+
+
+def display(datasets, dataset_names, familynames, cols, colnames):
+    matplotlib.rcParams['axes.prop_cycle'] = matplotlib.cycler(
+        'color', matplotlib.cm.tab20.colors)
+
+    fig, axes = pyplot.subplots(2, 3, sharex='all', sharey='all')
+    fig.suptitle("Lyness-Kaganove diagrams")
+
+    print("Lyness-Kaganove benchmark as described in Section 6.5 "
+          "of P. Gonnet's thesis.\n"
+          "Legend: success count (mean number of evaluations)\n")
+    print("family \ rtol\t{}\t\t{}\t\t{}".format(*colnames))
+
+    for (i, familyname), ((m, n), ax), *results in zip(
+            enumerate(familynames), np.ndenumerate(axes), *datasets):
+        if i:
+            print()
+        line = [f"f_{i}\t"]
+        for j, (nevals, quantiles, successratios) in enumerate(results):
+            for col in cols:
+                line.append(f"{successratios[col]} ({nevals[col]:.0f})")
+            print("\t".join(line))
+            line = ["\t"]
+
+            for k, (t, q) in enumerate(quantiles.items()):
+                if j == m == n == 0:
+                    label = f'{100*t:.0f}% quantile'
+                elif k == m == 0 and n == 1:
+                    label = dataset_names[j]
+                else:
+                    label = None
+                ax.loglog(q, nevals, 'o-', markevery=10,
+                          label=label)
+
+        ax.set_title(f'$f_{i}(x) = {familyname}$')
+        if m == 0 and n in (0, 1):
+            ax.legend()
+        if m == 1:
+            ax.set_xlabel("Relative true error")
+        if n == 0:
+            ax.set_ylabel("Mean number of evaluations")
+
+
+def main():
+    VERSION = 'vquad benchmark dump version 0'
+
+    #### Parse command line.
+    examples = """Examples:
+python3 -m vquad -
+python3 -m vquad -a quadpack -o quadpack
+python3 -m vquad quadpack vquad-original -
+"""
+    parser = argparse.ArgumentParser(
+        'python3 -m vquad',
+        description='Benchmark the integrator.',
+        formatter_class=argparse.RawDescriptionHelpFormatter,
+        epilog=examples)
+    parser.add_argument('datasets', nargs='*')
+    parser.add_argument('-a', '--algorithm', default='vquad',
+                        choices=['vquad', 'quadpack'])
+    parser.add_argument('-n', '--samples', type=int, default=200)
+    parser.add_argument('-o', '--output')
+    args = parser.parse_args()
+
+    # run_quadpack(args.samples, 10.0**np.linspace(0, -10, 101))
+
+    #### Setup benchmark.
+    if args.output or '-' in args.datasets:
+        if args.algorithm == 'vquad':
+            job = lk_job_vquad
+        elif args.algorithm == 'quadpack':
+            job = lk_job_quadpack
+        else:
+            raise RuntimeError('Unknown algorithm.')
+        computed, tobesaved = itertools.tee(run(
+            lk_runner, job,
+            args.samples, 10.0**np.linspace(0, -10, 101), [0.5, 0.95]))
+
+    datasets = []
+    dataset_names = []
+    #### Load any datasets.
+    for fname in args.datasets:
+        if fname == "-":
+            datasets.append(computed)
+            dataset_names.append('<{}>'.format(args.algorithm))
+            continue
+        with open(fname, 'rb') as f:
+            version, loaded = pickle.load(f)
+            if version != VERSION:
+                raise RuntimeError("Invalid version of file " + fname)
+            else:
+                datasets.append(loaded)
+                dataset_names.append(fname)
+
+    #### Launch computation and display.
+    if datasets:
+        display(datasets, dataset_names, (f[0] for f in families),
+                [30, 60, 90], ['1e-3', '1e-6', '1e-9'])
+
+    #### Save benchmark.
+    if args.output:
+        with open(args.output, 'wb') as f:
+            pickle.dump((VERSION, list(tobesaved)), f)
+
+    pyplot.show()
-- 
cgit v1.2.3-74-g4815