danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added chebfit_dct

Browse Source 
chebvandernd
chebfitnd
chebvalnd
chebgridnd
master
Per.Andreas.Brodtkorb 10 years ago
parent 34330797b3
commit b992e8c979
1 changed files with 438 additions and 74 deletions
Show Stats Download Patch File Download Diff File

512
pywafo/src/wafo/polynomial.py
Unescape Escape View File

@ -17,14 +17,16 @@
# Licence: LGPL # Licence: LGPL
# ------------------------------------------------------------------------- # -------------------------------------------------------------------------
# !/usr/bin/env python # !/usr/bin/env python
import warnings
from numpy.polynomial import polyutils as pu
from plotbackend import plotbackend as plt from plotbackend import plotbackend as plt
import numpy as np import numpy as np
from numpy.fft import fft, ifft
from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange, from numpy import (zeros, ones, zeros_like, array, asarray, newaxis, arange,
logical_or, any, pi, cos, round, diff, all, exp, atleast_1d, logical_or, any, pi, cos, round, diff, all, exp, atleast_1d,
where, extract, linalg, sign, concatenate, floor, isreal, where, extract, linalg, sign, concatenate, floor, isreal,
conj, remainder, linspace, sum) conj, remainder, linspace, sum, meshgrid, hstack)
from scipy.fftpack import dct, idct as _idct
from numpy.lib.polynomial import * # @UnusedWildImport from numpy.lib.polynomial import * # @UnusedWildImport
from scipy.misc.common import pade # @UnresolvedImport from scipy.misc.common import pade # @UnresolvedImport
__all__ = np.lib.polynomial.__all__ __all__ = np.lib.polynomial.__all__
@ -1171,64 +1173,101 @@ def chebfit(fun, n=10, a=-1, b=1, trace=False):
f = fun f = fun
n = len(f) n = len(f)
# N-1 # N-1
# c(k) = (2/N) sum w(n) f(n)*cos(pi*k*(2n+1)/(2N)), 0 <= k < N. # c[k] = (2/N) sum w[n] f[n]*cos(pi*k*(2n+1)/(2N)), 0 <= k < N.
# n=0 # n=0
# #
# w(0) = 0.5, w(n)=1 for n>0 # w[0] = 0.5, w[n]=1 for n>0
ck = dct(f[::-1]) / n ck = dct(f[::-1]) / n
ck[0] = ck[0] / 2. ck[0] = ck[0] / 2.
return ck[::-1] return ck[::-1]
def dct(x, n=None): def chebfit_dct(f, n=(10, ), domain=None):
""" """
Discrete Cosine Transform Fit Chebyshev series to N-dimensional function
so that f(x1, x2,..., xn) can be approximated by:
N-1 .. math:: f(x_1, x_2,...,x_n) = \\sum_{i,j,...k} c_i T_i(x_1)*...*c_k T_k(x_n) ,
y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
n=0
Examples where Tk is the k'th Chebyshev polynomial of the first kind.
--------
>>> import numpy as np
>>> x = np.arange(5)
>>> np.abs(x-idct(dct(x)))<1e-14
array([ True, True, True, True, True], dtype=bool)
>>> np.abs(x-dct(idct(x)))<1e-14
array([ True, True, True, True, True], dtype=bool)
Reference Parameters
--------- ----------
http://en.wikipedia.org/wiki/Discrete_cosine_transform f : callable
http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ function to approximate
""" n : list of integers, optional
number of base points (abscissas) used for each dimension.
Default n=10 (maximum 50)
domain : vector [a1,b1,a2,b2 ,..., an, bn], optional
defining the rectangle [a1, b1] x [a2, b2] x ...x [an, bn].
(default domain = (-1,1) * len(n))
Returns
-------
ck : ndarray
polynomial coefficients in Chebychev form.
x = atleast_1d(x) Examples
--------
Fit exponential function
if n is None: >>> import matplotlib.pyplot as plt
n = x.shape[-1] >>> domain = (0, 2)
>>> ck = chebfit_dct(np.exp, 7, domain)
>>> np.allclose(ck, [3.44152387e+00, 3.07252345e+00, 7.38000848e-01,
... 1.20520053e-01, 1.48805268e-02, 1.47579673e-03,
... 1.21719524e-04])
True
>>> x1 = map_to_interval(chebroot(9), *domain)
>>> ck1 = chebfit(np.exp(x1))
>>> np.allclose(ck1, [5.40019009e-07, 8.69418381e-06, 1.22261037e-04,
... 1.47582673e-03, 1.48805283e-02, 1.20520053e-01,
... 7.38000848e-01, 3.07252345e+00, 3.44152387e+00])
True
if x.shape[-1] < n: x = np.linspace(0,4)
n_shape = x.shape[:-1] + (n - x.shape[-1],) h = plt.plot(x, np.exp(x), 'r', x, chebvalnd(ck, x,ck,a,b), 'g.')
xx = hstack((x, zeros(n_shape))) h = plt.plot(x, np.exp(x), 'r', x, chebvalnd(ck1, x,ck1,a,b),'b.')
else: plt.close()
xx = x[..., :n]
real_x = all(isreal(xx)) See also
if (real_x and (remainder(n, 2) == 0)): --------
xp = 2 * fft(hstack((xx[..., ::2], xx[..., ::-2]))) chebval, chebvalnd
else:
xp = fft(hstack((xx, xx[..., ::-1])))
xp = xp[..., :n]
w = exp(-1j * arange(n) * pi / (2 * n)) Reference
---------
http://en.wikipedia.org/wiki/Chebyshev_nodes
http://mathworld.wolfram.com/ChebyshevApproximationFormula.html
y = xp * w W. Fraser (1965)
"A Survey of Methods of Computing Minimax and Near-Minimax Polynomial
Approximations for Functions of a Single Independent Variable"
Journal of the ACM (JACM), Vol. 12 , Issue 3, pp 295 - 314
"""
n = np.atleast_1d(n)
if np.any(n > 50):
warnings.warn('CHEBFIT should only be used for n<50')
if real_x: if hasattr(f, '__call__'):
return y.real if domain is None:
domain = (-1,1) * len(n)
domain = np.atleast_2d(domain).reshape((-1, 2))
xi = [map_to_interval(chebroot(ni), d[0], d[1])
for ni, d in zip(n, domain)]
Xi = np.meshgrid(*xi)
ck = f(*Xi)
else: else:
return y ck = f
n = ck.shape
ndim = len(n)
for i in range(ndim):
ck = dct(ck[...,::-1])
ck[..., 0] = ck[..., 0] / 2.
if i < ndim-1:
ck = np.rollaxis(ck, axis=-1)
return ck / np.product(n)
def idct(x, n=None): def idct(x, n=None):
@ -1256,35 +1295,7 @@ def idct(x, n=None):
http://en.wikipedia.org/wiki/Discrete_cosine_transform http://en.wikipedia.org/wiki/Discrete_cosine_transform
http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
""" """
return _idct(x, n=n, norm=None)*0.5/len(x)
x = atleast_1d(x)
if n is None:
n = x.shape[-1]
w = exp(1j * arange(n) * pi / (2 * n))
if x.shape[-1] < n:
n_shape = x.shape[:-1] + (n - x.shape[-1],)
xx = hstack((x, zeros(n_shape))) * w
else:
xx = x[..., :n] * w
real_x = all(isreal(x))
if (real_x and (remainder(n, 2) == 0)):
xx[..., 0] = xx[..., 0] * 0.5
yp = ifft(xx)
y = zeros(xx.shape, dtype=complex)
y[..., ::2] = yp[..., :n / 2]
y[..., ::-2] = yp[..., n / 2::]
else:
yp = ifft(hstack((xx, zeros_like(xx[..., 0]), conj(xx[..., :0:-1]))))
y = yp[..., :n]
if real_x:
return y.real
else:
return y
def _chebval(x, ck, kind=1): def _chebval(x, ck, kind=1):
@ -1792,7 +1803,7 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
smallest_devmax = BIG smallest_devmax = BIG
ncof = m + k + 1 ncof = m + k + 1
# % Number of points where function is evaluated, i.e. fineness of mesh # Number of points where function is evaluated, i.e. fineness of mesh
npt = NFAC * ncof npt = NFAC * ncof
if x is None: if x is None:
@ -1848,7 +1859,7 @@ def padefitlsq(fun, m, k, a=-1, b=1, trace=False, x=None, end_points=True):
wt = np.abs(ee) wt = np.abs(ee)
devmax = max(wt) devmax = max(wt)
mad = wt.mean() # % mean absolute deviation mad = wt.mean() # mean absolute deviation
# Save only the best coefficients found # Save only the best coefficients found
if (devmax <= smallest_devmax): if (devmax <= smallest_devmax):
@ -1940,9 +1951,362 @@ def test_docstrings():
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE) doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
def chebvandernd(deg, *xi):
"""Pseudo-Vandermonde matrix of given degrees.
Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
points `(x1, x2, ..., xn)`. If `l, m, n` are the given degrees in
`x1, x2, x3`, then The pseudo-Vandermonde matrix is defined by
.. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = T_i(x1)*T_j(x2)*T_k(x3),
where `0 <= i <= l`, `0 <= j <= m`, and `0 <= k <= n`. The leading
indices of `V` index the points `(x, y, z)` and the last index encodes
the degrees of the Chebyshev polynomials.
If ``V = chebvandernd([xdeg, ydeg, zdeg], x, y, z)``, then the columns
of `V` correspond to the elements of a 3-D coefficient array `c` of
shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
.. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
and ``np.dot(V, c.flat)`` and ``chebvalnd(c, x, y, z)`` will be the
same up to roundoff. This equivalence is useful both for least squares
fitting and for the evaluation of a large number of N-D Chebyshev
series of the same degrees and sample points.
Parameters
----------
deg : list of ints
List of maximum degrees of the form [x1_deg, x2_deg, ...,xn_deg].
x1, x2, ..., xn : array_like
Arrays of point coordinates, all of the same shape. The dtypes will
be converted to either float64 or complex128 depending on whether
any of the elements are complex. Scalars are converted to 1-D
arrays.
Returns
-------
vander : ndarray
The shape of the returned matrix is ``x1.shape + (order,)``, where
:math:`order = (deg[0]+1)*(deg([1]+1)*...*(deg[n-1]+1)`. The dtype will
be the same as the converted `x1`, `x2`, ... `xn`.
See Also
--------
chebvander, chebvalnd, chebfitnd
"""
ideg = [int(d) for d in deg]
is_valid = np.array([id == d and id >= 0 for id, d in zip(ideg, deg)])
if np.any(is_valid != 1):
raise ValueError("degrees must be non-negative integers")
ndim = len(xi)
if len(ideg)!=ndim:
raise ValueError('length of deg must be the same as number of dimensions')
xi = np.array(xi, copy=0) + 0.0
chebvander = np.polynomial.chebyshev.chebvander
shape0 = xi[0].shape
s0 = (1,) * ndim
vxi = [chebvander(x, d).reshape(shape0 + s0[:i] + (-1,) + s0[i + 1::])
for i, (d, x) in enumerate(zip(ideg, xi))]
v = reduce(np.multiply, vxi)
return v.reshape(v.shape[:-ndim] + (-1,))
def chebfitnd(xi, f, deg, rcond=None, full=False, w=None):
"""
Least squares fit of Chebyshev series to N-dimensional data.
Return the coefficients of a Chebyshev series of degree `deg` that is the
least squares fit to the data values `f` given at points `x1`, `x2`,..., `xn`
The fitted polynomial(s) are in the form
.. math:: p(x,y) = c_00 + c_11 * T_1(x)*T_1(y) + ..c_ij * T_i(x)*T_j(y).
+ c_nm * T_n(x)*T_m(y),
where `n`, `m` is `deg`.
Parameters
----------
xi: tuple
x1-, x2-,....xn-coordinates of the sample points.
f : array_like
function values at the sample points ``(x1[i], x2[i], ..., xn[i])``.
deg : list
Degrees of the fitting series in the x1, x2, ..., xn directions, respectively.
rcond : float, optional
Relative condition number of the fit. Singular values smaller than
this relative to the largest singular value will be ignored. The
default value is size(x1)*eps, where eps is the relative precision of
the float type, about 2e-16 in most cases.
full : bool, optional
Switch determining nature of return value. When it is False (the
default) just the coefficients are returned, when True diagnostic
information from the singular value decomposition is also returned.
w : array_like, optional
Weights. If not None, the contribution of each point
``(x1[i], x2[i], ..., xn[i])`` to the fit is weighted by `w[i]`.
Ideally the weights are chosen so that the errors of the products
``w[i]*f[i]`` all have the same variance. The default value is None.
Returns
-------
coef : ndarray, shape (M1, M2,..., Mn)
Chebyshev coefficients ordered from low to high.
[residuals, rank, singular_values, rcond] : list
These values are only returned if `full` = True
resid -- sum of squared residuals of the least squares fit
rank -- the numerical rank of the scaled Vandermonde matrix
sv -- singular values of the scaled Vandermonde matrix
rcond -- value of `rcond`.
For more details, see `linalg.lstsq`.
Warns
-----
RankWarning
The rank of the coefficient matrix in the least-squares fit is
deficient. The warning is only raised if `full` = False. The
warnings can be turned off by
>>> import warnings
>>> warnings.simplefilter('ignore', RankWarning)
See Also
--------
chebvalnd, chebgridnd
Notes
-----
The solution is the coefficients of the Chebyshev series `p` that
minimizes the sum of the weighted squared errors
.. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
where :math:`w_j` are the weights. This problem is solved by setting up
as the (typically) overdetermined matrix equation
.. math:: V(x, y) * c = w * y,
where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
coefficients to be solved for, `w` are the weights, and `y` are the
observed values. This equation is then solved using the singular value
decomposition of `V`.
If some of the singular values of `V` are so small that they are
neglected, then a `RankWarning` will be issued. This means that the
coefficient values may be poorly determined. Using a lower order fit
will usually get rid of the warning. The `rcond` parameter can also be
set to a value smaller than its default, but the resulting fit may be
spurious and have large contributions from roundoff error.
Fits using Chebyshev series are usually better conditioned than fits
using power series, but much can depend on the distribution of the
sample points and the smoothness of the data. If the quality of the fit
is inadequate splines may be a good alternative.
References
----------
.. [1] Wikipedia, "Curve fitting",
http://en.wikipedia.org/wiki/Curve_fitting
Examples
--------
"""
xi_ = np.array(xi, copy=0) + 0.0
z = np.array(f)
degrees = np.asarray(deg, dtype=int)
orders = degrees + 1
order = np.product(orders)
ndims = np.array([x.ndim for x in xi])
ndim = len(ndims)
sizes = np.array([x.size for x in xi])
if np.any(ndims!=ndim) or z.ndim!=ndim:
raise TypeError("expected %dD array for x1, x2,...,xn and f" % ndim)
if np.any(sizes == 0):
raise TypeError("expected non-empty vector for xi")
lhs = chebvandernd(degrees, *xi).reshape((-1, order))
rhs = z.ravel()
if w is not None:
w = np.asarray(w).ravel() + 0.0
if len(lhs) != len(w):
raise TypeError("expected x and w to have same length")
lhs = lhs * w
rhs = rhs * w
if rcond is None:
rcond = xi[0].size * np.finfo(x.dtype).eps
if issubclass(lhs.dtype.type, np.complexfloating):
scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(axis=0))
else:
scl = np.sqrt(np.square(lhs).sum(axis=0))
scl[scl == 0] = 1
# Solve the least squares problem.
c, resids, rank, s = np.linalg.lstsq(lhs/scl, rhs, rcond)
c = (c/scl).reshape(orders)
if rank != order and not full:
msg = "The fit may be poorly conditioned"
warnings.warn(msg, pu.RankWarning)
if full:
return c, [resids, rank, s, rcond]
else:
return c
def chebvalnd(c, *xi):
"""
Evaluate a N-D Chebyshev series at points (x1, x2, ..., xn).
This function returns the values:
.. math:: p(x1,x2,...,xn) = \\sum_{i,j,...,k} c_{i,j,...,k} * T_i(x1) * T_j(x2)*...* T_k(xn)
The parameters `x1`, `x2`, ...., `xn` are converted to arrays only if
they are tuples or a lists, otherwise they are treated as a scalars and
they must have the same shape after conversion. In either case, either
`x1`, `x2`, ..., `xn` or their elements must support multiplication and
addition both with themselves and with the elements of `c`.
If `c` has fewer than N dimensions, ones are implicitly appended to its
shape to make it N-D. The shape of the result will be c.shape[3:] +
x1.shape.
Parameters
----------
c : array_like
Array of coefficients ordered so that the coefficient of the term of
multi-degree i,j,...,k is contained in ``c[i,j,...,k]``. If `c` has
dimension greater than N the remaining indices enumerate multiple sets of
coefficients.
x1, x2,..., xn : array_like, compatible object
The N dimensional series is evaluated at the points
`(x1, x2,...,xn)`, where `x1`, `x2`,..., `xn` must have the same shape.
If any of `x1`, `x2`, ..., `xn` is a list or tuple, it is first converted
to an ndarray, otherwise it is left unchanged and if it isn't an
ndarray it is treated as a scalar.
Returns
-------
values : ndarray, compatible object
The values of the multidimensional polynomial on points formed with
triples of corresponding values from `x`, `y`, and `z`.
See Also
--------
chebval, chebgridnd, chebfitnd
"""
try:
xi = np.array(xi, copy=0)
except:
raise ValueError('x, y, z are incompatible')
chebval = np.polynomial.chebyshev.chebval
c = chebval(xi[0], c)
for x in xi[1:]:
c = chebval(x, c, tensor=False)
return c
def chebgridnd(c, *xi):
"""
Evaluate a N-D Chebyshev series on the Cartesian product of x1, x2,..., xn.
This function returns the values:
.. math:: p(a,b,...) = \\sum_{i,j,...} c_{i,j,...} * T_i(a) * T_j(b) * ...
where the points `(a, b, ...)` consist of all points formed by taking
`a` from `x1`, `b` from `x2`, and so on. The resulting points form
a grid with `x1` in the first dimension, `x2` in the second, and so on.
The parameters `x1`, `x2`, ... and `xn` are converted to arrays only if they
are tuples or a lists, otherwise they are treated as a scalars. In
either case, either `x1`, `x2`,... and `xn` or their elements must support
multiplication and addition both with themselves and with the elements
of `c`.
If `c` has fewer than N dimensions, ones are implicitly appended to
its shape to make it N-D. The shape of the result will be c.shape[3:] +
x1.shape + x2.shape + ... + xn.shape
Parameters
----------
c : array_like
Array of coefficients ordered so that the coefficients for terms of
degree i,j are contained in ``c[i,j]``. If `c` has dimension
greater than two the remaining indices enumerate multiple sets of
coefficients.
x1, x2,..., xn : ndarray, compatible object
1-D arrays representing the coordinates of a grid.
The N dimensional series is evaluated at the points in the
Cartesian product of `x1`, `x2`, ... and `xn`. If `xi`, is a
list or tuple, it is first converted to an ndarray, otherwise it is
left unchanged and, if it isn't an ndarray, it is treated as a
scalar.
Returns
-------
values : ndarray, compatible object
The values of the N dimensional polynomial at points in the Cartesian
product of `x1`, `x2`, ... and `xn`.
See Also
--------
chebval, chebvalnd, chebfitnd
"""
chebval = np.polynomial.chebyshev.chebval
for x in xi:
c = chebval(x, c)
return c
def test_chebfit1d():
n = 63
x = chebroot(n=64, kind=1)
def f(x):
return np.exp(-x**2)
z = f(x)
c = chebfit(f, n=64)[::-1]
xi = np.linspace(-1, 1, 151)
zi = np.polynomial.chebyshev.chebval(xi, c)
#plt.plot(xi, zi,'.', xi, f(xi))
plt.semilogy(xi, np.abs(zi-f(xi)))
plt.show('hold')
def test_chebfit2d():
n = 3
xorder, yorder = n-1, n-1
x = chebroot(n=n, kind=1)
xgrid, ygrid = meshgrid(x, x)
def f(x, y):
return np.exp(-x**2-6*y**2)
zgrid = f(xgrid, ygrid)
#v2d = np.polynomial.chebyshev.chebvander2d(xgrid, ygrid, [xorder,yorder]).reshape((-1, (xorder+1)*(yorder+1)))
#coeff, residuals, rank, s = np.linalg.lstsq(v2d, zgrid.ravel())
#doeff = coeff.reshape(xorder+1,yorder+1)
dcoeff2 = chebfitnd((xgrid, ygrid), zgrid, [xorder,yorder])
dcoeff = chebfit_dct(f, n=(xorder+1,yorder+1))
xi = np.linspace(-1, 1, 151)
Xi,Yi = np.meshgrid(xi, xi)
Zi = f(Xi, Yi)
zzi = chebvalnd(dcoeff, Xi, Yi)
devi = Zi - zzi
# plot residuals
#zz = np.polynomial.chebyshev.chebval2d(xgrid, ygrid, dcoeff)
zz = chebvalnd(dcoeff, xgrid, ygrid)
dev = zgrid - zz
#plt.spy(np.abs(dcoeff)>1e-13)
plt.contourf(xgrid, ygrid, np.abs(dev))
#plt.contourf(Xi, Yi, np.abs(devi))
plt.colorbar()
# plt.semilogy(np.abs(devi.ravel()))
plt.show('hold')
if __name__ == '__main__': if __name__ == '__main__':
if False: # True: # if False: # True: #
main() main()
else: else:
test_docstrings() test_chebfit2d()
# test_docstrings()
# test_polydeg() # test_polydeg()

Write Preview
Loading…
Cancel
Save