Skip to content

Navigation Menu

Sign in
Appearance settings

Search code, repositories, users, issues, pull requests...

Provide feedback

We read every piece of feedback, and take your input very seriously.

Saved searches

Use saved searches to filter your results more quickly
Appearance settings

Commit 2953284

Browse filesBrowse files
pbrodpbrod
committed
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.
1 parent 7ca1d1a commit 2953284
Copy full SHA for 2953284

File tree

2 files changed

+205
-30
lines changed
Filter options

2 files changed

+205
-30
lines changed

‎numpy/polynomial/chebyshev.py

Copy file name to clipboardExpand all lines: numpy/polynomial/chebyshev.py
+191-30Lines changed: 191 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -1777,27 +1777,166 @@ def chebroots(c):
17771777
return r
17781778

17791779

1780-
def chebinterpolate(func, deg, args=()):
1780+
def __dct(x, n=None):
1781+
"""
1782+
Discrete Cosine Transform type II
1783+
N-1
1784+
y[k] = 2 * sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
1785+
n=0
1786+
Examples
1787+
--------
1788+
>>> import numpy as np
1789+
>>> x = np.arange(5)
1790+
>>> x1 = np.arange(6)
1791+
>>> y = x + 1j*x
1792+
>>> y1 = x1 + 1j*x1
1793+
1794+
>>> np.allclose(__dct(x), [20. , -9.95959314, 0. , -0.89805595, 0. ])
1795+
True
1796+
1797+
>>> np.allclose(__dct(x1), [30. , -14.41953704, 0. , -1.41421356, 0. , -0.27740142])
1798+
True
1799+
1800+
>>> np.allclose(__dct(y), [20. +20.j , -9.95959314 -9.95959314j,
1801+
... 0. +0.j , -0.89805595 -0.89805595j,
1802+
... 0. +0.j ])
1803+
True
1804+
1805+
>>> np.allclose(__dct(y1), [ 30. +30.j , -14.41953704-14.41953704j,
1806+
... 0. +0.j , -1.41421356 -1.41421356j,
1807+
... 0. +0.j , -0.27740142 -0.27740142j])
1808+
True
1809+
1810+
>>> np.allclose(x, __idct(__dct(x)))
1811+
True
1812+
1813+
>>> np.allclose(x, __dct(__idct(x)))
1814+
True
1815+
1816+
References
1817+
----------
1818+
.. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
1819+
Makhoul, `IEEE Transactions on acoustics, speech and signal
1820+
processing` vol. 28(1), pp. 27-34,
1821+
:doi:`10.1109/TASSP.1980.1163351` (1980).
1822+
.. [2] Wikipedia, "Discrete cosine transform",
1823+
https://en.wikipedia.org/wiki/Discrete_cosine_transform
1824+
"""
1825+
fft = np.fft.fft
1826+
x = np.atleast_1d(x)
1827+
1828+
if n is None:
1829+
n = x.shape[-1]
1830+
1831+
if x.shape[-1] < n:
1832+
n_shape = x.shape[:-1] + (n - x.shape[-1],)
1833+
xx = np.hstack((x, np.zeros(n_shape)))
1834+
else:
1835+
xx = x[..., :n]
1836+
1837+
real_x = np.all(np.isreal(xx))
1838+
if (real_x and (np.remainder(n, 2) == 0)):
1839+
xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2])))
1840+
else:
1841+
xp = fft(np.hstack((xx, xx[..., ::-1])))[..., :n]
1842+
1843+
w = np.exp(-1j * np.arange(n) * np.pi / (2 * n))
1844+
1845+
y = xp * w
1846+
1847+
if real_x:
1848+
return y.real
1849+
return y
1850+
1851+
1852+
def __idct(x, n=None):
1853+
"""
1854+
Inverse Discrete Cosine Transform type II
1855+
N-1
1856+
y[k] = 1/N sum w[n]*x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
1857+
n=0
1858+
w[0] = 1/2
1859+
w[n] = 1 for n>0
1860+
1861+
Examples
1862+
--------
1863+
>>> import numpy as np
1864+
>>> x = np.arange(5)
1865+
>>> np.allclose(x, __idct(__dct(x)))
1866+
True
1867+
1868+
>>> np.allclose(x, __dct(__idct(x)))
1869+
True
1870+
1871+
References
1872+
----------
1873+
.. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
1874+
Makhoul, `IEEE Transactions on acoustics, speech and signal
1875+
processing` vol. 28(1), pp. 27-34,
1876+
:doi:`10.1109/TASSP.1980.1163351` (1980).
1877+
.. [2] Wikipedia, "Discrete cosine transform",
1878+
https://en.wikipedia.org/wiki/Discrete_cosine_transform
1879+
"""
1880+
1881+
ifft = np.fft.ifft
1882+
x = np.atleast_1d(x)
1883+
1884+
if n is None:
1885+
n = x.shape[-1]
1886+
1887+
w = np.exp(1j * np.arange(n) * np.pi / (2 * n))
1888+
1889+
if x.shape[-1] < n:
1890+
n_shape = x.shape[:-1] + (n - x.shape[-1],)
1891+
xx = np.hstack((x, np.zeros(n_shape))) * w
1892+
else:
1893+
xx = x[..., :n] * w
1894+
1895+
real_x = np.all(np.isreal(x))
1896+
if (real_x and (np.remainder(n, 2) == 0)):
1897+
xx[..., 0] = xx[..., 0] * 0.5
1898+
yp = ifft(xx)
1899+
y = np.zeros(xx.shape, dtype=complex)
1900+
y[..., ::2] = yp[..., :n // 2]
1901+
y[..., ::-2] = yp[..., n // 2::]
1902+
else:
1903+
yp = ifft(np.hstack((xx,
1904+
np.zeros_like(xx[..., 0]),
1905+
np.conj(xx[..., :0:-1]))))
1906+
y = yp[..., :n]
1907+
1908+
if real_x:
1909+
return y.real
1910+
return y
1911+
1912+
1913+
def chebinterpolate(func, deg, args=(), domain=None):
17811914
"""Interpolate a function at the Chebyshev points of the first kind.
17821915
1783-
Returns the Chebyshev series that interpolates `func` at the Chebyshev
1784-
points of the first kind in the interval [-1, 1]. The interpolating
1785-
series tends to a minmax approximation to `func` with increasing `deg`
1786-
if the function is continuous in the interval.
1916+
Returns the Chebyshev series that interpolates the multi-dimensional
1917+
function `func` at the Chebyshev points of the first kind scaled and
1918+
shifted to the `domain`. The interpolating series tends to a minmax
1919+
approximation of `func` with increasing `deg` if the function is
1920+
continuous in the input domain.
17871921
17881922
.. versionadded:: 1.14.0
17891923
17901924
Parameters
17911925
----------
17921926
func : function
1793-
The function to be approximated. It must be a function of a single
1794-
variable of the form ``f(x, a, b, c...)``, where ``a, b, c...`` are
1795-
extra arguments passed in the `args` parameter.
1796-
deg : int
1797-
Degree of the interpolating polynomial
1927+
The function to be approximated. It must be a function of the form
1928+
``f(x1, x2,..., xn, a, b, c...)``, where ``x1, x2,..., xn`` are the
1929+
n-dimensional coordinates and ``a, b, c...`` are extra arguments
1930+
passed in the `args` parameter.
1931+
deg : list of int, size ndim
1932+
Degree of the interpolating polynomial used for each dimension.
17981933
args : tuple, optional
17991934
Extra arguments to be used in the function call. Default is no extra
18001935
arguments.
1936+
domain : array-like, shape ndim x 2, optional
1937+
defines the domain ``[a1, b1] x [a2, b2] x ... x [an, bn]`` over
1938+
which `func` is interpolated. The default is None, in which case
1939+
the domain is [-1,1] * len(deg))
18011940
18021941
Returns
18031942
-------
@@ -1807,11 +1946,19 @@ def chebinterpolate(func, deg, args=()):
18071946
18081947
Examples
18091948
--------
1949+
>>> import numpy as np
18101950
>>> import numpy.polynomial.chebyshev as C
1811-
>>> C.chebfromfunction(lambda x: np.tanh(x) + 0.5, 8)
1812-
array([ 5.00000000e-01, 8.11675684e-01, -9.86864911e-17,
1813-
-5.42457905e-02, -2.71387850e-16, 4.51658839e-03,
1814-
2.46716228e-17, -3.79694221e-04, -3.26899002e-16])
1951+
>>> c1 = C.chebinterpolate(lambda x: np.tanh(x) + 0.5, 8)
1952+
>>> x = C.chebpts1(9)
1953+
>>> np.all(np.abs(C.chebval(x, c1) - np.tanh(x)-0.5) < 1e-15)
1954+
True
1955+
1956+
Two-dimensional interpolation to machine accuracy
1957+
>>> x1,x2 = np.meshgrid(x,x)
1958+
>>> c2 = C.chebinterpolate(lambda x,y: np.tanh(x+y) + 0.5, deg=(8, 8))
1959+
>>> np.all(np.abs(C.chebval2d(x1, x2, c2) - np.tanh(x1+x2)-0.5) < 1e-15)
1960+
True
1961+
18151962
18161963
Notes
18171964
-----
@@ -1824,24 +1971,41 @@ def chebinterpolate(func, deg, args=()):
18241971
zero. For instance, if the function is even then the coefficients of the
18251972
terms of odd degree in the result can be set to zero.
18261973
1974+
References
1975+
----------
1976+
.. [1] 'A Survey of Methods of Computing Minimax and Near-Minimax Polynomial
1977+
Approximations for Functions of a Single Independent Variable' by
1978+
W. Fraser, `Journal of the ACM` (JACM),
1979+
vol. 12(3), pp. 295-314, (1965)
1980+
.. [2] MathWorld, a Wolfram Web Resource, 'Chebyshev Approximation Formula',
1981+
https://mathworld.wolfram.com/ChebyshevApproximationFormula.html
1982+
.. [3] Wikipedia, 'Chebyshev nodes',
1983+
http://en.wikipedia.org/wiki/Chebyshev_nodes
1984+
18271985
"""
1828-
deg = np.asarray(deg)
1986+
deg = np.atleast_1d(deg)
18291987

18301988
# check arguments.
1831-
if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
1989+
if deg.ndim != 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
18321990
raise TypeError("deg must be an int")
1833-
if deg < 0:
1991+
if np.any(deg < 0):
18341992
raise ValueError("expected deg >= 0")
1835-
18361993
order = deg + 1
1837-
xcheb = chebpts1(order)
1838-
yfunc = func(xcheb, *args)
1839-
m = chebvander(xcheb, deg)
1840-
c = np.dot(m.T, yfunc)
1841-
c[0] /= order
1842-
c[1:] /= 0.5*order
1994+
ndim = len(order)
18431995

1844-
return c
1996+
if domain is None:
1997+
domain = chebdomain.tolist() * ndim
1998+
domain = np.atleast_2d(domain).reshape((-1, 2))
1999+
x_i = [pu.mapdomain(chebpts1(o), chebdomain, d)
2000+
for o, d in zip(order, domain)]
2001+
c = func(*np.meshgrid(*x_i), *args)
2002+
2003+
for i in range(ndim):
2004+
c = __dct(c[..., ::-1])
2005+
c[..., 0] /= 2.
2006+
if i < ndim - 1:
2007+
c = np.moveaxis(c, source=-1, destination=0)
2008+
return c / np.product(order)
18452009

18462010

18472011
def chebgauss(deg):
@@ -2060,13 +2224,10 @@ def interpolate(cls, func, deg, domain=None, args=()):
20602224
20612225
Notes
20622226
-----
2063-
See `numpy.polynomial.chebfromfunction` for more details.
2227+
See `numpy.polynomial.chebyshev.chebinterpolate` for more details.
20642228
20652229
"""
2066-
if domain is None:
2067-
domain = cls.domain
2068-
xfunc = lambda x: func(pu.mapdomain(x, cls.window, domain), *args)
2069-
coef = chebinterpolate(xfunc, deg)
2230+
coef = chebinterpolate(func, deg, args=args, domain=domain)
20702231
return cls(coef, domain=domain)
20712232

20722233
# Virtual properties

‎numpy/polynomial/tests/test_chebyshev.py

Copy file name to clipboardExpand all lines: numpy/polynomial/tests/test_chebyshev.py
+14Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -506,6 +506,20 @@ def powx(x, p):
506506
c = cheb.chebinterpolate(powx, deg, (p,))
507507
assert_almost_equal(cheb.chebval(x, c), powx(x, p), decimal=12)
508508

509+
def test_2d_approximation(self):
510+
511+
def powxy(x, y, p):
512+
return (x+y)**p
513+
514+
x = np.linspace(-1, 1, 10)
515+
x1, x2 = np.meshgrid(x, x)
516+
for deg in range(0, 10):
517+
for p in range(0, deg + 1):
518+
c = cheb.chebinterpolate(powxy, (deg, deg), (p,))
519+
assert_almost_equal(cheb.chebval2d(x1, x2, c),
520+
powxy(x1, x2, p),
521+
decimal=12)
522+
509523

510524
class TestCompanion:
511525

0 commit comments

Comments
0 (0)
Morty Proxy This is a proxified and sanitized view of the page, visit original site.