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

Added improved Sheater-Jones plugin estimate of the smoothing parameter

Browse Source 
Replaced dct and idct with a call to the ones in scipy.fftpack.
Added n-dimensional  version dctn and idctn to dctpack.py
master
Per.Andreas.Brodtkorb 13 years ago
parent 676d30b822
commit 214d467939
3 changed files with 700 additions and 274 deletions
Show Stats Download Patch File Download Diff File

785
pywafo/src/wafo/dctpack.py
Unescape Escape View File

@ -1,112 +1,150 @@
import numpy as np import numpy as np
from scipy.fftpack import dct as _dct
from scipy.fftpack import idct as _idct
__all__ = ['dct', 'idct', 'dctn', 'idctn'] __all__ = ['dct', 'idct', 'dctn', 'idctn']
def dct(x, n=None):
"""
Discrete Cosine Transform
N-1 def dct(x, type=2, n=None, axis=-1, norm='ortho'):
y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. '''
n=0 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.
Examples See Also
-------- --------
>>> import numpy as np idct
>>> 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 Notes
--------- -----
http://en.wikipedia.org/wiki/Discrete_cosine_transform For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/ MATLAB ``dct(x)``.
"""
fft = np.fft.fft
x = np.atleast_1d(x)
if n is None: There are theoretically 8 types of the DCT, only the first 3 types are
n = x.shape[-1] implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the'
Inverse DCT generally refers to DCT type 3.
if x.shape[-1] < n: type I
n_shape = x.shape[:-1] + (n - x.shape[-1],) ~~~~~~
xx = np.hstack((x, np.zeros(n_shape))) There are several definitions of the DCT-I; we use the following
else: (for ``norm=None``)::
xx = x[..., :n]
real_x = np.all(np.isreal(xx)) N-2
if (real_x and (np.remainder(n, 2) == 0)): y[k] = x[0] + (-1)**k x[N-1] + 2 * sum x[n]*cos(pi*k*n/(N-1))
xp = 2 * fft(np.hstack((xx[..., ::2], xx[..., ::-2]))) n=1
else:
xp = fft(np.hstack((xx, xx[..., ::-1])))
xp = xp[..., :n]
w = np.exp(-1j * np.arange(n) * np.pi / (2 * n)) Only None is supported as normalization mode for DCT-I. Note also that the
DCT-I is only supported for input size > 1
y = xp * w type II
~~~~~~~
There are several definitions of the DCT-II; we use the following
(for ``norm=None``)::
if real_x:
return y.real
else:
return y
def idct(x, n=None):
"""
Inverse Discrete Cosine Transform
N-1 N-1
x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N. y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
n=0 n=0
w(0) = 1/2 If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
w(n) = 1 for n>0
Examples f = sqrt(1/(4*N)) if k = 0,
-------- f = sqrt(1/(2*N)) otherwise.
>>> 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 Which makes the corresponding matrix of coefficients orthonormal
--------- (``OO' = Id``).
http://en.wikipedia.org/wiki/Discrete_cosine_transform
http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
"""
ifft = np.fft.ifft type III
x = np.atleast_1d(x) ~~~~~~~~
if n is None: There are several definitions, we use the following
n = x.shape[-1] (for ``norm=None``)::
w = np.exp(1j * np.arange(n) * np.pi / (2 * n)) N-1
y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N.
n=1
if x.shape[-1] < n: or, for ``norm='ortho'`` and 0 <= k < N::
n_shape = x.shape[:-1] + (n - x.shape[-1],)
xx = np.hstack((x, np.zeros(n_shape))) * w N-1
else: y[k] = x[0] / sqrt(N) + sqrt(1/N) * sum x[n]*cos(pi*(k+0.5)*n/N)
xx = x[..., :n] * w n=1
real_x = np.all(np.isreal(x)) The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
if (real_x and (np.remainder(n, 2) == 0)): to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
xx[..., 0] = xx[..., 0] * 0.5 the orthonormalized DCT-II.
yp = ifft(xx)
y = np.zeros(xx.shape, dtype=complex) References
y[..., ::2] = yp[..., :n / 2] ----------
y[..., ::-2] = yp[..., n / 2::]
else: http://en.wikipedia.org/wiki/Discrete_cosine_transform
yp = ifft(np.hstack((xx, np.zeros_like(xx[..., 0]), np.conj(xx[..., :0:-1]))))
y = yp[..., :n]
if real_x: 'A Fast Cosine Transform in One and Two Dimensions', by J. Makhoul, `IEEE
return y.real 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: else:
return y 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.
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.
def dctn(y, axis=None, w=None): 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 _idct(farr(x), type, n, axis, norm)
def dctn(x, type=2, axis=None, norm='ortho'):
''' '''
DCTN N-D discrete cosine transform. DCTN N-D discrete cosine transform.
@ -146,171 +184,380 @@ def dctn(y, axis=None, w=None):
idctn, dct, idct idctn, dct, idct
''' '''
y = np.atleast_1d(y) y = np.atleast_1d(x)
shape0 = y.shape shape0 = y.shape
if axis is None: if axis is None:
y = y.squeeze() # Working across singleton dimensions is useless y = y.squeeze() # Working across singleton dimensions is useless
dimy = y.ndim ndim = y.ndim
if dimy==1: isvector = max(shape0)==y.size
y = np.atleast_2d(y) if isvector:
y = y.T if ndim==1:
# Some modifications are required if Y is a vector y = np.atleast_2d(y)
# if isvector(y): y = y.T
# if y.shape[0]==1: elif y.shape[0]==1:
# if axis==0: if axis==0:
# return y, None return x
# elif axis==1: elif axis==1:
# axis=0 axis=0
# y = y.T y = y.T
# elif axis==1: elif axis==1:
# return y, None return y
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(): 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: else:
for dim in range(dimy): y = np.asfarray(y)
y = shiftdim(y,1) 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: if axis is not None and dim!=axis:
y = shiftdim(y, 1)
continue continue
siz = y.shape y = _dct(y, type, norm=norm)
n = siz[-1] return y.reshape(shape0)
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): def idctn(x, type=2, axis=None, norm='ortho'):
''' y = np.atleast_1d(x)
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: <a
href="matlab:web('http://www.biomecardio.com')">www.BiomeCardio.com</a>
----------
[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 shape0 = y.shape
if axis is None: if axis is None:
y = y.squeeze() # Working across singleton dimensions is useless y = y.squeeze() # Working across singleton dimensions is useless
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
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(): 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: else:
for dim in range(dimy): y = np.asfarray(y)
y = shiftdim(y,1) 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: if axis is not None and dim!=axis:
#y = shiftdim(y, 1)
continue continue
siz = y.shape y = _idct(y, type, norm=norm)
n = siz[-1] return y.reshape(shape0)
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) #def dct(x, n=None):
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]) # Discrete Cosine Transform
y = y[:,I] #
# N-1
y = y.reshape(siz) # y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
# n=0
#
y = y.reshape(shape0); # Examples
return y, w # --------
# >>> 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: <a
# href="matlab:web('http://www.biomecardio.com')">www.BiomeCardio.com</a>
#
# ----------
# [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:] return x[first:]
def shiftdim(x, n=None): 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: 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)) return x.reshape(no_leading_ones(x.shape))
elif n>=0: 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)) return x.transpose(np.roll(range(x.ndim), -n))
else: else:
# When n is negative, shiftdim shifts the dimensions
# to the right and pads with singletons.
return x.reshape((1,)*-n+x.shape) return x.reshape((1,)*-n+x.shape)
def test_dctn(): def test_dctn():
a = np.arange(12).reshape((3,-1)) a = np.arange(12) #.reshape((3,-1))
#y = dct(a) print('a = ', a)
#x = idct(y) print(' ')
#print(y) y = dct(a)
#print(x) x = idct(y)
print('y = dct(a)')
print(a) print(y)
yn = dctn(a)[0] print('x = idct(y)')
xn = idctn(yn)[0] print(x)
print(' ')
# 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(yn)
print('xn = idctn(yn)')
print(xn) print(xn)
print(' ')
# yn1 = dctn1(a)
# xn1 = idctn1(yn1)
# print('yn1 = dctn1(a)')
# print(yn1)
# print('xn1 = idctn1(yn)')
# print(xn1)

143
pywafo/src/wafo/kdetools.py
Unescape Escape View File

@ -15,9 +15,12 @@ from misc import tranproc #, trangood
from numpy import pi, sqrt, atleast_2d, exp, newaxis #@UnresolvedImport from numpy import pi, sqrt, atleast_2d, exp, newaxis #@UnresolvedImport
from scipy import interpolate, linalg, sparse from scipy import interpolate, linalg, sparse
from scipy.special import gamma from scipy.special import gamma
from scipy.optimize import brentq
from wafo.misc import meshgrid from wafo.misc import meshgrid
from wafo.wafodata import WafoData from wafo.wafodata import WafoData
from dctpack import dct
import copy import copy
import numpy as np import numpy as np
import scipy import scipy
@ -1034,7 +1037,7 @@ class Kernel(object):
'Density estimation for statistics and data analysis' 'Density estimation for statistics and data analysis'
Chapman and Hall, pp 31, 103, 175 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.kernel = _MKERNEL_DICT[name[:4]]
#self.name = self.kernel.__name__.replace('mkernel_', '').title() #self.name = self.kernel.__name__.replace('mkernel_', '').title()
try: try:
@ -1363,6 +1366,99 @@ class Kernel(object):
h[dim] = h1 h[dim] = h1
#end % for dim loop #end % for dim loop
return h 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): def hstt(self, data, h0=None, inc=128, maxit=100, releps=0.01, abseps=0.0):
'''HSTT Scott-Tapia-Thompson estimate of smoothing parameter. '''HSTT Scott-Tapia-Thompson estimate of smoothing parameter.
@ -2461,7 +2557,7 @@ def kde_demo1():
pylab.axis([x.min(), x.max(), 0, 0.5]) pylab.axis([x.min(), x.max(), 0, 0.5])
def kde_demo2(): 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 KDEDEMO2 shows that the transformation KDE is a better estimate for
Rayleigh distributed data around 0 than the ordinary KDE. Rayleigh distributed data around 0 than the ordinary KDE.
@ -2469,10 +2565,10 @@ def kde_demo2():
import scipy.stats as st import scipy.stats as st
data = st.rayleigh.rvs(scale=1, size=300) 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) kde = KDE(data)
f = kde(output='plot', title='Ordinary KDE') f = kde(output='plot', title='Ordinary KDE (hs=%g)' % kde.hs)
pylab.figure(0) pylab.figure(0)
f.plot() f.plot()
@ -2481,7 +2577,7 @@ def kde_demo2():
#plotnorm((data).^(L2)) % gives a straight line => L2 = 0.5 reasonable #plotnorm((data).^(L2)) % gives a straight line => L2 = 0.5 reasonable
tkde = TKDE(data, L2=0.5) 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) pylab.figure(1)
ft.plot() ft.plot()
@ -2490,7 +2586,7 @@ def kde_demo2():
pylab.figure(0) pylab.figure(0)
def kde_demo3(): 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 KDEDEMO3 shows that the transformation KDE is a better estimate for
Rayleigh distributed data around 0 than the ordinary KDE. Rayleigh distributed data around 0 than the ordinary KDE.
@ -2521,6 +2617,9 @@ def kde_demo3():
pylab.figure(0) pylab.figure(0)
def kde_demo4(hs=None, fast=False): def kde_demo4(hs=None, fast=False):
'''
'''
N = 100 N = 100
#ei = np.random.normal(loc=0, scale=0.075, size=(N,)) #ei = np.random.normal(loc=0, scale=0.075, size=(N,))
ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525, ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525,
@ -2569,10 +2668,40 @@ def kde_demo4(hs=None, fast=False):
print(kreg.tkde.tkde.inv_hs) print(kreg.tkde.tkde.inv_hs)
print(kreg.tkde.tkde.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(): def test_docstrings():
import doctest import doctest
doctest.testmod() doctest.testmod()
if __name__ == '__main__': if __name__ == '__main__':
#test_docstrings() #test_docstrings()
kde_demo4() kde_demo2()

10
pywafo/src/wafo/test/test_kdetools.py
Unescape Escape View File

@ -212,6 +212,16 @@ def test_smooth_params():
>>> gauss.hldpi(data) >>> gauss.hldpi(data)
array([ 0.1732289 , 0.33159097, 0.3107633 ]) 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(): def test_gridcount_1D():
''' '''

Write Preview
Loading…
Cancel
Save