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 cb5261f

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 numpy#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 numpy#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 cb5261f
Copy full SHA for cb5261f

File tree

2 files changed

+203
-28
lines changed
Filter options

2 files changed

+203
-28
lines changed

‎numpy/polynomial/chebyshev.py

Copy file name to clipboardExpand all lines: numpy/polynomial/chebyshev.py
+189-28Lines changed: 189 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1777,27 +1777,165 @@ 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+
>>> x = np.arange(5)
1789+
>>> x1 = np.arange(6)
1790+
>>> y = x + 1j*x
1791+
>>> y1 = x1 + 1j*x1
1792+
1793+
>>> np.allclose(__dct(x), [20. , -9.95959314, 0. , -0.89805595, 0. ])
1794+
True
1795+
1796+
>>> np.allclose(__dct(x1),
1797+
... [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+
>>> x = np.arange(5)
1864+
>>> np.allclose(x, __idct(__dct(x)))
1865+
True
1866+
1867+
>>> np.allclose(x, __dct(__idct(x)))
1868+
True
1869+
1870+
References
1871+
----------
1872+
.. [1] 'A Fast Cosine Transform in One and Two Dimensions', by J.
1873+
Makhoul, `IEEE Transactions on acoustics, speech and signal
1874+
processing` vol. 28(1), pp. 27-34,
1875+
:doi:`10.1109/TASSP.1980.1163351` (1980).
1876+
.. [2] Wikipedia, "Discrete cosine transform",
1877+
https://en.wikipedia.org/wiki/Discrete_cosine_transform
1878+
"""
1879+
1880+
ifft = np.fft.ifft
1881+
x = np.atleast_1d(x)
1882+
1883+
if n is None:
1884+
n = x.shape[-1]
1885+
1886+
w = np.exp(1j * np.arange(n) * np.pi / (2 * n))
1887+
1888+
if x.shape[-1] < n:
1889+
n_shape = x.shape[:-1] + (n - x.shape[-1],)
1890+
xx = np.hstack((x, np.zeros(n_shape))) * w
1891+
else:
1892+
xx = x[..., :n] * w
1893+
1894+
real_x = np.all(np.isreal(x))
1895+
if (real_x and (np.remainder(n, 2) == 0)):
1896+
xx[..., 0] = xx[..., 0] * 0.5
1897+
yp = ifft(xx)
1898+
y = np.zeros(xx.shape, dtype=complex)
1899+
y[..., ::2] = yp[..., :n // 2]
1900+
y[..., ::-2] = yp[..., n // 2::]
1901+
else:
1902+
yp = ifft(np.hstack((xx,
1903+
np.zeros_like(xx[..., 0]),
1904+
np.conj(xx[..., :0:-1]))))
1905+
y = yp[..., :n]
1906+
1907+
if real_x:
1908+
return y.real
1909+
return y
1910+
1911+
1912+
def chebinterpolate(func, deg, args=(), domain=None):
17811913
"""Interpolate a function at the Chebyshev points of the first kind.
17821914
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.
1915+
Returns the Chebyshev series that interpolates the multi-dimensional
1916+
function `func` at the Chebyshev points of the first kind scaled and
1917+
shifted to the `domain`. The interpolating series tends to a minmax
1918+
approximation of `func` with increasing `deg` if the function is
1919+
continuous in the input domain.
17871920
17881921
.. versionadded:: 1.14.0
17891922
17901923
Parameters
17911924
----------
17921925
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
1926+
The function to be approximated. It must be a function of the form
1927+
``f(x1, x2,..., xn, a, b, c...)``, where ``x1, x2,..., xn`` are the
1928+
n-dimensional coordinates and ``a, b, c...`` are extra arguments
1929+
passed in the `args` parameter.
1930+
deg : list of int, size ndim
1931+
Degree of the interpolating polynomial used for each dimension.
17981932
args : tuple, optional
17991933
Extra arguments to be used in the function call. Default is no extra
18001934
arguments.
1935+
domain : array-like, shape ndim x 2, optional
1936+
defines the domain ``[a1, b1] x [a2, b2] x ... x [an, bn]`` over
1937+
which `func` is interpolated. The default is None, in which case
1938+
the domain is [-1,1] * len(deg))
18011939
18021940
Returns
18031941
-------
@@ -1808,10 +1946,17 @@ def chebinterpolate(func, deg, args=()):
18081946
Examples
18091947
--------
18101948
>>> 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])
1949+
>>> c1 = C.chebinterpolate(lambda x: np.tanh(x) + 0.5, 8)
1950+
>>> x = C.chebpts1(9)
1951+
>>> np.all(np.abs(C.chebval(x, c1) - np.tanh(x)-0.5) < 1e-15)
1952+
True
1953+
1954+
Two-dimensional interpolation to machine accuracy
1955+
>>> x1,x2 = np.meshgrid(x,x)
1956+
>>> c2 = C.chebinterpolate(lambda x,y: np.tanh(x+y) + 0.5, deg=(8, 8))
1957+
>>> np.all(np.abs(C.chebval2d(x1, x2, c2) - np.tanh(x1+x2)-0.5) < 1e-15)
1958+
True
1959+
18151960
18161961
Notes
18171962
-----
@@ -1824,24 +1969,41 @@ def chebinterpolate(func, deg, args=()):
18241969
zero. For instance, if the function is even then the coefficients of the
18251970
terms of odd degree in the result can be set to zero.
18261971
1972+
References
1973+
----------
1974+
.. [1] 'A Survey of Methods of Computing Minimax and Near-Minimax
1975+
Polynomial Approximations for Functions of a Single Independent
1976+
Variable' by W. Fraser, `Journal of the ACM` (JACM),
1977+
vol. 12(3), pp. 295-314, (1965)
1978+
.. [2] MathWorld, 'Chebyshev Approximation Formula',
1979+
https://mathworld.wolfram.com/ChebyshevApproximationFormula.html
1980+
.. [3] Wikipedia, 'Chebyshev nodes',
1981+
http://en.wikipedia.org/wiki/Chebyshev_nodes
1982+
18271983
"""
1828-
deg = np.asarray(deg)
1984+
deg = np.atleast_1d(deg)
18291985

18301986
# check arguments.
1831-
if deg.ndim > 0 or deg.dtype.kind not in 'iu' or deg.size == 0:
1987+
if deg.ndim != 1 or deg.dtype.kind not in 'iu' or deg.size == 0:
18321988
raise TypeError("deg must be an int")
1833-
if deg < 0:
1989+
if np.any(deg < 0):
18341990
raise ValueError("expected deg >= 0")
1835-
18361991
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
1992+
ndim = len(order)
18431993

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

18462008

18472009
def chebgauss(deg):
@@ -2060,13 +2222,12 @@ def interpolate(cls, func, deg, domain=None, args=()):
20602222
20612223
Notes
20622224
-----
2063-
See `numpy.polynomial.chebfromfunction` for more details.
2225+
See `numpy.polynomial.chebyshev.chebinterpolate` for more details.
20642226
20652227
"""
20662228
if domain is None:
20672229
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.