Compile with arbitrary array argument

Question

I need solve the task of classification of the points regarding their position. So I use such a function:

selector[mas_, reg_] := Select[mas, 
                      RegionMember[
                           Rectangle[{First@reg - st/2, Last@reg - st/2}, {First@reg + st/2, Last@reg + st/2}], 
                           #[[2]]] &];

It checks if second elements of mas entries are inside the rectangle defined by point reg. Finally it returns the list of right mas entries. There are ~4000 rectangles with ~250000 points and it evidently works too long.

I try to make the compiled version of such function:

  cs=Compile[{mas,_Real,4},{reg,_Real,2},
       Select[mas, 
             RegionMember[
                 Rectangle[{First@reg - st/2, Last@reg - st/2}, {First@reg + st/2, Last@reg + st/2}],
                 #[[2]]] &]

]

This function produce an error like "The first argument should be a 4 rank tensor". My array has a complicated form like mas={{1,{123.4,654.3},{12.1,21.2}},{...},...}. Of course it is not a good matrix but what is a problem with an arbitrary shape array? I was trying make the packed array from the initial one but it does not help.

The question is how to compile a function with arbitrary shape array of Real arguments?

Answer 1

You should not try to compile the brute-force algorithm. Instead you may use hierarchical data structures that allow the lookups in nearly linear time.

First read off all points pts= mas[[All,2]]. Then use Nearest to do the lookup (it employs a k-d-Tree). Simulating some data:

m = 4000;
n = 4000;
st = 0.5;
pts = RandomReal[{-1, 1}, {m, 2}];
rect = RandomReal[{-1 + st, 1 - st}, {n, 2}];
result = Nearest[pts -> "Index", rect, {\[Infinity], st/2}, DistanceFunction -> ChessboardDistance];

Now result[[k]] contains all the indices of the points in the k-th rectangle (and hence the indixec for the rows of your original array).

Here a test for correctness:

k = 13;

Graphics[{
  Point[pts],
  EdgeForm[Red], FaceForm[],
  Rectangle[rect[[k]] - st/2, rect[[k]] + st/2],
  Green, Point[pts[[result[[k]]]]]
  }]

For m = 250000 this takes two seconds on my machine.