# Copyright (c) 2017, The Chancellor, Masters and Scholars of the University # of Oxford, and the Chebfun Developers. All rights reserved. # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are met: # * Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # * Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the distribution. # * Neither the name of the University of Oxford nor the names of its # contributors may be used to endorse or promote products derived from # this software without specific prior written permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND # ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED # WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR # ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES # (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; # LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND # ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. from math import factorial import numpy as np from numpy.testing import assert_allclose, assert_equal, assert_array_less import pytest import scipy from scipy.interpolate import AAA, FloaterHormannInterpolator, BarycentricInterpolator TOL = 1e4 * np.finfo(np.float64).eps UNIT_INTERVAL = np.linspace(-1, 1, num=1000) PTS = np.logspace(-15, 0, base=10, num=500) PTS = np.concatenate([-PTS[::-1], [0], PTS]) @pytest.mark.parametrize("method", [AAA, FloaterHormannInterpolator]) @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128]) def test_dtype_preservation(method, dtype): rtol = np.finfo(dtype).eps ** 0.75 * 100 if method is FloaterHormannInterpolator: rtol *= 100 rng = np.random.default_rng(59846294526092468) z = np.linspace(-1, 1, dtype=dtype) r = method(z, np.sin(z)) z2 = rng.uniform(-1, 1, size=100).astype(dtype) assert_allclose(r(z2), np.sin(z2), rtol=rtol) assert r(z2).dtype == dtype if method is AAA: assert r.support_points.dtype == dtype assert r.support_values.dtype == dtype assert r.errors.dtype == z.real.dtype assert r.weights.dtype == dtype assert r.poles().dtype == np.result_type(dtype, 1j) assert r.residues().dtype == np.result_type(dtype, 1j) assert r.roots().dtype == np.result_type(dtype, 1j) @pytest.mark.parametrize("method", [AAA, FloaterHormannInterpolator]) @pytest.mark.parametrize("dtype", [np.int16, np.int32, np.int64]) def test_integer_promotion(method, dtype): z = np.arange(10, dtype=dtype) r = method(z, z) assert r.weights.dtype == np.result_type(dtype, 1.0) if method is AAA: assert r.support_points.dtype == np.result_type(dtype, 1.0) assert r.support_values.dtype == np.result_type(dtype, 1.0) assert r.errors.dtype == np.result_type(dtype, 1.0) assert r.poles().dtype == np.result_type(dtype, 1j) assert r.residues().dtype == np.result_type(dtype, 1j) assert r.roots().dtype == np.result_type(dtype, 1j) assert r(z).dtype == np.result_type(dtype, 1.0) class TestAAA: def test_input_validation(self): with pytest.raises(ValueError, match="same size"): AAA([0], [1, 1]) with pytest.raises(ValueError, match="1-D"): AAA([[0], [0]], [[1], [1]]) with pytest.raises(ValueError, match="finite"): AAA([np.inf], [1]) with pytest.raises(TypeError): AAA([1], [1], max_terms=1.0) with pytest.raises(ValueError, match="greater"): AAA([1], [1], max_terms=-1) @pytest.mark.thread_unsafe def test_convergence_error(self): with pytest.warns(RuntimeWarning, match="AAA failed"): AAA(UNIT_INTERVAL, np.exp(UNIT_INTERVAL), max_terms=1) # The following tests are based on: # https://github.com/chebfun/chebfun/blob/master/tests/chebfun/test_aaa.m def test_exp(self): f = np.exp(UNIT_INTERVAL) r = AAA(UNIT_INTERVAL, f) assert_allclose(r(UNIT_INTERVAL), f, atol=TOL) assert_equal(r(np.nan), np.nan) assert np.isfinite(r(np.inf)) m1 = r.support_points.size r = AAA(UNIT_INTERVAL, f, rtol=1e-3) assert r.support_points.size < m1 def test_tan(self): f = np.tan(np.pi * UNIT_INTERVAL) r = AAA(UNIT_INTERVAL, f) assert_allclose(r(UNIT_INTERVAL), f, atol=10 * TOL, rtol=1.4e-7) assert_allclose(np.min(np.abs(r.roots())), 0, atol=3e-10) assert_allclose(np.min(np.abs(r.poles() - 0.5)), 0, atol=TOL) # Test for spurious poles (poles with tiny residue are likely spurious) assert np.min(np.abs(r.residues())) > 1e-13 def test_short_cases(self): # Computed using Chebfun: # >> format long # >> [r, pol, res, zer, zj, fj, wj, errvec] = aaa([1 2], [0 1]) z = np.array([0, 1]) f = np.array([1, 2]) r = AAA(z, f, rtol=1e-13) assert_allclose(r(z), f, atol=TOL) assert_allclose(r.poles(), 0.5) assert_allclose(r.residues(), 0.25) assert_allclose(r.roots(), 1/3) assert_equal(r.support_points, z) assert_equal(r.support_values, f) assert_allclose(r.weights, [0.707106781186547, 0.707106781186547]) assert_equal(r.errors, [1, 0]) # >> format long # >> [r, pol, res, zer, zj, fj, wj, errvec] = aaa([1 0 0], [0 1 2]) z = np.array([0, 1, 2]) f = np.array([1, 0, 0]) r = AAA(z, f, rtol=1e-13) assert_allclose(r(z), f, atol=TOL) assert_allclose(np.sort(r.poles()), np.sort([1.577350269189626, 0.422649730810374])) assert_allclose(np.sort(r.residues()), np.sort([-0.070441621801729, -0.262891711531604])) assert_allclose(np.sort(r.roots()), np.sort([2, 1])) assert_equal(r.support_points, z) assert_equal(r.support_values, f) assert_allclose(r.weights, [0.577350269189626, 0.577350269189626, 0.577350269189626]) assert_equal(r.errors, [1, 1, 0]) def test_scale_invariance(self): z = np.linspace(0.3, 1.5) f = np.exp(z) / (1 + 1j) r1 = AAA(z, f) r2 = AAA(z, (2**311 * f).astype(np.complex128)) r3 = AAA(z, (2**-311 * f).astype(np.complex128)) assert_equal(r1(0.2j), 2**-311 * r2(0.2j)) assert_equal(r1(1.4), 2**311 * r3(1.4)) def test_log_func(self): rng = np.random.default_rng(1749382759832758297) z = rng.standard_normal(10000) + 3j * rng.standard_normal(10000) def f(z): return np.log(5 - z) / (1 + z**2) r = AAA(z, f(z)) assert_allclose(r(0), f(0), atol=TOL) def test_infinite_data(self): z = np.linspace(-1, 1) r = AAA(z, scipy.special.gamma(z)) assert_allclose(r(0.63), scipy.special.gamma(0.63), atol=1e-15) def test_nan(self): x = np.linspace(0, 20) with np.errstate(invalid="ignore"): f = np.sin(x) / x r = AAA(x, f) assert_allclose(r(2), np.sin(2) / 2, atol=1e-15) def test_residues(self): x = np.linspace(-1.337, 2, num=537) r = AAA(x, np.exp(x) / x) ii = np.flatnonzero(np.abs(r.poles()) < 1e-8) assert_allclose(r.residues()[ii], 1, atol=1e-15) r = AAA(x, (1 + 1j) * scipy.special.gamma(x)) ii = np.flatnonzero(abs(r.poles() - (-1)) < 1e-8) assert_allclose(r.residues()[ii], -1 - 1j, atol=1e-15) # The following tests are based on: # https://github.com/complexvariables/RationalFunctionApproximation.jl/blob/main/test/interval.jl @pytest.mark.parametrize("func,atol,rtol", [(lambda x: np.abs(x + 0.5 + 0.01j), 5e-13, 1e-7), (lambda x: np.sin(1/(1.05 - x)), 2e-13, 1e-7), (lambda x: np.exp(-1/(x**2)), 3.5e-13, 0), (lambda x: np.exp(-100*x**2), 8e-13, 0), (lambda x: np.exp(-10/(1.2 - x)), 1e-14, 0), (lambda x: 1/(1+np.exp(100*(x + 0.5))), 2e-13, 1e-7), (lambda x: np.abs(x - 0.95), 1e-6, 1e-7)]) def test_basic_functions(self, func, atol, rtol): with np.errstate(divide="ignore"): f = func(PTS) assert_allclose(AAA(UNIT_INTERVAL, func(UNIT_INTERVAL))(PTS), f, atol=atol, rtol=rtol) def test_poles_zeros_residues(self): def f(z): return (z+1) * (z+2) / ((z+3) * (z+4)) r = AAA(UNIT_INTERVAL, f(UNIT_INTERVAL)) assert_allclose(np.sum(r.poles() + r.roots()), -10, atol=1e-12) def f(z): return 2/(3 + z) + 5/(z - 2j) r = AAA(UNIT_INTERVAL, f(UNIT_INTERVAL)) assert_allclose(r.residues().prod(), 10, atol=1e-8) r = AAA(UNIT_INTERVAL, np.sin(10*np.pi*UNIT_INTERVAL)) assert_allclose(np.sort(np.abs(r.roots()))[18], 0.9, atol=1e-12) def f(z): return (z - (3 + 3j))/(z + 2) r = AAA(UNIT_INTERVAL, f(UNIT_INTERVAL)) assert_allclose(r.poles()[0]*r.roots()[0], -6-6j, atol=1e-12) @pytest.mark.parametrize("func", [lambda z: np.zeros_like(z), lambda z: z, lambda z: 1j*z, lambda z: z**2 + z, lambda z: z**3 + z, lambda z: 1/(1.1 + z), lambda z: 1/(1 + 1j*z), lambda z: 1/(3 + z + z**2), lambda z: 1/(1.01 + z**3)]) def test_polynomials_and_reciprocals(self, func): assert_allclose(AAA(UNIT_INTERVAL, func(UNIT_INTERVAL))(PTS), func(PTS), atol=2e-13) # The following tests are taken from: # https://github.com/macd/BaryRational.jl/blob/main/test/test_aaa.jl def test_spiral(self): z = np.exp(np.linspace(-0.5, 0.5 + 15j*np.pi, num=1000)) r = AAA(z, np.tan(np.pi*z/2)) assert_allclose(np.sort(np.abs(r.poles()))[:4], [1, 1, 3, 3], rtol=9e-7) @pytest.mark.thread_unsafe def test_spiral_cleanup(self): z = np.exp(np.linspace(-0.5, 0.5 + 15j*np.pi, num=1000)) # here we set `rtol=0` to force froissart doublets, without cleanup there # are many spurious poles with pytest.warns(RuntimeWarning): r = AAA(z, np.tan(np.pi*z/2), rtol=0, max_terms=60, clean_up=False) n_spurious = np.sum(np.abs(r.residues()) < 1e-14) with pytest.warns(RuntimeWarning): assert r.clean_up() >= 1 # check there are less potentially spurious poles than before assert np.sum(np.abs(r.residues()) < 1e-14) < n_spurious # check accuracy assert_allclose(r(z), np.tan(np.pi*z/2), atol=6e-12, rtol=3e-12) class TestFloaterHormann: def runge(self, z): return 1/(1 + z**2) def scale(self, n, d): return (-1)**(np.arange(n) + d) * factorial(d) def test_iv(self): with pytest.raises(ValueError, match="`x`"): FloaterHormannInterpolator([[0]], [0], d=0) with pytest.raises(ValueError, match="`y`"): FloaterHormannInterpolator([0], 0, d=0) with pytest.raises(ValueError, match="dimension"): FloaterHormannInterpolator([0], [[1, 1], [1, 1]], d=0) with pytest.raises(ValueError, match="finite"): FloaterHormannInterpolator([np.inf], [1], d=0) with pytest.raises(ValueError, match="`d`"): FloaterHormannInterpolator([0], [0], d=-1) with pytest.raises(ValueError, match="`d`"): FloaterHormannInterpolator([0], [0], d=10) with pytest.raises(TypeError): FloaterHormannInterpolator([0], [0], d=0.0) # reference values from Floater and Hormann 2007 page 8. @pytest.mark.parametrize("d,expected", [ (0, [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]), (1, [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1]), (2, [1, 3, 4, 4, 4, 4, 4, 4, 4, 3, 1]), (3, [1, 4, 7, 8, 8, 8, 8, 8, 7, 4, 1]), (4, [1, 5, 11, 15, 16, 16, 16, 15, 11, 5, 1]) ]) def test_uniform_grid(self, d, expected): # Check against explicit results on an uniform grid x = np.arange(11) r = FloaterHormannInterpolator(x, 0.0*x, d=d) assert_allclose(r.weights.ravel()*self.scale(x.size, d), expected, rtol=1e-15, atol=1e-15) @pytest.mark.parametrize("d", range(10)) def test_runge(self, d): x = np.linspace(0, 1, 51) rng = np.random.default_rng(802754237598370893) xx = rng.uniform(0, 1, size=1000) y = self.runge(x) h = x[1] - x[0] r = FloaterHormannInterpolator(x, y, d=d) tol = 10*h**(d+1) assert_allclose(r(xx), self.runge(xx), atol=1e-10, rtol=tol) # check interpolation property assert_equal(r(x), self.runge(x)) def test_complex(self): x = np.linspace(-1, 1) z = x + x*1j r = FloaterHormannInterpolator(z, np.sin(z), d=12) xx = np.linspace(-1, 1, num=1000) zz = xx + xx*1j assert_allclose(r(zz), np.sin(zz), rtol=1e-12) def test_polyinterp(self): # check that when d=n-1 FH gives a polynomial interpolant x = np.linspace(0, 1, 11) xx = np.linspace(0, 1, 1001) y = np.sin(x) r = FloaterHormannInterpolator(x, y, d=x.size-1) p = BarycentricInterpolator(x, y) assert_allclose(r(xx), p(xx), rtol=1e-12, atol=1e-12) @pytest.mark.parametrize("y_shape", [(2,), (2, 3, 1), (1, 5, 6, 4)]) @pytest.mark.parametrize("xx_shape", [(100), (10, 10)]) def test_trailing_dim(self, y_shape, xx_shape): x = np.linspace(0, 1) y = np.broadcast_to( np.expand_dims(np.sin(x), tuple(range(1, len(y_shape) + 1))), x.shape + y_shape ) r = FloaterHormannInterpolator(x, y) rng = np.random.default_rng(897138947238097528091759187597) xx = rng.random(xx_shape) yy = np.broadcast_to( np.expand_dims(np.sin(xx), tuple(range(xx.ndim, len(y_shape) + xx.ndim))), xx.shape + y_shape ) rr = r(xx) assert rr.shape == xx.shape + y_shape assert_allclose(rr, yy, rtol=1e-6) def test_zeros(self): x = np.linspace(0, 10, num=100) r = FloaterHormannInterpolator(x, np.sin(np.pi*x)) err = np.abs(np.subtract.outer(r.roots(), np.arange(11))).min(axis=0) assert_array_less(err, 1e-5) def test_no_poles(self): x = np.linspace(-1, 1) r = FloaterHormannInterpolator(x, 1/x**2) p = r.poles() mask = (p.real >= -1) & (p.real <= 1) & (np.abs(p.imag) < 1.e-12) assert np.sum(mask) == 0