# Computing closest set of points to each point in large set — running out of memory (arrays unpacking)

I have a large set of points in 3D and I'm trying to identify all the points that lie within a certain distance of each point. Then using this data store the vector between the pairs of points.

I have about 1M points and this results in millions of vectors.

It seems that when I Map Nearest over my set of points it is unpacking the array (because it becomes ragged?) and this appears to take up all my memory. At the end I Flatten the array to a specific depth and can manually pack it again and I estimate it should fit in memory.

closepointsvectors finds all the points within 1 unit of point x and then returns the vector between the point x and the "close" points. I then Map this over all the points in the dataset. This generates many many vectors.

Here is a minimal (non)working example.

data = RandomReal[10, {10^6, 3}];
nearestfunction = Nearest[data];
closepointsvectors[x_] := (x - #) & /@ nearestfunction[x, {Infinity, 1}]

(*lets try on a small subset of the data*)
smalltest = Flatten[Map[closepointsvectors, data[[1 ;; 10000]]],2]; // AbsoluteTiming

(*{1.5870908, Null}*)

DeveloperPackedArrayQ@data

(*True*)

DeveloperPackedArrayQ@smalltest

(*False*)

ByteCount[smalltest]

(*40818664*)

packedsmalltest = DeveloperToPackedArray@smalltest;
DeveloperPackedArrayQ@packedsmalltest

(*True*)

ByteCount[packedsmalltest]

(*10204824*)


Is there a way to reduce memory usage, possibly keeping the array packed the whole time. Can I pad the array with zeroes as I go to prevent it from becoming ragged? I tried to take advantage of compile but it doesn't appear to work with Nearest.

-
I don't have M available at the moment. It might be helpful to use On["Packing"] to determine exactly where the unpacking is occurring as a start. –  Andy Ross Nov 30 '12 at 5:08
Did I miss the definition of nearestfunction? –  image_doctor Nov 30 '12 at 10:20
Is it the unpacking of the original array that causes memory to be exhausted or the generation of the set of closest vectors? If the latter you could always write out subsets of the closest vectors to disk as you work through the original data. –  image_doctor Nov 30 '12 at 10:31
@image_doctor, fixed sorry. –  s0rce Nov 30 '12 at 15:55

I suggest packing the output row-by-row as you create it:

data = RandomReal[10, {10000, 3}];
nf = Nearest[data];
f[x_] := DeveloperToPackedArray[(x - #) & /@ nf[x, {Infinity, 1}]]
smalltest = f /@ data;

ByteCount[smalltest]
(* 10611288 *)


As you can see, the memory usage is kept low - smalltest is an unpacked array of packed arrays.

However, you want the final output flattened, which will unpack the elements:

ByteCount[Flatten[smalltest]]
(* 32239408 *)


The workaround is to use Join to do the flattening:

smalltest = Join @@ Join @@ smalltest;

ByteCount[smalltest]
(* 9211384 *)

DeveloperPackedArrayQ[smalltest]
(* True *)

-