3
$\begingroup$

I have a 3x3 matrix (252^3) of data (densities) and I want to compute a correlation xi from nearest neighbours as well as a sum involving a check whether the density is in a certain bin. The computation works in principle (so for subsamples, when step size is about 10) but is painfully slow. I am not an expert on handling data with Mathematica, so I am pretty sure there is some solution that works way better than my first try. For the correlation function I think that basically need an analogue of AbsoluteCorrelationFunction for my matrix.

xi=
   1/3*Sum[data[[i,j,k]] (data[[i + 1,j,k]] + 
          data[[i,j + 1,k]] + 
          data[[i,j,k + 1]]), {i, 1, 251, 1}, {j, 1, 
        251, 1}, {k, 1, 251, 1}]/26^3 - 1;

For the following quantity I am not sure if there is another faster Mathematica routine. I need to check whether the density is in a bin of width deltarho around rho and if yes add the values of the neighbouring densities

brho[rho_]:=(Sum[If[Abs[data[[i,j,k]] - rho] < 
      deltarho/2, (data[[i + 1,j,k]] + 
       data[[i,j + 1,k]] + 
       data[[i,j,k + 1]]), 0], {i, 1, 251, 1}, {j, 1, 
     251, 1}, {k, 1, 251, 1}]/
   Sum[Boole[
     Abs[data[[i,j,k]] - rho] < deltarho/2], {i, 1, 
     251, 1}, {j, 1, 251, 1}, {k, 1, 251, 1}] - 1)

I tried replacing Sum with ParallelSum but at least for the subsample it makes things slower.

Are there any ideas how to speed this up?

Thanks a lot in advance!

edit: improved readability & corrected typo, thanks to comment of Jason

$\endgroup$
  • $\begingroup$ Just first impressions - to improve readability in the code for yourself and others, you can replace things like data[[i]][[j]][[k]] with data[[i,j,k]] and , {i, 1, 251, 1}, {j, 1, 251, 1}, {k, 1, 251, 1}] with ,{i, 251}, {j, 251}, {k, 251} but this is just a style point and doesn't address your main question. $\endgroup$ – Jason B. Apr 15 '16 at 9:35
  • $\begingroup$ I tried running your definition of xi on a random matrix dims = {251, 251, 251}; data = RandomReal[1, dims]; and it takes a long time. You could get the same thing virtually instantly by a judicious use of Total and RotateLeft and Drop $\endgroup$ – Jason B. Apr 15 '16 at 9:52
  • $\begingroup$ Also, Can I assume that {i, 1, 251, 10} is a typo in the brho code? It's the only 10 in there so I assume it should be a 1 $\endgroup$ – Jason B. Apr 15 '16 at 10:06
5
$\begingroup$

Well, basically your code is designed with a less efficient algorithm.

Bear in mind that Mathematica generally treats a matrix the "same" as a number, so don't spend time on computation conducted at the number level.

Then we can convert your code to be running at the matrix level, where the speed is much increased.

len = 20;
data = RandomReal[{1, 10}, {len, len, len}];
xi=1/3*Sum[data[[i,j,k]] (data[[i+1,j,k]]+data[[i,j+1,k]]+data[[i,j,k+1]]),{i,len-1},{j,len-1},{k,len-1}]/26^3-1;//AbsoluteTiming
(*{0.030215, Null}*)
xxi=Total@Flatten@(data[[1;;len-1,1;;len-1,1;;len-1]](data[[2;;len,1;;len-1,1;;len-1]]+data[[1;;len-1,2;;len,1;;len-1]]+data[[1;;len-1,1;;len-1,2;;len]]))/3/26^3-1;//AbsoluteTiming
(*{0.00039, Null}*)
xi == xxi
(*True*)

Now we increase len to 252 as required by your scenario:

len = 252;
data = RandomReal[{1, 10}, {len, len, len}];
xxi=Total@Flatten@(data[[1;;len-1,1;;len-1,1;;len-1]](data[[2;;len,1;;len-1,1;;len-1]]+data[[1;;len-1,2;;len,1;;len-1]]+data[[1;;len-1,1;;len-1,2;;len]]))/3/26^3-1;//AbsoluteTiming
(*{0.653753, Null}*)

See, the performance is not bad.

I think now you can improve the second part of your problem with the idea above.

Edit: As for the second part of your problem, there is no need to perform testing over each element: a holistic testing works as well:

Clear[brho]
brho[data_List,rho_,deltarho_]:=Module[{c,sum},
c=UnitStep[deltarho/2-Abs[data[[1;;-2,1;;-2,1;;-2]]-rho]];
sum=data[[2;;-1,1;;-2,1;;-2]]+data[[1;;-2,2;;-1,1;;-2]]+data[[1;;-2,1;;-2,2;;-1]];
Total@Flatten[c*sum]/Total@Flatten[c]
] 

Give it a try:

brho[data,10,3]

which will give the answer.

$\endgroup$
  • $\begingroup$ Thanks so much, you are great guys! This runs in 2s now :-) I hope I will be able to generalize this ideas also to the other functions I want to compute and do the error estimation by using subsamples. $\endgroup$ – Masterplan Apr 16 '16 at 11:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.