From fe5dd3354ad2b5b48c104dca358a9c5e3119bcd1 Mon Sep 17 00:00:00 2001 From: Christoph Groth Date: Sat, 2 Dec 2017 13:23:38 +0100 Subject: initial commit --- vquad/test/__init__.py | 0 vquad/test/test_core.py | 132 ++++++++++++++++++++++++++++++++++++++++++++++ vquad/test/test_tables.py | 53 +++++++++++++++++++ 3 files changed, 185 insertions(+) create mode 100644 vquad/test/__init__.py create mode 100644 vquad/test/test_core.py create mode 100644 vquad/test/test_tables.py (limited to 'vquad/test') diff --git a/vquad/test/__init__.py b/vquad/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/vquad/test/test_core.py b/vquad/test/test_core.py new file mode 100644 index 0000000..70b0c10 --- /dev/null +++ b/vquad/test/test_core.py @@ -0,0 +1,132 @@ +# 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 numpy.testing import assert_allclose + +from .. import core + +def f0(x): + return x * np.sin(1/x) * np.sqrt(abs(1 - x)) + + +def f7(x): + return x**-0.5 + + +def f24(x): + return np.floor(np.exp(x)) + + +def f21(x): + y = 0 + for i in range(1, 4): + y += 1 / np.cosh(20**i * (x - 2 * i / 10)) + return y + + +def f63(x, alpha, beta): + return abs(x - beta) ** alpha + + +def F63(x, alpha, beta): + return (x - beta) * abs(x - beta) ** alpha / (alpha + 1) + + +def fdiv(x): + return abs(x - 0.987654321) ** -1.1 + + +def f_one_with_nan(x): + x = np.asarray(x) + result = np.ones(x.shape) + result[x == 0] = np.inf + return result + + +def test_downdate(level=3): + vals = np.abs(core.tbls.nodes[level]) + vals[1::2] = np.nan + c_downdated = core._calc_coeffs(vals, level) + + level -= 1 + vals = np.abs(core.tbls.nodes[level]) + c = core._calc_coeffs(vals, level) + + assert_allclose(c_downdated[:len(c)], c, rtol=0, atol=1e-9) + + +def test_integration(): + old_settings = np.seterr(all='ignore') + + igral, err, nr_points = core.vquad(f0, 0, 3, 1e-5) + assert_allclose(igral, 1.98194117954329, 1e-15) + assert_allclose(err, 1.9563545589988155e-05, 1e-10) + assert nr_points == 1129 + + igral, err, nr_points = core.vquad(f7, 0, 1, 1e-6) + assert_allclose(igral, 1.9999998579359648, 1e-15) + assert_allclose(err, 1.8561437334964041e-06, 1e-10) + assert nr_points == 693 + + igral, err, nr_points = core.vquad(f24, 0, 3, 1e-3) + assert_allclose(igral, 17.664696186312934, 1e-15) + assert_allclose(err, 0.017602618074957457, 1e-10) + assert nr_points == 4519 + + igral, err, nr_points = core.vquad(f21, 0, 1, 1e-3) + assert_allclose(igral, 0.16310022131213361, 1e-15) + assert_allclose(err, 0.00011848806384952786, 1e-10) + assert nr_points == 191 + + igral, err, nr_points = core.vquad(f_one_with_nan, -1, 1, 1e-12) + assert_allclose(igral, 2, 1e-15) + assert_allclose(err, 2.4237853822937613e-15, 1e-7) + assert nr_points == 33 + + try: + igral, err, nr_points = core.vquad(fdiv, 0, 1, 1e-6) + except core.DivergentIntegralError as e: + assert e.igral == np.inf + assert e.err is None + assert e.nr_points == 431 + + np.seterr(**old_settings) + + +def test_analytic(n=200): + def f(x): + return f63(x, alpha, beta) + + def F(x): + return F63(x, alpha, beta) + + old_settings = np.seterr(all='ignore') + + np.random.seed(123) + params = np.empty((n, 2)) + params[:, 0] = np.linspace(-0.5, -1.5, n) + params[:, 1] = np.random.random_sample(n) + + false_negatives = 0 + false_positives = 0 + + for alpha, beta in params: + try: + igral, err, nr_points = core.vquad(f, 0, 1, 1e-3) + except core.DivergentIntegralError: + assert alpha < -0.8 + false_negatives += alpha > -1 + else: + if alpha <= -1: + false_positives += 1 + else: + igral_exact = F(1) - F(0) + assert alpha < -0.7 or abs(igral - igral_exact) < err + + assert false_negatives < 0.05 * n + assert false_positives < 0.05 * n + + np.seterr(**old_settings) diff --git a/vquad/test/test_tables.py b/vquad/test/test_tables.py new file mode 100644 index 0000000..79b8ffe --- /dev/null +++ b/vquad/test/test_tables.py @@ -0,0 +1,53 @@ +# 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. + +from fractions import Fraction as Frac +from itertools import combinations +import numpy as np +from numpy.testing import assert_allclose + +from .. import tables + + +def test_legendre(): + legs = tables.legendre(11) + comparisons = [(legs[0], [1], 1), + (legs[1], [0, 1], 1), + (legs[10], [-63, 0, 3465, 0, -30030, 0, + 90090, 0, -109395, 0, 46189], 256)] + for a, b, div in comparisons: + for c, d in zip(a, b): + assert c * div == d + + +def test_scalar_product(n=33): + legs = tables.legendre(n) + selection = [0, 5, 7, n-1] + for i in selection: + for j in selection: + assert (tables.scalar_product(legs[i], legs[j]) + == ((i == j) and Frac(2, 2*i + 1))) + + +def simple_newton(n): + """Slower than 'newton()' and prone to numerical error.""" + nodes = -np.cos(np.arange(n) / (n-1) * np.pi) + return [sum(np.prod(-np.asarray(sel)) + for sel in combinations(nodes, n - d)) + for d in range(n + 1)] + + +def test_newton(): + assert_allclose(tables.newton(9), simple_newton(9), atol=1e-15) + + +def test_newton_legendre(level=1): + legs = [np.array(leg, float) + for leg in tables.legendre(tables.sizes[level] + 1)] + result = np.zeros(len(legs[-1])) + for factor, leg in zip(tables.newton_coeffs[level], legs): + factor *= np.sqrt((2*len(leg) - 1) / 2) + result[:len(leg)] += factor * leg + assert_allclose(result, tables.newton(tables.sizes[level]), rtol=1e-15) -- cgit v1.2.3-74-g4815