[SciPy-user] 3D density calculation
David Huard
david.huard at gmail.com
Thu Jun 28 13:28:36 EDT 2007
2007/6/28, Chris Lee <c.j.lee at tnw.utwente.nl>:
>
> 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.
>
Good.
> 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?
>
No, it's much worse than that. It's like ordering [1,2,3,4,5,6] to a (2x3)
matrix or a (3x2) matrix
[[1,2,3]
[4,5,6]]
[[1,2],
[3,4],
[5,6]]
One is not the transpose of the other. Sorry about that. You could cut and
paste the fixed code and import histogramdd from a local module instead of
from numpy.
David
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 at 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 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
> > >
> >
> > _______________________________________________
> > 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
>
>
>
> _______________________________________________
> 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/20070628/8341f699/attachment.html>
More information about the SciPy-User
mailing list