[SciPy-user] 3D density calculation

[Chris Lee]

Hi David,

histrogramdd does exactly what I want and seems to be fast enough.  
However, I think it may have an irregular bug

if I have some data with a shape (x,4) where x varies from call to call 
and I specify 4 bin numbers bins=(binx, biny, binz, bint) and then call

hist, edges = histrogramdd(data, bins=(binx, biny, binz, bint)))

then most of the time hist will have the shape (binx, biny, binz, bint) 
but sometimes it will return a different shape e.g., (bint, binz, biny, 
binx) and the edges array order will be different again e.g., (x, y, t, 
z). 

Here is the chunk of code that generates this (with appropriate an 
appropriate non null data of course)

densityHistogram, edges = numpy.histogramdd(data, 
bins=(xBins,yBins,zBins,tBins))
print densityHistogram.shape
SHGPhotons = n.asarray(densityHistogram*0.01, n.int32)
totalPhotons = SHGPhotons.sum()
print SHGPhotons.shape
if totalPhotons>0:
    idxSHG = n.nonzero(SHGPhotons)
    #print idxSHG
    xEdges = edges[0]
    yEdges = edges[1]
    zEdges = edges[2]
    tEdges = edges[3]
    #print zEdges
    print xEdges.shape, yEdges.shape, zEdges.shape, tEdges.shape

Here is a typical (non error generating output)
20 20 1 1  <- number of bins in order
(20, 20, 1, 1) <- shape of histogram
(20, 20, 1, 1) <- shape of an array generated from the histrogram
(21,) (21,) (2,) (2,) <- shape of the inidividual edge arrays in order

here is the version that puts out an error
21 20 2 3 <- bins in order
(3, 2, 21, 20) <- histogram is backwards
(3, 2, 21, 20) <- as is the array generated from it
(22,) (21,) (3,) (4,) <- edge arrays are in order though

So, is this an error, or am I assuming too much about how histogramdd 
operates?

I will start checking the order in the program, but this will never be 
perfect since the number of bins could be the same for multiple axis.

Thanks for any advice
Cheers
Chris



David Huard wrote:
> Hi Chris,
>
> Have you tried numpy.histogramdd ? If its still too slow, I have a 
> fortran implementation on the back burner. I could try to finish it 
> quickly and send you a preliminary version.
>
> Other thought: the kernel density estimator scipy.stats.gaussian_kde
>
> David
>
> 2007/6/17, Bernhard Voigt <bernhard.voigt at gmail.com 
> <mailto:bernhard.voigt at gmail.com>>:
>
>     Hi Chris!
>
>     you could try a grid of unit cells that cover your phase space
>     (x,y,z,t). Count the number of photons per unit cell of your
>     initial configuration and track photons leaving and entering a
>     particular cell. A dictionary with a tuple of x,y,z,t coordinates
>     obtained from integer division of the x,y,z,t coordinates could
>     serve as keys.
>
>     Example for 2-D:
>
>     from numpy import *
>     # phase space in x,y
>     x = arange(-100,100.1,.1)
>     y = arange(-100,100.1,.1)
>     # cell dimension in both dimensions the same
>     GRID_WIDTH=7.5
>
>     # computes the grid key from x,y coordinates
>     def gridKey(x,y):
>         '''return the a tuple of x,y integer divided by GRID_WIDHT'''
>         return (int(x // GRID_WIDTH), int(y // GRID_WIDTH))
>
>     # setup your grid dictionary
>     gridLowX, gridHighX = gridKey(min(x), max(x))
>     gridLowY, gridHighY = gridKey(min(y), max(y))
>     keys = [(i,j) for i in xrange(gridLowX, gridHighX + 1) \
>                 for j in xrange(gridLowY, gridHighY + 1)]
>     grid = dict().fromkeys(keys, 0)
>
>     # random photons
>     photons = random.uniform(-100.,100., (100000,2))
>
>     # count photons in each grid cell
>     for p in photons:
>         grid[gridKey(*p)] += 1
>
>     #########################################
>     # in your simulation you have to keep track of where your photons
>     # are going to...
>     # (the code below won't run, it's just an example)
>     #########################################
>     oldKey = gridKey(photon)
>     propagate(photon)  # changes x,y coordinates of photon
>     newKey = gridKey(photon)
>     if oldKey != newKey:
>         grid[oldKey] -= 1
>         grid[newKey] += 1
>
>     I hope this helps! Bernhard
>
>
>     On 6/15/07, * Chris Lee* < c.j.lee at tnw.utwente.nl
>     <mailto:c.j.lee at tnw.utwente.nl>> wrote:
>
>         Hi everyone,
>
>         I was hoping this list could point me in the direction of a more
>         efficient solution to a problem I have.
>
>         I have 4 vectors: x, y, z, and t that are about 1 million in
>         length
>         that describe the positions of photons.  As my simulation
>         progresses
>         it updates the positions so x, y, z, and t change by an
>         unknown (and
>         unknowable) amount every update.
>
>         This worked very well for its original purpose but now I need to
>         calculate the photon density change over time.  Currently
>         after each
>         update, I iterate over time slices, x slices, and y slices and
>         then
>         make an histogram of z which I then stitch together to create a
>         density.  However, this becomes very slow as the photons
>         spread out
>         in space and time.
>
>         Does anyone know how to take such a large vector set and return a
>         density efficiently?
>
>
>         _______________________________________________
>         SciPy-user mailing list
>         SciPy-user at scipy.org <mailto:SciPy-user at scipy.org>
>         http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
>
>
>     _______________________________________________
>     SciPy-user mailing list
>     SciPy-user at scipy.org <mailto:SciPy-user at scipy.org>
>     http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>   
-------------- next part --------------
A non-text attachment was scrubbed...
Name: c.j.lee.vcf
Type: text/x-vcard
Size: 174 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20070628/21287dd4/attachment.vcf>