[SciPy-user] 3D density calculation

[Chris Lee]

Thank you both for your reply.  Sorry for my lack of response I was at a 
conference and the webmail interface only gave me read only access :(

histogramdd looks like it might do it.  I had thought about consecutive 
1D histograms but realised that it would not work.

I had thought that I couldn't use a fixed grid because I require a large 
total volume and also to be able to resolve small distances.  However, 
the David's idea looks pretty good so I will try that as well.

Thank you both for your help.

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/20070625/98a53efa/attachment.vcf>