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@gmail.com <mailto:bernhard.voigt@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@tnw.utwente.nl <mailto:c.j.lee@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@scipy.org <mailto:SciPy-user@scipy.org> http://projects.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org <mailto:SciPy-user@scipy.org> http://projects.scipy.org/mailman/listinfo/scipy-user
------------------------------------------------------------------------
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user