Skip to main content

performance tuning - Compiling Map over expression that yields a ragged array


I'm trying to speed up a function that looks in the neighborhood of each 3D point in a large dataset and finds all the points within 1 unit in each direction, x, y, z.


I've started by using Select to find the points around a specified point and then Map this function over the dataset.



First, lets make up some "data" and then define the findpoints function.


data = RandomReal[10, {10^4, 3}];

findpoints[data_, pt_]:= Block[{x = pt[[1]], y = pt[[2]], z = pt[[3]]},
Select[data,
x - 1 <= #[[1]] <= x + 1 && y - 1 <= #[[2]] <= y + 1 &&
z - 1 <= #[[3]] <= z + 1 &]
]

Now lets map it over the random data.



Map[findpoints[data, #] &, data]; // AbsoluteTiming

This is quite slow, I gave up waiting, so I started by compiling the findpoint function.


findpointsC = 
Compile[{{data, _Real, 2}, {pt, _Real, 1}},
Block[{x = pt[[1]], y = pt[[2]], z = pt[[3]]},
Select[data,
x - 1 <= #[[1]] <= x + 1 && y - 1 <= #[[2]] <= y + 1 &&
z - 1 <= #[[3]] <= z + 1 &]
], CompilationTarget -> "C", RuntimeAttributes -> {Listable},

Parallelization -> True, RuntimeOptions -> "Speed"];

output= Map[findpointsC[data, #] &, data]; // AbsoluteTiming

(*{3.7542147, Null}*)

Much faster, great.


Now I figured since compiling was so beneficial here and in other functions I found that by compiling the Map function I can gain some more speed, lets compile the whole thing.


mapfindpointsC = 
Compile[{{data, _Real, 2}}, Map[findpointsC[data, #] &, data],

CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed"];

mapfindpointsC[data]

(*CompiledFunction::cflist: Nontensor object generated; proceeding with uncompiled evaluation. >>*)

However, as per the error message, in this case the output is ragged.


TensorQ[output]


(*False*)

I've tried to flatten the output from map but nothing I can think of works. Is it possible to compile this expression.



Answer



Please compare


Map[findpointsC[data, #] &, data]; // AbsoluteTiming
(* {3.567939, Null} *)

and


findpointsC[data, data] ; // AbsoluteTiming

(* {0.841954, Null} *)

And of course is


Map[findpointsC[data, #] &, data] == findpointsC[data, data] 
(* True *)

The explanation is simple: You typed RuntimeAttributes->{Listable} and Parallelization but you have not used it.


Explanation for your error message


The problem with your Map arises because Compiled functions can only work with rectangular arrays. Since every point has (most likely) a different number of points in its neighborhood, the result array is ragged, or not rectangular.


It's not enough to use Flatten because even then, you still create such a ragged array. One (slower) way out of this is to collect all points in a linear list, but then you loose the information, which point-set was created by which point. This cannot be the goal, if the solution is as simple as described above. Anyway, as you might have seen we had a topic here about Internal`Bag which can be used to collect an unknown number of elements. Please see the discussion in the thread



Internal`Bag inside Compile


Based on this, you could make an all in one compiled function which I only give for reference and not for usage. It is much slower than the Listable approach. I use the already compiled function findpointsC by using With and giving the option


CompilationOptions -> {"InlineCompiledFunctions" -> True}

I have replaced your Map with a Do loop which runs over all points and finds the nearest neighbors. To insert such a list of points into an Internal`Bag we have to Flatten it which can be done with the 3rd parameter of Internal`StuffBag. To reshape the flat list at the end into a list of points I use Partition


mapfindpointsC =
With[{fpC = findpointsC},
Compile[{{data, _Real, 2}},
Module[{result = Internal`Bag[Most[{.0}]]},
Do[Internal`StuffBag[result, fpC[data, pt], 2], {pt, data}];

Partition[Internal`BagPart[result, All], 3]
], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True, RuntimeOptions -> "Speed",
CompilationOptions -> {"InlineCompiledFunctions" -> True}]
]

mapfindpointsC[data]; // AbsoluteTiming

(* {3.831499, Null} *)

Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...