I'm working on a Monte Carlo code and suddenly I've
found that sorting the results into cells and summing
the contents of each cell is taking far longer than
the actual computation. I'd guess that this type of
sorting problem is well known, and before I spend much
time on it, I'd like to ask for suggestions.
I have large numbers of photon bundles bouncing
around in a slab, and when they emerge they have a
specific intensity, but are also tagged with three
coordinates: the location on the surface where they
exit, and the angle with which they emerge.
So the problem looks like this:
a -- intensity array, size = np (nphotons > 1e5)
a0, a1, a2 -- tags for position & angle, all size np
f0, f1, f2 -- boundaries of the cells into which I want
to sort the the photons.
Thus, b0=. f0 I. a0 will be the index which shows into
which f0 interval each photon in "a" falls. The number
of boundaries might be ~40, so there will be 40^3=64000
cells to which the exiting photons might be assigned.
Here's something that works, but is not J-ish: it's
looping and takes far to long:
b=. >(f0 I. a0);(f1 I. a1);(f2 I. a2)
n0=. >: #f0
n1=. >: #f1
n2=. >: #f2
S=. (n0;n1;n2) Gg b;a
where
Gg=: 4 : 0
'd1 d2 d3'=. x
'bin Iout'=. y
S=. (d1,d2,d3)$0 NB. empty array for results
i1=. i2=. i3=. 0
while. i1<(d1+1) do.
while. i2<(d2+1) do.
while. i3<(d3+1) do.
ii=. i1,i2,i3
idx=. I. */ ii=bin
if. ($idx)>0 do. NB. many bins may be empty
item=. +/ idx {Iout NB. sum all in bin
S=. item (<ii)} S NB. put sum into output array
end.
i3=. >:i3
end.
i3=. 0
i2=. >: i2
end.
i2=. 0
i1=. >:i1
end.
S
)
The idea is to make "np" as large as possible so we can
do the scattering computations efficiently; thus making
an array of size n0*n1*n2*np is out of the question.
I'd guess that experts on this forum would know how to
approach this type of problem. I'd be grateful for any
hints.
Patrick
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
