[SciPy-user] 3D density calculation

[David Huard]

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>:
>
> 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> 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
> > 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 --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20070618/98816140/attachment.html>