Looking for suggestions to speed up my code

Question

What i have to do is simple: estimate distances between a huge number (halonumb, roughly 130000) of objects in space then count how many couples of objects have distances in 10 different range. Range's extrema are here notated as binsEX. Range index is n, object index are i and j. I did some test and basic math and it looks like i need something around 25 days in order to do this. That's too much. Any suggestions about how to speed up my code? I thought about using mathematica's parallel computing capabilities but i'm not used to the argument by any means and ended up with actual times longer then those from single kernel computation.

Do[

If[ i < j,

  Do[

   If[DataDistance[i, j] > binsEX[[n]] && 
  DataDistance[i, j] < binsEX[[n + 1]],

     DDbins[[n]] = DDbins[[n]] + 1;
   , {n, 10}
    ]
  ]

 , {i, halonumb}, {j, halonumb}]

Answer 1

Look up BinCounts and DistanceMatrix.

Example

Here are 20,000 points in the unit square. As I understand, you want the histogram of their pairwise distances.

pts = RandomReal[1, {20000, 2}];

DistanceMatrix[pts] would calculate pairwise distances. But the matrix it would produce will not fit in memory (for 20000 it might, but for 100000 it will not). So we need to compute the distance matrix block by block. Let one block be of size blockSize by blockSize.

blockSize = 1000;

parts = Partition[pts, blockSize, blockSize, {1, 1}, {}];

res = ConstantArray[0, 10];
Do[
  res += BinCounts[
    Flatten@DistanceMatrix[p1, p2], {0, 1, 0.1}],
  {p1, parts}, {p2, parts}
] // AbsoluteTiming

res[[1]] -= Length[pts];
res /= 2;

This takes 7.5 seconds on my machine. For your size of data, it would take 5-6 minutes. If you have higher-than-2 dimensional data, it would take longer, but not significantly so.

The last two line of the code correct for counting the zero-distances on the diagonal of the matrix, and for double counting all other distances.

Parallel version

We can speed this up somewhat using parallelization. Here's how:

LaunchKernels[]

DistributeDefinitions[parts]

res = Total@ParallelMap[
   BinCounts[Flatten@DistanceMatrix[parts[[First[#]]], parts[[Last[#]]]], {0, 1, 0.1}] &,
   Tuples[Range@Length[parts], {2}]
  ]; // AbsoluteTiming

res[[1]] -= Length[pts];
res /= 2;

A trick I used here was to distribute parts only once, and then index into it. A more straightforward solution would send a block of data with each evaluation. This would result in a lot of data transfer overhead.

With a block size of 5000 and 130000 points, this took 220 s on my 4-core laptop. It is not a huge speedup over 5-6 minutes (which I have extrapolated, not measured), but it is a useful speedup nevertheless. Once reason for the modest performance increase is that DistanceMatrix already makes use of all CPU cores.

Answer 2

As far as i understood you want to fill a histogram of distances with the number of point pairs which match the binned distance.

You can reduce the neccessary calculations by a huge amount if you discriminate your points to spacial bins and calculate the distances just for the points in the bins and for the bins themselve.

Especially if your bin size is constant you just can set ub spacial bins with a certain size, fill them with your points and then just calculate the distance of the points in the specific spacial bins. Afterwards you calculate the distance between the spacial bins and fill the corresponding bin of your final histogram with the number of points in both spacial bins.