|
|
|
@ -981,7 +981,17 @@ def accum(accmap, a, func=None, size=None, fill_value=0, dtype=None):
|
|
|
|
|
out[s] = func(vals[s])
|
|
|
|
|
|
|
|
|
|
return out
|
|
|
|
|
def bitget(int_type, offset):
|
|
|
|
|
'''
|
|
|
|
|
Returns the value of the bit at the offset position in int_type.
|
|
|
|
|
|
|
|
|
|
Example
|
|
|
|
|
-------
|
|
|
|
|
>>> bitget(5, np.r_[0:4])
|
|
|
|
|
array([1, 0, 1, 0])
|
|
|
|
|
'''
|
|
|
|
|
mask = (1 << offset)
|
|
|
|
|
return (int_type & mask) != 0
|
|
|
|
|
def gridcount(data, X):
|
|
|
|
|
'''
|
|
|
|
|
GRIDCOUNT D-dimensional histogram using linear binning.
|
|
|
|
@ -1016,11 +1026,13 @@ def gridcount(data, X):
|
|
|
|
|
>>> import numpy as np
|
|
|
|
|
>>> import wafo.kdetools as wk
|
|
|
|
|
>>> import pylab as plb
|
|
|
|
|
>>> N = 500;
|
|
|
|
|
>>> N = 20;
|
|
|
|
|
>>> data = np.random.rayleigh(1,N)
|
|
|
|
|
>>> x = np.linspace(0,max(data)+1,50)
|
|
|
|
|
>>> dx = x[1]-x[0]
|
|
|
|
|
|
|
|
|
|
>>> c = wk.gridcount(data,x)
|
|
|
|
|
|
|
|
|
|
>>> h = plb.plot(x,c,'.') # 1D histogram
|
|
|
|
|
>>> h1 = plb.plot(x,c/dx/N) # 1D probability density plot
|
|
|
|
|
>>> np.trapz(x,c/dx/N)
|
|
|
|
@ -1043,7 +1055,6 @@ def gridcount(data, X):
|
|
|
|
|
if d != d1:
|
|
|
|
|
raise ValueError('Dimension 0 of data and X do not match.')
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
dx = np.diff(x[:, :2], axis=1)
|
|
|
|
|
xlo = x[:, 0]
|
|
|
|
|
xup = x[:, -1]
|
|
|
|
@ -1061,47 +1072,46 @@ def gridcount(data, X):
|
|
|
|
|
abs = np.abs
|
|
|
|
|
if d == 1:
|
|
|
|
|
x.shape = (-1,)
|
|
|
|
|
c = (accum(binx, (x[binx + 1] - data), size=[inc, ]) +
|
|
|
|
|
accum(binx, (data - x[binx]), size=[inc, ])) / w
|
|
|
|
|
c = (accum(binx, (x[binx + 1] - dat), size=[inc, ]) +
|
|
|
|
|
accum(binx, (dat - x[binx]), size=[inc, ])) / w
|
|
|
|
|
elif d == 2:
|
|
|
|
|
b2 = binx[1]
|
|
|
|
|
b1 = binx[0]
|
|
|
|
|
c_ = np.c_
|
|
|
|
|
stk = np.vstack
|
|
|
|
|
c = (accum(c_[b1, b2] , abs(np.prod(stk([X[0, b1 + 1], X[1, b2 + 1]]) - data, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 + 1, b2] , abs(np.prod(stk([X[0, b1], X[1, b2 + 1]]) - data, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 , b2 + 1], abs(np.prod(stk([X[0, b1 + 1], X[1, b2]]) - data, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 + 1, b2 + 1], abs(np.prod(stk([X[0, b1], X[1, b2]]) - data, axis=0)), size=[inc, inc])) / w
|
|
|
|
|
c = (accum(c_[b1, b2] , abs(np.prod(stk([X[0, b1 + 1], X[1, b2 + 1]]) - dat, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 + 1, b2] , abs(np.prod(stk([X[0, b1], X[1, b2 + 1]]) - dat, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 , b2 + 1], abs(np.prod(stk([X[0, b1 + 1], X[1, b2]]) - dat, axis=0)), size=[inc, inc]) +
|
|
|
|
|
accum(c_[b1 + 1, b2 + 1], abs(np.prod(stk([X[0, b1], X[1, b2]]) - dat, axis=0)), size=[inc, inc])) / w
|
|
|
|
|
|
|
|
|
|
else: # % d>2
|
|
|
|
|
raise ValueError('Not implemented for d>2')
|
|
|
|
|
|
|
|
|
|
Nc = csiz.prod()
|
|
|
|
|
c = np.zeros((Nc, 1))
|
|
|
|
|
c = np.zeros((Nc,))
|
|
|
|
|
|
|
|
|
|
fact2 = inc * np.arange(d)
|
|
|
|
|
fact1 = csiz.cumprod() / inc
|
|
|
|
|
fact2 = np.asarray(np.reshape(inc * np.arange(d), (d, -1)), dtype=int)
|
|
|
|
|
fact1 = np.asarray(np.reshape(csiz.cumprod() / inc, (d, -1)), dtype=int)
|
|
|
|
|
#fact1 = fact1(ones(n,1),:);
|
|
|
|
|
# for ir in xrange(2**(d-1)):
|
|
|
|
|
# bt0[:,:,1] = bitget(ir,1:d)
|
|
|
|
|
# bt0[:,:,2] = 1-bt0[:,:,1]
|
|
|
|
|
# for ix in range(2):
|
|
|
|
|
# one = mod(ix,2)+1;
|
|
|
|
|
# two = mod(ix+1,2)+1;
|
|
|
|
|
# # Convert to linear index (faster than sub2ind)
|
|
|
|
|
# b1 = sum((binx + bt0(ones(n,1),:,one)-1).*fact1,2)+1; #%linear index to c
|
|
|
|
|
# bt2 = bt0(:,:,two) + fact2;
|
|
|
|
|
# b2 = binx + bt2(ones(n,1),:); #% linear index to X
|
|
|
|
|
#
|
|
|
|
|
# c = c + accum(b1,abs(prod(X(b2)-data,2)),[Nc,1]);
|
|
|
|
|
# #c = c + accum([b1,ones(n,1)],abs(prod(X(b2)-data,2)),[Nc,1]);
|
|
|
|
|
# #[len,bin,val] = bincount(b1,abs(prod(X(b2)-data,2)));
|
|
|
|
|
# #c(bin) = c(bin)+val;
|
|
|
|
|
#
|
|
|
|
|
# #end
|
|
|
|
|
# #end
|
|
|
|
|
# c = reshape(c/w,csiz);
|
|
|
|
|
#end
|
|
|
|
|
if d == 2: #% make sure c is stored in the same way as meshgrid
|
|
|
|
|
bt0 = [0, 0]
|
|
|
|
|
X1 = X.ravel()
|
|
|
|
|
for ir in xrange(2 ** (d - 1)):
|
|
|
|
|
bt0[0] = np.reshape(bitget(ir, np.arange(d)), (d, -1))
|
|
|
|
|
bt0[1] = 1 - bt0[0]
|
|
|
|
|
for ix in xrange(2):
|
|
|
|
|
one = np.mod(ix, 2)
|
|
|
|
|
two = np.mod(ix + 1, 2)
|
|
|
|
|
# Convert to linear index
|
|
|
|
|
b1 = np.sum((binx + bt0[one]) * fact1, axis=0) #linear index to c
|
|
|
|
|
bt2 = bt0[two] + fact2
|
|
|
|
|
b2 = binx + bt2 # linear index to X
|
|
|
|
|
c += accum(b1, abs(np.prod(X1[b2] - dat, axis=0)), size=(Nc,))
|
|
|
|
|
|
|
|
|
|
c = np.reshape(c / w, csiz)
|
|
|
|
|
# TODO: check that the flipping of axis is correct
|
|
|
|
|
T = range(d); T[-2],T[-1] = T[-1], T[-2]
|
|
|
|
|
c = c.transpose(*T)
|
|
|
|
|
|
|
|
|
|
if d == 2: # make sure c is stored in the same way as meshgrid
|
|
|
|
|
c = c.T
|
|
|
|
|
elif d == 3:
|
|
|
|
|
c = c.transpose(1, 0, 2)
|
|
|
|
@ -1130,6 +1140,7 @@ def test_gridcount():
|
|
|
|
|
data = np.random.rayleigh(1, size=(2, N))
|
|
|
|
|
x = np.linspace(0, max(data.ravel()) + 1, 10)
|
|
|
|
|
X = np.vstack((x, x))
|
|
|
|
|
|
|
|
|
|
dx = x[1] - x[0]
|
|
|
|
|
c = wk.gridcount(data, X)
|
|
|
|
|
h = plb.contourf(x, x, c)
|
|
|
|
|