Fluorescence Correlation Spectroscopy Simulation

I am simulating Fluorescence Correlation Spectroscopy which basically involves tracking the random motion of particles in a box with periodic boundary conditions and then calculating their intensity using a 3D Gaussian function where the particles are at their maximum intensity at the center of the box. Finally, the intensities of the particles are then summed up for each time step. The problem I am having is manipulating a very large list of lists when calculating the intensity of each particle at each time step.

Timing[Clear["Global*"]
ClearSystemCache[];
n = 10.;
k = 1.38*10.^-23.;
T = 293.;
Eta = 1.*10.^-3.;
Deltat = 500000.;
time = 10.*10.^-6.;
Taud = (Omegar^2./(4.*d));
Omegar = 200.*10.^-9.;
Omegaz = 5.*Omegar;
Io = 1;
boundary = 5.*10.^-6.;
initial = boundary/2.;
step = Sqrt[2.*d*time];
RandomWalk[x_] :=
Accumulate[
Join[{RandomReal[{-initial, initial}, 3]},
RandomVariate[NormalDistribution[0, step], {x, 3}]]];
p = Table[
Mod[RandomWalk[Deltat], boundary, -initial], {i, n}];
particleintensity =
Table[Io*Exp[(-2.*(p[[i, t, 1]]^2. + p[[i, t, 2]]^2.))/Omegar^2. - (2. (p[[i, t,
3]]^2.))/Omegaz^2.], {t, 1, Deltat + 1}, {i, n}];
int = Total[particleintensity, {2}];
ListPlot[int, AxesOrigin -> {0., 0.}, AxesLabel -> {"t", "I"},
PlotRange -> Full, PlotRangeClipping -> False, PlotStyle -> Red]]

{87.2093,"Plot"}


In the example above, I have 10 particles for 500,000 time steps and it takes a considerable amount of time. The majority of the calculation time occurs at the end when calculating the variable particleintensity. I am thinking that there has to be a faster way to do this that utilizes less memory. Any suggestions would be greatly appreciated.

• Not an answer to your question per se ... but ... if you are doing work on Fluorescence Fluctuation Spectroscopy, you might be interested in a paper that uses Mathematica and mathStatica to analyse this topic by: Muller, Joachim D. (2004), Cumulant Analysis in Fluorescence Fluctuation Spectroscopy, Biophysical Journal, Volume 86, 3981–3992 Nov 4, 2012 at 16:21
• What happens if you use particleintensity = Io*Exp[-2 Table[With[{v = p[[i, t]]}, Norm[Take[v, 2]]^2/Omegar^2 + (Last[v]^2)/Omegaz^2], {t, Deltat + 1}, {i, n}]]? Nov 4, 2012 at 16:24
• @J.M. It is actually about 10% slower for me.
– Kane
Nov 4, 2012 at 20:44
• @Kane, have you tried to compile the Table?
– user21
Nov 5, 2012 at 0:32
• @Kane this is rather interesting. Is this molecular dynamics? Nov 5, 2012 at 14:29

There's no need to iterate over i and t to build particleintensity, just do the whole dot product in one go. The Exp moves outside the dot product (so we use it once as a Listable function rather than 5 million times individually).

particleintensity = Io Exp[(p^2).{-2/Omegar^2, -2/Omegar^2, -2/Omegaz^2}];
int = Total[particleintensity];


This runs about 250 times faster than the original Table version.

Update

Since the final result int is obtained by summing over n, there is no need to hold all the data in memory at once. The $n$ elements of p and particleintensity can be computed sequentially, accumulating the values of particleintensity at each step. This will result in lower memory usage:

int = 0;
Do[
p = Mod[RandomWalk[Deltat], boundary, -initial];
particleintensity = Io Exp[(p^2).{-2/Omegar^2, -2/Omegar^2, -2/Omegaz^2}];
int += particleintensity,
{n}]

• This has been working great. However, I have found that I need Deltat to be at least 1000000 and n to be at least 50. Using these parameters, the calculation of particle intensity still takes a considerable amount of time. Realistically, I only need to calculate the particle intensity when particles are close to the origin, say within a distance of Omegaz or Omegar. If I could only calculate particle intensity when particles are close to the origin, I think that would seriously reduce my calculation time. Any thoughts?
– Kane
Nov 18, 2012 at 15:48
• @Kane, for me the slowest part of the calculation now is the generation of the random numbers. Are you still finding the computation of particleintensity to be the bottleneck? Nov 19, 2012 at 13:01
• Yes, when i set n to 50 and Deltat to 10^6 'particleintensity' is the slowest bit. However, this is not the case when n is set to 5 and Deltat is 10^6. I need a way to only calculate particle intensity for particles that are close to the origin.
– Kane
Nov 19, 2012 at 16:05
• @Kane, with n=50 and Deltat=10^6 I get timings of 13s for calculating p and 2.5s for calculating particleintensity, so I can't reproduce your problem. It could be memory related - I'll post an edit with lower memory usage. Nov 19, 2012 at 22:23
• Yes, I suspect that it is a memory issue. I am working with 4GB of memory and a 2GHz Intel Core 2 Duo processor.
– Kane
Nov 20, 2012 at 3:07

I managed to double the speed by simply moving some multiplications and divisions...

AbsoluteTiming[
or2 = -2./Omegar^2.;
oz2 = 2./Omegaz^2.;
p = p^2.;
Io*Table[
Exp[(p[[i, t, 1]] + p[[i, t, 2]])*or2 - p[[i, t, 3]]*oz2],
{t, 1, Deltat + 1}, {i, n}];
]


A dot product is 30 to 40% faster still.

AbsoluteTiming[
or2 = -2./Omegar^2.;
oz2 = 2./Omegaz^2.;
v = {or2, or2, -oz2};
p = p^2.;
Io*Table[Exp[p[[i, t]].v],{t,1,Deltat+1},{i,n}];
]
`
• This is about 50% faster for me too. Thanks.
– Kane
Nov 4, 2012 at 20:57
• That's even better! Thanks.
– Kane
Nov 5, 2012 at 16:03