From 214d467939b6251ca984c201a8dccc10c0dff582 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Thu, 24 Nov 2011 15:39:41 +0000 Subject: [PATCH] Added improved Sheater-Jones plugin estimate of the smoothing parameter Replaced dct and idct with a call to the ones in scipy.fftpack. Added n-dimensional version dctn and idctn to dctpack.py --- pywafo/src/wafo/dctpack.py | 821 +++++++++++++++++--------- pywafo/src/wafo/kdetools.py | 143 ++++- pywafo/src/wafo/test/test_kdetools.py | 10 + 3 files changed, 700 insertions(+), 274 deletions(-) diff --git a/pywafo/src/wafo/dctpack.py b/pywafo/src/wafo/dctpack.py index de9aa72..79dbd70 100644 --- a/pywafo/src/wafo/dctpack.py +++ b/pywafo/src/wafo/dctpack.py @@ -1,112 +1,150 @@ import numpy as np +from scipy.fftpack import dct as _dct +from scipy.fftpack import idct as _idct __all__ = ['dct', 'idct', 'dctn', 'idctn'] -def dct(x, n=None): - """ - Discrete Cosine Transform - N-1 - y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. - n=0 - - Examples +def dct(x, type=2, n=None, axis=-1, norm='ortho'): + ''' + Return the Discrete Cosine Transform of arbitrary type sequence x. + + Parameters + ---------- + x : array_like + The input array. + type : {1, 2, 3}, optional + Type of the DCT (see Notes). Default type is 2. + n : int, optional + Length of the transform. + axis : int, optional + Axis over which to compute the transform. + norm : {None, 'ortho'}, optional + Normalization mode (see Notes). Default is 'ortho'. + + Returns + ------- + y : ndarray of real + The transformed input array. + + See Also -------- - >>> 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 - --------- - http://en.wikipedia.org/wiki/Discrete_cosine_transform - http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ - """ - 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]))) - xp = xp[..., :n] - - w = np.exp(-1j * np.arange(n) * np.pi / (2 * n)) - - y = xp * w - - if real_x: - return y.real - else: - return y - -def idct(x, n=None): - """ - Inverse Discrete Cosine Transform - + idct + + Notes + ----- + For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to + MATLAB ``dct(x)``. + + There are theoretically 8 types of the DCT, only the first 3 types are + implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the' + Inverse DCT generally refers to DCT type 3. + + type I + ~~~~~~ + There are several definitions of the DCT-I; we use the following + (for ``norm=None``):: + + N-2 + y[k] = x[0] + (-1)**k x[N-1] + 2 * sum x[n]*cos(pi*k*n/(N-1)) + n=1 + + Only None is supported as normalization mode for DCT-I. Note also that the + DCT-I is only supported for input size > 1 + + type II + ~~~~~~~ + There are several definitions of the DCT-II; we use the following + (for ``norm=None``):: + + N-1 - x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. - n=0 - - w(0) = 1/2 - w(n) = 1 for n>0 - - Examples - -------- - >>> 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 - --------- + y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. + n=0 + + If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`:: + + f = sqrt(1/(4*N)) if k = 0, + f = sqrt(1/(2*N)) otherwise. + + Which makes the corresponding matrix of coefficients orthonormal + (``OO' = Id``). + + type III + ~~~~~~~~ + + There are several definitions, we use the following + (for ``norm=None``):: + + N-1 + y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N. + n=1 + + or, for ``norm='ortho'`` and 0 <= k < N:: + + N-1 + y[k] = x[0] / sqrt(N) + sqrt(1/N) * sum x[n]*cos(pi*(k+0.5)*n/N) + n=1 + + The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up + to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of + the orthonormalized DCT-II. + + References + ---------- + http://en.wikipedia.org/wiki/Discrete_cosine_transform - http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ - """ - - 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::] + + '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, http://dx.doi.org/10.1109/TASSP.1980.1163351 (1980). + ''' + farr = np.asfarray + if np.iscomplex(x).any(): + return _dct(farr(x.real), type, n, axis, norm) + 1j*_dct(farr(x.imag), type, n, axis, norm) else: - yp = ifft(np.hstack((xx, np.zeros_like(xx[..., 0]), np.conj(xx[..., :0:-1])))) - y = yp[..., :n] + return _dct(farr(x), type, n, axis, norm) +def idct(x, type=2, n=None, axis=-1, norm='ortho'): + ''' + Return the Inverse Discrete Cosine Transform of an arbitrary type sequence. - if real_x: - return y.real + Parameters + ---------- + x : array_like + The input array. + type : {1, 2, 3}, optional + Type of the DCT (see Notes). Default type is 2. + n : int, optional + Length of the transform. + axis : int, optional + Axis over which to compute the transform. + norm : {None, 'ortho'}, optional + Normalization mode (see Notes). Default is 'ortho'. + + Returns + ------- + y : ndarray of real + The transformed input array. + + See Also + -------- + dct + + Notes + ----- + For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to + matlab ``idct(x)``. + + 'The' IDCT is the IDCT of type 2, which is the same as DCT of type 3. + + IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type 3, + and IDCT of type 3 is the DCT of type 2. For the definition of these types, + see `dct`. + ''' + farr = np.asarray + if np.iscomplex(x).any(): + return _idct(farr(x.real), type, n, axis, norm) + 1j*_idct(farr(x.imag), type, n, axis, norm) else: - return y - -def dctn(y, axis=None, w=None): + return _idct(farr(x), type, n, axis, norm) +def dctn(x, type=2, axis=None, norm='ortho'): ''' DCTN N-D discrete cosine transform. @@ -145,172 +183,381 @@ def dctn(y, axis=None, w=None): -------- idctn, dct, idct ''' - - y = np.atleast_1d(y) - shape0 = y.shape - + y = np.atleast_1d(x) + shape0 = y.shape if axis is None: y = y.squeeze() # Working across singleton dimensions is useless - dimy = y.ndim - if dimy==1: - y = np.atleast_2d(y) - y = y.T - # Some modifications are required if Y is a vector -# if isvector(y): -# if y.shape[0]==1: -# if axis==0: -# return y, None -# elif axis==1: -# axis=0 -# y = y.T -# elif axis==1: -# return y, None - - if w is None: - w = [0,] * dimy - for dim in range(dimy): - if axis is not None and dim!=axis: - continue - n = (dimy==1)*y.size + (dimy>1)*shape0[dim] - #w{dim} = exp(1i*(0:n-1)'*pi/2/n); - w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n)) - - # --- DCT algorithm --- + ndim = y.ndim + isvector = max(shape0)==y.size + if isvector: + if ndim==1: + y = np.atleast_2d(y) + y = y.T + elif y.shape[0]==1: + if axis==0: + return x + elif axis==1: + axis=0 + y = y.T + elif axis==1: + return y + if np.iscomplex(y).any(): - y = np.complex(dctn(np.real(y),axis,w),dctn(np.imag(y),axis,w)) + y = dctn(y.real, type, axis, norm) + 1j*dctn(y.imag, type, axis, norm) else: - for dim in range(dimy): - y = shiftdim(y,1) + y = np.asfarray(y) + for dim in range(ndim): + y = y.transpose(np.roll(range(y.ndim), -1)) + #y = shiftdim(y,1) if axis is not None and dim!=axis: - y = shiftdim(y, 1) continue - siz = y.shape - n = siz[-1] - y = y[...,np.r_[0:n:2, 2*int(n//2)-1:0:-2]] - y = y.reshape((-1,n)) - y = y*np.sqrt(2*n); - y = (np.fft.ifft(y, n=n, axis=1) * w[dim]).real - y[:,0] = y[:,0]/np.sqrt(2) - y = y.reshape(siz) - - #end - #end - - return y.reshape(shape0), w - -def idctn(y, axis=None, w=None): - ''' - IDCTN N-D inverse discrete cosine transform. - X = IDCTN(Y) inverts the N-D DCT transform, returning the original - array if Y was obtained using Y = DCTN(X). - - IDCTN(X,DIM) applies the IDCTN operation across the dimension DIM. - - Class Support - ------------- - Input array can be numeric or logical. The returned array is of class - double. - - Reference - --------- - Narasimha M. et al, On the computation of the discrete cosine - transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936. - - Example - ------- - RGB = imread('autumn.tif'); - I = rgb2gray(RGB); - J = dctn(I); - imshow(log(abs(J)),[]), colormap(jet), colorbar - - The commands below set values less than magnitude 10 in the DCT matrix - to zero, then reconstruct the image using the inverse DCT. - - J(abs(J)<10) = 0; - K = idctn(J); - figure, imshow(I) - figure, imshow(K,[0 255]) - - See also - -------- - dctn, idct, dct - - -- Damien Garcia -- 2009/04, revised 2009/11 - website: www.BiomeCardio.com - - ---------- - [Y,W] = IDCTN(X,DIM,W) uses and returns the weights which are used by - the program. If IDCTN is required for several large arrays of same - size, the weights can be reused to make the algorithm faster. A typical - syntax is the following: - w = []; - for k = 1:10 - [y{k},w] = idctn(x{k},[],w); - end - The weights (w) are calculated during the first call of IDCTN then - reused in the next calls. - ''' - - y = np.atleast_1d(y) - shape0 = y.shape + y = _dct(y, type, norm=norm) + return y.reshape(shape0) +def idctn(x, type=2, axis=None, norm='ortho'): + y = np.atleast_1d(x) + shape0 = y.shape if axis is None: y = y.squeeze() # Working across singleton dimensions is useless - - dimy = y.ndim - if dimy==1: - y = np.atleast_2d(y) - y = y.T - # Some modifications are required if Y is a vector -# if isvector(y): -# if y.shape[0]==1: -# if axis==0: -# return y, None -# elif axis==1: -# axis=0 -# y = y.T -# elif axis==1: -# return y, None -# - - - if w is None: - w = [0,] * dimy - for dim in range(dimy): - if axis is not None and dim!=axis: - continue - n = (dimy==1)*y.size + (dimy>1)*shape0[dim] - #w{dim} = exp(1i*(0:n-1)'*pi/2/n); - w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n)) - # --- IDCT algorithm --- + ndim = y.ndim + isvector = max(shape0)==y.size + if isvector: + if ndim==1: + y = np.atleast_2d(y) + y = y.T + elif y.shape[0]==1: + if axis==0: + return x + elif axis==1: + axis=0 + y = y.T + elif axis==1: + return y + if np.iscomplex(y).any(): - y = np.complex(idctn(np.real(y),axis,w),idctn(np.imag(y),axis,w)) + y = idctn(y.real, type, axis, norm) + 1j*idctn(y.imag, type, axis, norm) else: - for dim in range(dimy): - y = shiftdim(y,1) + y = np.asfarray(y) + for dim in range(ndim): + y = y.transpose(np.roll(range(y.ndim), -1)) + #y = shiftdim(y,1) if axis is not None and dim!=axis: - #y = shiftdim(y, 1) continue - siz = y.shape - n = siz[-1] - - y = y.reshape((-1,n)) * w[dim] - y[:,0] = y[:,0]/np.sqrt(2) - y = (np.fft.ifft(y, n=n, axis=1)).real - y = y * np.sqrt(2*n) - - I = np.empty(n,dtype=int) - I.put(np.r_[0:n:2],np.r_[0:int(n//2)+np.remainder(n,2)]) - I.put(np.r_[1:n:2],np.r_[n-1:int(n//2)-1:-1]) - y = y[:,I] - - y = y.reshape(siz) - - - y = y.reshape(shape0); - return y, w + y = _idct(y, type, norm=norm) + return y.reshape(shape0) + +#def dct(x, n=None): +# """ +# Discrete Cosine Transform +# +# N-1 +# y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. +# n=0 +# +# Examples +# -------- +# >>> 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 +# --------- +# http://en.wikipedia.org/wiki/Discrete_cosine_transform +# http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ +# """ +# 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]))) +# xp = xp[..., :n] +# +# w = np.exp(-1j * np.arange(n) * np.pi / (2 * n)) +# +# y = xp * w +# +# if real_x: +# return y.real +# else: +# return y +# +#def idct(x, n=None): +# """ +# Inverse Discrete Cosine Transform +# +# N-1 +# x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. +# n=0 +# +# w(0) = 1/2 +# w(n) = 1 for n>0 +# +# Examples +# -------- +# >>> 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 +# --------- +# http://en.wikipedia.org/wiki/Discrete_cosine_transform +# http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ +# """ +# +# 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 +# else: +# return y +# +#def dctn(y, axis=None, w=None): +# ''' +# DCTN N-D discrete cosine transform. +# +# Y = DCTN(X) returns the discrete cosine transform of X. The array Y is +# the same size as X and contains the discrete cosine transform +# coefficients. This transform can be inverted using IDCTN. +# +# DCTN(X,axis) applies the DCTN operation across the dimension axis. +# +# Class Support +# ------------- +# Input array can be numeric or logical. The returned array is of class +# double. +# +# Reference +# --------- +# Narasimha M. et al, On the computation of the discrete cosine +# transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936. +# +# Example +# ------- +# RGB = imread('autumn.tif'); +# I = rgb2gray(RGB); +# J = dctn(I); +# imshow(log(abs(J)),[]), colormap(jet), colorbar +# +# The commands below set values less than magnitude 10 in the DCT matrix +# to zero, then reconstruct the image using the inverse DCT. +# +# J(abs(J)<10) = 0; +# K = idctn(J); +# figure, imshow(I) +# figure, imshow(K,[0 255]) +# +# See also +# -------- +# idctn, dct, idct +# ''' +# +# y = np.atleast_1d(y) +# shape0 = y.shape +# +# +# if axis is None: +# y = y.squeeze() # Working across singleton dimensions is useless +# dimy = y.ndim +# if dimy==1: +# y = np.atleast_2d(y) +# y = y.T +# # Some modifications are required if Y is a vector +## if isvector(y): +## if y.shape[0]==1: +## if axis==0: +## return y, None +## elif axis==1: +## axis=0 +## y = y.T +## elif axis==1: +## return y, None +# +# if w is None: +# w = [0,] * dimy +# for dim in range(dimy): +# if axis is not None and dim!=axis: +# continue +# n = (dimy==1)*y.size + (dimy>1)*shape0[dim] +# #w{dim} = exp(1i*(0:n-1)'*pi/2/n); +# w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n)) +# +# # --- DCT algorithm --- +# if np.iscomplex(y).any(): +# y = dctn(np.real(y),axis,w) + 1j*dctn(np.imag(y),axis,w) +# else: +# for dim in range(dimy): +# y = shiftdim(y,1) +# if axis is not None and dim!=axis: +# #y = shiftdim(y, 1) +# continue +# siz = y.shape +# n = siz[-1] +# y = y[...,np.r_[0:n:2, 2*int(n//2)-1:0:-2]] +# y = y.reshape((-1,n)) +# y = y*np.sqrt(2*n); +# y = (np.fft.ifft(y, n=n, axis=1) * w[dim]).real +# y[:,0] = y[:,0]/np.sqrt(2) +# y = y.reshape(siz) +# +# #end +# #end +# +# return y.reshape(shape0), w +# +#def idctn(y, axis=None, w=None): +# ''' +# IDCTN N-D inverse discrete cosine transform. +# X = IDCTN(Y) inverts the N-D DCT transform, returning the original +# array if Y was obtained using Y = DCTN(X). +# +# IDCTN(X,DIM) applies the IDCTN operation across the dimension DIM. +# +# Class Support +# ------------- +# Input array can be numeric or logical. The returned array is of class +# double. +# +# Reference +# --------- +# Narasimha M. et al, On the computation of the discrete cosine +# transform, IEEE Trans Comm, 26, 6, 1978, pp 934-936. +# +# Example +# ------- +# RGB = imread('autumn.tif'); +# I = rgb2gray(RGB); +# J = dctn(I); +# imshow(log(abs(J)),[]), colormap(jet), colorbar +# +# The commands below set values less than magnitude 10 in the DCT matrix +# to zero, then reconstruct the image using the inverse DCT. +# +# J(abs(J)<10) = 0; +# K = idctn(J); +# figure, imshow(I) +# figure, imshow(K,[0 255]) +# +# See also +# -------- +# dctn, idct, dct +# +# -- Damien Garcia -- 2009/04, revised 2009/11 +# website: www.BiomeCardio.com +# +# ---------- +# [Y,W] = IDCTN(X,DIM,W) uses and returns the weights which are used by +# the program. If IDCTN is required for several large arrays of same +# size, the weights can be reused to make the algorithm faster. A typical +# syntax is the following: +# w = []; +# for k = 1:10 +# [y{k},w] = idctn(x{k},[],w); +# end +# The weights (w) are calculated during the first call of IDCTN then +# reused in the next calls. +# ''' +# +# y = np.atleast_1d(y) +# shape0 = y.shape +# +# if axis is None: +# y = y.squeeze() # Working across singleton dimensions is useless +# +# dimy = y.ndim +# if dimy==1: +# y = np.atleast_2d(y) +# y = y.T +# # Some modifications are required if Y is a vector +## if isvector(y): +## if y.shape[0]==1: +## if axis==0: +## return y, None +## elif axis==1: +## axis=0 +## y = y.T +## elif axis==1: +## return y, None +## +# +# +# if w is None: +# w = [0,] * dimy +# for dim in range(dimy): +# if axis is not None and dim!=axis: +# continue +# n = (dimy==1)*y.size + (dimy>1)*shape0[dim] +# #w{dim} = exp(1i*(0:n-1)'*pi/2/n); +# w[dim] = np.exp(1j * np.arange(n) * np.pi / (2 * n)) +# # --- IDCT algorithm --- +# if np.iscomplex(y).any(): +# y = np.complex(idctn(np.real(y),axis,w),idctn(np.imag(y),axis,w)) +# else: +# for dim in range(dimy): +# y = shiftdim(y,1) +# if axis is not None and dim!=axis: +# #y = shiftdim(y, 1) +# continue +# siz = y.shape +# n = siz[-1] +# +# y = y.reshape((-1,n)) * w[dim] +# y[:,0] = y[:,0]/np.sqrt(2) +# y = (np.fft.ifft(y, n=n, axis=1)).real +# y = y * np.sqrt(2*n) +# +# I = np.empty(n,dtype=int) +# I.put(np.r_[0:n:2],np.r_[0:int(n//2)+np.remainder(n,2)]) +# I.put(np.r_[1:n:2],np.r_[n-1:int(n//2)-1:-1]) +# y = y[:,I] +# +# y = y.reshape(siz) +# +# +# y = y.reshape(shape0); +# return y, w @@ -323,33 +570,73 @@ def no_leading_ones(x): return x[first:] def shiftdim(x, n=None): + ''' + Shift dimensions + + Parameters + ---------- + x : array + n : int + + Notes + ----- + Shiftdim is handy for functions that intend to work along the first + non-singleton dimension of the array. + + If n is None returns the array with the same number of elements as X but + with any leading singleton dimensions removed. + + When n is positive, shiftdim shifts the dimensions to the left and wraps + the n leading dimensions to the end. + + When n is negative, shiftdim shifts the dimensions to the right and pads + with singletons. + + See also + -------- + reshape, squeeze + ''' if n is None: - # returns the array B with the same number of - # elements as X but with any leading singleton - # dimensions removed. return x.reshape(no_leading_ones(x.shape)) elif n>=0: - # When n is positive, shiftdim shifts the dimensions - # to the left and wraps the n leading dimensions to the end. return x.transpose(np.roll(range(x.ndim), -n)) else: - # When n is negative, shiftdim shifts the dimensions - # to the right and pads with singletons. return x.reshape((1,)*-n+x.shape) def test_dctn(): - a = np.arange(12).reshape((3,-1)) - #y = dct(a) - #x = idct(y) - #print(y) - #print(x) + a = np.arange(12) #.reshape((3,-1)) + print('a = ', a) + print(' ') + y = dct(a) + x = idct(y) + print('y = dct(a)') + print(y) + print('x = idct(y)') + print(x) + print(' ') - print(a) - yn = dctn(a)[0] - xn = idctn(yn)[0] +# y1 = dct1(a) +# x1 = idct1(y) +# print('y1 = dct1(a)') +# print(y1) +# print('x1 = idct1(y)') +# print(x1) +# print(' ') + yn = dctn(a) + xn = idctn(yn) + print('yn = dctn(a)') print(yn) + print('xn = idctn(yn)') print(xn) + print(' ') + +# yn1 = dctn1(a) +# xn1 = idctn1(yn1) +# print('yn1 = dctn1(a)') +# print(yn1) +# print('xn1 = idctn1(yn)') +# print(xn1) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 3e2d11e..64a3f80 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -15,9 +15,12 @@ from misc import tranproc #, trangood from numpy import pi, sqrt, atleast_2d, exp, newaxis #@UnresolvedImport from scipy import interpolate, linalg, sparse from scipy.special import gamma +from scipy.optimize import brentq from wafo.misc import meshgrid from wafo.wafodata import WafoData +from dctpack import dct + import copy import numpy as np import scipy @@ -1034,7 +1037,7 @@ class Kernel(object): 'Density estimation for statistics and data analysis' Chapman and Hall, pp 31, 103, 175 ''' - def __init__(self, name, fun='hns'): + def __init__(self, name, fun='hisj'): #'hns'): self.kernel = _MKERNEL_DICT[name[:4]] #self.name = self.kernel.__name__.replace('mkernel_', '').title() try: @@ -1363,6 +1366,99 @@ class Kernel(object): h[dim] = h1 #end % for dim loop return h + def hisj(self, data, inc=128): + ''' + HISJ Improved Sheather-Jones estimate of smoothing parameter. + + Unlike many other implementations, this one is immune to problems + caused by multimodal densities with widely separated modes. The + estimation does not deteriorate for multimodal densities, because + it do not assume a parametric model for the data. + + Parameters + ---------- + data - a vector of data from which the density estimate is constructed; + n - the number of mesh points used in the uniform discretization + + Returns + ------- + bandwidth - the optimal bandwidth + + Reference + --------- + Kernel density estimation via diffusion + Z. I. Botev, J. F. Grotowski, and D. P. Kroese (2010) + Annals of Statistics, Volume 38, Number 5, pages 2916-2957. + ''' + A = np.atleast_2d(data) + d, n = A.shape + + # R= int(mkernel(x)^2), mu2= int(x^2*mkernel(x)) + mu2, R, unusedRdd = self.stats() + STEconstant = R / (mu2 ** (2) * n) + + amin = A.min(axis=1) # Find the minimum value of A. + amax = A.max(axis=1) # Find the maximum value of A. + arange = amax - amin # Find the range of A. + + #% xa holds the x 'axis' vector, defining a grid of x values where + #% the k.d. function will be evaluated. + + ax1 = amin - arange / 8.0 + bx1 = amax + arange / 8.0 + + kernel2 = Kernel('gauss') + mu2, R, unusedRdd = kernel2.stats() + STEconstant2 = R / (mu2 ** (2) * n) + + def fixed_point(t, N, I, a2): + ''' this implements the function t-zeta*gamma^[L](t)''' + + prod = np.prod + L = 7 + logI = np.log(I) + f = 2 * pi**(2 * L) * (a2 * exp(L * logI -I * pi ** 2 * t)).sum() + for s in range(L - 1, 1, -1): + K0 = prod(np.r_[1:2*s:2]) / sqrt(2 * pi); + const = (1 + (1. / 2) ** (s + 1. / 2)) / 3 + time = (2 * const * K0 / N / f) ** (2. / (3 + 2 * s)) + f = 2 * pi ** (2 * s) * (a2 * exp(s * logI-I * pi ** 2 * time)).sum() + return t - (2 * N * sqrt(pi) * f) ** (-2. / 5) + + h = np.empty(d) + for dim in range(d): + ax = ax1[dim] + bx = bx1[dim] + xa = np.linspace(ax, bx, inc) + R = bx-ax + + c = gridcount(A[dim], xa) + N = len(set(A[dim])) + a = dct(c/c.sum(), norm=None) + + #% now compute the optimal bandwidth^2 using the referenced method + I = np.asfarray(np.arange(1, inc))**2 + a2 = (a[1:]/2)**2 + fun = lambda t: fixed_point(t, N, I, a2) + # use fzero to solve the equation t=zeta*gamma^[5](t) + try: + t_star = brentq(fun, a=0, b=.1) + except: + t_star = 0.28*N**(-2./5) + + # smooth the discrete cosine transform of initial data using t_star + # a_t = a*exp(-np.arange(inc)**2*pi**2*t_star/2) + # now apply the inverse discrete cosine transform + #density = idct(a_t)/R; + + # take the rescaling of the data into account + bandwidth = sqrt(t_star)*R + + # Kernel other than Gaussian scale bandwidth + h[dim] = bandwidth * (STEconstant / STEconstant2) ** (1.0 / 5) + #end % for dim loop + return h + def hstt(self, data, h0=None, inc=128, maxit=100, releps=0.01, abseps=0.0): '''HSTT Scott-Tapia-Thompson estimate of smoothing parameter. @@ -2461,7 +2557,7 @@ def kde_demo1(): pylab.axis([x.min(), x.max(), 0, 0.5]) def kde_demo2(): - '''Demonstrate the difference between transformation- and ordinary-KDE + '''Demonstrate the difference between transformation- and ordinary-KDE in 1D KDEDEMO2 shows that the transformation KDE is a better estimate for Rayleigh distributed data around 0 than the ordinary KDE. @@ -2469,10 +2565,10 @@ def kde_demo2(): import scipy.stats as st data = st.rayleigh.rvs(scale=1, size=300) - x = np.linspace(1.5e-3, 5, 55); + x = np.linspace(1.5e-2, 5, 55); kde = KDE(data) - f = kde(output='plot', title='Ordinary KDE') + f = kde(output='plot', title='Ordinary KDE (hs=%g)' % kde.hs) pylab.figure(0) f.plot() @@ -2481,7 +2577,7 @@ def kde_demo2(): #plotnorm((data).^(L2)) % gives a straight line => L2 = 0.5 reasonable tkde = TKDE(data, L2=0.5) - ft = tkde(x, output='plot', title='Transformation KDE') + ft = tkde(x, output='plot', title='Transformation KDE (hs=%g)' % tkde.tkde.hs) pylab.figure(1) ft.plot() @@ -2490,7 +2586,7 @@ def kde_demo2(): pylab.figure(0) def kde_demo3(): - '''Demonstrate the difference between and ordinary-KDE + '''Demonstrate the difference between transformation and ordinary-KDE in 2D KDEDEMO3 shows that the transformation KDE is a better estimate for Rayleigh distributed data around 0 than the ordinary KDE. @@ -2521,6 +2617,9 @@ def kde_demo3(): pylab.figure(0) def kde_demo4(hs=None, fast=False): + ''' + + ''' N = 100 #ei = np.random.normal(loc=0, scale=0.075, size=(N,)) ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525, @@ -2568,6 +2667,36 @@ def kde_demo4(hs=None, fast=False): print(kreg.tkde.tkde.inv_hs) print(kreg.tkde.tkde.hs) + +def kde_demo5(N=50): + '''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior + for multimodal distributions + + KDEDEMO5 shows that the improved Sheather-Jones plug-in smoothing is a better + compared to normal reference rules (in this case the hns) + ''' + import scipy.stats as st + + data = np.hstack((st.norm.rvs(loc=5, scale=1, size=(N,)), st.norm.rvs(loc=-5, scale=1, size=(N,)))) + + #x = np.linspace(1.5e-3, 5, 55) + + kde = KDE(data, kernel=Kernel('gauss', 'hns')) + f = kde(output='plot', title='Ordinary KDE', plotflag=1) + + + + kde1 = KDE(data, kernel=Kernel('gauss', 'hisj')) + f1 = kde1(output='plot', label='Ordinary KDE', plotflag=1) + + pylab.figure(0) + f.plot('r', label='hns=%g' % kde.hs) + #pylab.figure(2) + f1.plot('b', label='hisj=%g' % kde1.hs) + x = np.linspace(-4,4) + for loc in [-5,5]: + pylab.plot(x + loc, st.norm.pdf(x, 0, scale=1)/2, 'k:', label='True density') + pylab.legend() def test_docstrings(): import doctest @@ -2575,4 +2704,4 @@ def test_docstrings(): if __name__ == '__main__': #test_docstrings() - kde_demo4() \ No newline at end of file + kde_demo2() \ No newline at end of file diff --git a/pywafo/src/wafo/test/test_kdetools.py b/pywafo/src/wafo/test/test_kdetools.py index 4a603d1..bdd5180 100644 --- a/pywafo/src/wafo/test/test_kdetools.py +++ b/pywafo/src/wafo/test/test_kdetools.py @@ -212,6 +212,16 @@ def test_smooth_params(): >>> gauss.hldpi(data) array([ 0.1732289 , 0.33159097, 0.3107633 ]) + + + >>> gauss.hisj(data) + array([ 0.29400043, 0.74277133, 0.51251583]) + + >>> data = np.array([0.753557920000000, 0.727791940000000, 0.941491690000000, + ... 0.078411190000000, 2.322918870000000, 1.104199950000000, 0.770551140000000, + ... 0.602882730000000, 1.368836350000000, 1.747543260000000, 1.095475610000000, + ... 1.016711330000000, 0.732111430000000, 0.618917190000000, 0.759034870000000, + ... 1.891946900000000, 0.724338080000000, 1.929730940000000, 0.447498380000000, 1.365084520000000]) ''' def test_gridcount_1D(): '''