Fluorescence Correlation Spectroscopy Simulation

Question

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.;
 radius = 5.*10.^-9.;
 k = 1.38*10.^-23.;
 T = 293.;
 Eta = 1.*10.^-3.;
 d = (k*T)/(6.*\[Pi]*Eta*radius);
 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.

Answer 1

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}]

Answer 2

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}];
]