[SciPy-user] 3D density calculation
Chris Lee
c.j.lee@tnw.utwente...
Thu Jun 28 12:14:02 CDT 2007
That would be the problem. I just tested it on a machine with an up
to date version of numpy and the problem went away.
Out of curiosity (and perhaps practically since I have no power to
upgrade the other machine I was using) is the buggy output always a
transpose of the correct output?
Cheers
Chris
On Jun 28, 2007, at 5:31 PM, David Huard wrote:
> Hi Chris,
>
> It's a bug alright, but I believe it has been fixed. What version
> of numpy are you using? If you're not using the latest release,
> please try it out and see if it still bugs. The bug had to do with
> the reordering of the flattened array into a N-D array, which is
> bin length dependent.
>
> Cheers,
> David
>
> 2007/6/28, Chris Lee <c.j.lee@tnw.utwente.nl>:
> 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@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
> >
>
> _______________________________________________
> SciPy-user mailing list
> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-user/attachments/20070628/e3b9f393/attachment.html
More information about the SciPy-user
mailing list