From 1b9e8d365fc45473183af149b630980862690600 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Tue, 7 Dec 2021 10:43:38 +0100 Subject: [PATCH] ENH: Extend chebinterpolate to interpolate a multidimensional function Currently there are no functions in chebyshev.py that do N-dimensional fitting for N>1 This submission is a attempt to fill the above mentioned gap and complements the efforts in #14151. The multidimensional interpolation method implemented here uses the 1D Discrete Cosine Transform (DCT) instead of the vandermonde method. This is a reimplementation of the chebfitfunc function found in #6071. The accuracy of the new interpolation method is better than the old one especially for polynomials of high degrees. Also added domain as input to chebinterpolate in order to interpolate the function over any input domain. Added test_2d_approximation for testing interpolation of a 2D function. Replaced all references to the missing chebfromfunction with chebinterpolate. Made doctest of chebinterpolate more robust. --- numpy/polynomial/chebyshev.py | 217 ++++++++++++++++++++--- numpy/polynomial/tests/test_chebyshev.py | 14 ++ 2 files changed, 203 insertions(+), 28 deletions(-) diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py index 5c595bcf668d..90011de99ffc 100644 --- a/numpy/polynomial/chebyshev.py +++ b/numpy/polynomial/chebyshev.py @@ -1777,27 +1777,165 @@ def chebroots(c): return r -def chebinterpolate(func, deg, args=()): +def __dct(x, n=None): + """ + Discrete Cosine Transform type II + N-1 + y[k] = 2 * sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. + n=0 + Examples + -------- + >>> x = np.arange(5) + >>> x1 = np.arange(6) + >>> y = x + 1j*x + >>> y1 = x1 + 1j*x1 + + >>> np.allclose(__dct(x), [20. , -9.95959314, 0. , -0.89805595, 0. ]) + True + + >>> np.allclose(__dct(x1), + ... [30. , -14.41953704, 0. , -1.41421356, 0. , -0.27740142]) + True + + >>> np.allclose(__dct(y), [20. +20.j, -9.95959314 -9.95959314j, + ... 0. +0.j, -0.89805595 -0.89805595j, + ... 0. +0.j ]) + True + + >>> np.allclose(__dct(y1), [ 30. +30.j , -14.41953704-14.41953704j, + ... 0. +0.j , -1.41421356 -1.41421356j, + ... 0. +0.j , -0.27740142 -0.27740142j]) + True + + >>> np.allclose(x, __idct(__dct(x))) + True + + >>> np.allclose(x, __dct(__idct(x))) + True + + References + ---------- + .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J. + Makhoul, `IEEE Transactions on acoustics, speech and signal + processing` vol. 28(1), pp. 27-34, + :doi:`10.1109/TASSP.1980.1163351` (1980). + .. [2] Wikipedia, "Discrete cosine transform", + https://en.wikipedia.org/wiki/Discrete_cosine_transform + """ + fft = np.fft.fft + x = np.atleast_1d(x) + + if n is None: + n = x.shape[-1] + + if x.shape[-1] < n: + n_shape = x.shape[:-1] + (n - x.shape[-1],) + xx = np.hstack((x, np.zeros(n_shape))) + else: + xx = x[..., :n] + + real_x = np.all(np.isreal(xx)) + if (real_x and (np.remainder(n, 2) == 0)): + xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2]))) + else: + xp = fft(np.hstack((xx, xx[..., ::-1])))[..., :n] + + w = np.exp(-1j * np.arange(n) * np.pi / (2 * n)) + + y = xp * w + + if real_x: + return y.real + return y + + +def __idct(x, n=None): + """ + Inverse Discrete Cosine Transform type II + N-1 + y[k] = 1/N sum w[n]*x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. + n=0 + w[0] = 1/2 + w[n] = 1 for n>0 + + Examples + -------- + >>> x = np.arange(5) + >>> np.allclose(x, __idct(__dct(x))) + True + + >>> np.allclose(x, __dct(__idct(x))) + True + + References + ---------- + .. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J. + Makhoul, `IEEE Transactions on acoustics, speech and signal + processing` vol. 28(1), pp. 27-34, + :doi:`10.1109/TASSP.1980.1163351` (1980). + .. [2] Wikipedia, "Discrete cosine transform", + https://en.wikipedia.org/wiki/Discrete_cosine_transform + """ + + ifft = np.fft.ifft + x = np.atleast_1d(x) + + if n is None: + n = x.shape[-1] + + w = np.exp(1j * np.arange(n) * np.pi / (2 * n)) + + if x.shape[-1] < n: + n_shape = x.shape[:-1] + (n - x.shape[-1],) + xx = np.hstack((x, np.zeros(n_shape))) * w + else: + xx = x[..., :n] * w + + real_x = np.all(np.isreal(x)) + if (real_x and (np.remainder(n, 2) == 0)): + xx[..., 0] = xx[..., 0] * 0.5 + yp = ifft(xx) + y = np.zeros(xx.shape, dtype=complex) + y[..., ::2] = yp[..., :n // 2] + y[..., ::-2] = yp[..., n // 2::] + else: + yp = ifft(np.hstack((xx, + np.zeros_like(xx[..., 0]), + np.conj(xx[..., :0:-1])))) + y = yp[..., :n] + + if real_x: + return y.real + return y + + +def chebinterpolate(func, deg, args=(), domain=None): """Interpolate a function at the Chebyshev points of the first kind. - Returns the Chebyshev series that interpolates `func` at the Chebyshev - points of the first kind in the interval [-1, 1]. The interpolating - series tends to a minmax approximation to `func` with increasing `deg` - if the function is continuous in the interval. + Returns the Chebyshev series that interpolates the multi-dimensional + function `func` at the Chebyshev points of the first kind scaled and + shifted to the `domain`. The interpolating series tends to a minmax + approximation of `func` with increasing `deg` if the function is + continuous in the input domain. .. versionadded:: 1.14.0 Parameters ---------- func : function - The function to be approximated. It must be a function of a single - variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are - extra arguments passed in the `args` parameter. - deg : int - Degree of the interpolating polynomial + The function to be approximated. It must be a function of the form + ``f(x1, x2,..., xn, a, b, c...)``, where ``x1, x2,..., xn`` are the + n-dimensional coordinates and ``a, b, c...`` are extra arguments + passed in the `args` parameter. + deg : list of int, size ndim + Degree of the interpolating polynomial used for each dimension. args : tuple, optional Extra arguments to be used in the function call. Default is no extra arguments. + domain : array-like, shape ndim x 2, optional + defines the domain ``[a1, b1] x [a2, b2] x ... x [an, bn]`` over + which `func` is interpolated. The default is None, in which case + the domain is [-1,1] * len(deg)) Returns ------- @@ -1808,10 +1946,17 @@ def chebinterpolate(func, deg, args=()): Examples -------- >>> import numpy.polynomial.chebyshev as C - >>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8) - array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17, - -5.42457905e-02, -2.71387850e-16, 4.51658839e-03, - 2.46716228e-17, -3.79694221e-04, -3.26899002e-16]) + >>> c1 = C.chebinterpolate(lambda x: np.tanh(x) + 0.5, 8) + >>> x = C.chebpts1(9) + >>> np.allclose(C.chebval(x, c1), np.tanh(x)+0.5)) + True + + Two-dimensional interpolation to machine accuracy + >>> x1,x2 = np.meshgrid(x,x) + >>> c2 = C.chebinterpolate(lambda x,y: np.tanh(x+y) + 0.5, deg=(8, 8)) + >>> np.allclose(C.chebval2d(x1, x2, c2), np.tanh(x1+x2)+0.5)) + True + Notes ----- @@ -1824,24 +1969,41 @@ def chebinterpolate(func, deg, args=()): zero. For instance, if the function is even then the coefficients of the terms of odd degree in the result can be set to zero. + References + ---------- + .. [1] 'A Survey of Methods of Computing Minimax and Near-Minimax + Polynomial Approximations for Functions of a Single Independent + Variable' by W. Fraser, `Journal of the ACM` (JACM), + vol. 12(3), pp. 295-314, (1965) + .. [2] MathWorld, 'Chebyshev Approximation Formula', + https://mathworld.wolfram.com/ChebyshevApproximationFormula.html + .. [3] Wikipedia, 'Chebyshev nodes', + http://en.wikipedia.org/wiki/Chebyshev_nodes + """ - deg = np.asarray(deg) + deg = np.atleast_1d(deg) # check arguments. - if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0: + if deg.ndim != 1 or deg.dtype.kind not in 'iu' or deg.size == 0: raise TypeError("deg must be an int") - if deg < 0: + if np.any(deg < 0): raise ValueError("expected deg >= 0") - order = deg + 1 - xcheb = chebpts1(order) - yfunc = func(xcheb, *args) - m = chebvander(xcheb, deg) - c = np.dot(m.T, yfunc) - c[0] /= order - c[1:] /= 0.5*order + ndim = len(order) - return c + if domain is None: + domain = chebdomain.tolist() * ndim + domain = np.atleast_2d(domain).reshape((-1, 2)) + x_i = [pu.mapdomain(chebpts1(o), chebdomain, d) + for o, d in zip(order, domain)] + c = func(*np.meshgrid(*x_i), *args) + + for i in range(ndim): + c = __dct(c[..., ::-1]) + c[..., 0] /= 2. + if i < ndim - 1: + c = np.moveaxis(c, source=-1, destination=0) + return c / np.product(order) def chebgauss(deg): @@ -2060,13 +2222,12 @@ def interpolate(cls, func, deg, domain=None, args=()): Notes ----- - See `numpy.polynomial.chebfromfunction` for more details. + See `numpy.polynomial.chebyshev.chebinterpolate` for more details. """ if domain is None: domain = cls.domain - xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args) - coef = chebinterpolate(xfunc, deg) + coef = chebinterpolate(func, deg, args=args, domain=domain) return cls(coef, domain=domain) # Virtual properties diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py index 2f54bebfdb27..24c4906e6c65 100644 --- a/numpy/polynomial/tests/test_chebyshev.py +++ b/numpy/polynomial/tests/test_chebyshev.py @@ -506,6 +506,20 @@ def powx(x, p): c = cheb.chebinterpolate(powx, deg, (p,)) assert_almost_equal(cheb.chebval(x, c), powx(x, p), decimal=12) + def test_2d_approximation(self): + + def powxy(x, y, p): + return (x+y)**p + + x = np.linspace(-1, 1, 10) + x1, x2 = np.meshgrid(x, x) + for deg in range(0, 10): + for p in range(0, deg + 1): + c = cheb.chebinterpolate(powxy, (deg, deg), (p,)) + assert_almost_equal(cheb.chebval2d(x1, x2, c), + powxy(x1, x2, p), + decimal=12) + class TestCompanion: