Re: Fast calculation of pair correlation function
- To: mathgroup at smc.vnet.net
- Subject: [mg98635] Re: [mg98626] Fast calculation of pair correlation function
- From: Sseziwa Mukasa <mukasa at jeol.com>
- Date: Wed, 15 Apr 2009 04:57:44 -0400 (EDT)
- References: <200904141019.GAA07830@smc.vnet.net>
On Apr 14, 2009, at 6:19 AM, markus wrote:
> Hi,
>
> I have a quite large list (length ~ 10,000) of coordinates (positions
> of particles from a simulation), e.g. data={{x1,y1,z1},
> {x2,y2,z2},...}, and I am looking for fast way to calculate the "pair
> correlation" or "radial distribution function", which tells me, for
> each value of r, how many particles are separated by a distance r (+-
> epsilon).
> Usually, this can be calculated by something like:
>
> Do[Do[histogr[[ Ceiling[Norm[ data[[i]]-data[[j]] ] / delta] ]], {j, i
> +1, Length[data]}], {i, 1, Length[data]-1}]
>
> where histogr is just a list whose index corresponds to the distance
> r.
>
> Unfortunately I have found that for lists of length >1,000, this way
> of calculating is very slow (it can take ~minutes), compared to the
> equivalent C code, which proceeds in a few seconds. I have also tried
> to "Compile" the function, but the speed does not increase...
>
> Does anybody know some fast way of calculating the pair distribution
> function in Mathematica?
The code you provide doesn't compute the count, did you want to
increment the value of histogr[[i]]?
At any rate it seems faster to compute work recursively along the
diagonals of the matrix like this:
In[172]:=
data=RandomReal[{-1.,1.},{1000,3}];
histogr=Range[1,35];
delta=0.1;
Timing[Do[Do[histogr[[Ceiling[Norm[data[[i]]-data[[j]]]/delta]]],{j,i
+1,Length[data]}],{i,1,Length[data]-1}];]
Block[{reclim=$RecursionLimit,ret},$RecursionLimit=Infinity;ret=Timing
[{histogr[[#]]&/@Ceiling[Norm[#1-#2]/delta],Sequence@@If[Length[#1]
>1,#0[Most[#1],Rest[#2]],{}]}&[Most[data],Rest[data]];];
$RecursionLimit=reclim;ret]
Out[175]= {7.35285,Null}
Out[176]= {0.25028,Null}
You can also significantly speed up your code by reducing the number
of times the same object is indexed and using Outer to thread the
subtraction over your list, this method returns the data in the same
order as your original function:
In[219]:=
data=RandomReal[{-1.,1.},{1000,3}];
histogr=Range[1,35];
delta=0.1;
Timing[Do[Do[histogr[[Ceiling[Norm[data[[i]]-data[[j]]]/delta]]],{j,i
+1,Length[data]}],{i,1,Length[data]-1}];]
Timing[Do[Block[{x={data[[i]]}},histogr[[#]]&/@Ceiling[Norm[Flatten
[Outer[Subtract,x,data[[Range[i+1,Length[data]]]],1],1]]/delta]],{i,
1,Length[data]-1}];]
Out[222]= {7.30233,Null}
Out[223]= {0.807352,Null}
- References:
- Fast calculation of pair correlation function
- From: markus <markusg.phys@googlemail.com>
- Fast calculation of pair correlation function