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

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.