Why is there no PositionFunction in Mathematica?

Question

Too often I have seen the programs of inexperienced users greatly slowed by using Position in an iterative fashion, when far faster would have been to compute a look-up table for positions beforehand.

Mathematica provides such functionality for Nearest and Interpolation "out of the box" with the syntax of function[data] along with dedicated functions NearestFunction and InterpolatingFunction.

There is no equivalent with Position[data] producing a PositionFunction object. I observe that Position appears to be used more often than either of the other functions so this seems like a regrettable omission. Two complications I can think of are:

• the levelspec of Position

• the handling of patterns

Based on this I ask:

1. What other complications are there for creating a PositionFunction?

2. How best can such a function be implemented?

3. How can the utility and performance of the function be maximized?

4. All limitations considered would such a function be valuable?
   Why is there no such functionality in Mathematica?

Answer 1

I see no mention of the new-in-10 PositionIndex in the other answers, which takes a list (or association) of values and returns a 'reverse lookup' that maps from values in the list to the positions where they occur:

In[1]:= index = PositionIndex[{a, b, c, a, c, a}]

Out[1]= <|a -> {1, 4, 6}, b -> {2}, c -> {3, 5}|>

It doesn't take a level spec yet (though I do want to add that).

In any case, the returned association is already a function in some sense because of the function application way of doing key lookup in associations.

So in the above example you could write index[a] to get a list of places where a occurred.

Answer 2

This is perhaps a place to start:

position[expr_, level_: 1] :=
 With[{positionData =
    SortBy[
      #[[1, 1]] -> #[[All, 2]] & /@
       GatherBy[Extract[expr, #, Verbatim] -> # & /@ Position[expr, _, level], First],
      Min[Length /@ #[[2]]] &
    ] // Dispatch},
  Replace[#, positionData] &
 ]

The second argument controls the depth of the indexing. An example:

f = position[x^2 + y^2 + q_^r_, 3];

f[2]

{{1, 2}, {2, 2}}

Due to the use of Verbatim patterns are matched literally, which deviates from normal Position behavior:

f[_]

{{3, 1, 2}, {3, 2, 2}}

Answer 3

In the function I propose, I build an association with keys that are not supposed to be evaluated. There are some issues with this, see this answer by Taliesin, with the following quote.

generally this just sounds like a dangerous and confusing game to play, to me.

I think the function presented in this answer deals with the complications you mention reasonably well. It uses an option to set the levelspec. To see how patterns are handled, see the section Verbatim.

Concerning point (1),(2) and (3): Before there were a lot of additional complications. But I now that we have Association we no longer have to deal with those. Making this work with held expressions is just a matter of being thorough in surrounding expressions with Hold and Unevaluated. My intuition is also that Association should have better performance than a Dispatch table or something similar. An Association should be unbeatable in terms of how long it takes to look up a particular sub-expression. But maybe we should do a proper comparison.

ClearAll[positionIndexGeneral]
Options[positionIndexGeneral] = {Heads -> True};
SetAttributes[positionIndexGeneral, HoldAll];
positionIndexGeneral[expr_, lev_: {1,Infinity}, OptionsPattern[]] := 
 Module[{subExprs, positions, len, together, gathered, hGathered, 
   gatheredSubExprs, gatheredPos},
  subExprs = 
   Level[Unevaluated@expr, lev, Hold, Heads -> OptionValue[Heads]];
  positions = 
   Position[Unevaluated@expr, _, lev, Heads -> OptionValue[Heads]];
  len = subExprs // Length;
  together = Transpose[{List @@ Hold /@ subExprs, positions}];
  gathered = GatherBy[together, First];
  hGathered = Hold@Evaluate@gathered;
  gatheredSubExprs = hGathered[[All, All, 1, 1, 1]];
  gatheredPos = gathered[[All, All, 2]];
  AssociationThread @@ {Unevaluated @@ gatheredSubExprs, gatheredPos}]

Example:

a = 3;
positionAssoc = 
 positionIndexGeneral[{a, 2, {3, 4, a}}]
positionAssoc[Unevaluated[a]]

{{1},{3,3}}

corresponding to

Position[Unevaluated@{a, 2, {3, 4, a}}, Unevaluated[a]]

{{1},{3,3}}

Verbatim

Note that in general we are simulating how Position works with Verbatim.

positionAssoc = positionIndexGeneral[{a, 2, {3, 4, a_}}]
positionAssoc[a_]

{{3,3}}

Corresponding to

Position[Unevaluated@{a, 2, {3, 4, a_}}, Verbatim[a_]]

{{3,3}}

To simulate how Position works without Verbatim in this way is probably not very useful. There are infinitely many patterns against which an expression can be tested, so of course we cannot make a big lookup table. For a very specific pattern like List | Hold we might make some specialised code that looks up both List and Hold in the association.

Timing

My function can kind of compete with a specialised function by Mr.Wiz in the 1D case, and of course it dwarfs the built in PositionIndex for large data.

f[x_] := AssociationThread @@ {Hold[
       Unevaluated[x]][[1, {1}, #[[All, 1]]]], #} &@
  GatherBy[Range@Length@x, Hold[x][[{1}, #]] &]

Now let's make some data and compare

data = RandomInteger[999, 1*^5];
(jacobGen = 
    positionIndexGeneral[Evaluate@data, {1, 1}, Heads -> False]) // 
  Timing // First
(mma1D = PositionIndex[data]) // Timing // First
(wiz1D = f[data]) // Timing // First
Position[data, 115] === jacobGen[115] === List /@ wiz1D[115] === 
 List /@ mma1D[115]

0.214873
0.174309
0.164100
True

data = RandomInteger[10, 1*^5];
(jacobGen = 
    positionIndexGeneral[Evaluate@data, {1, 1}, Heads -> False]) // 
  Timing // First
(mma1D = PositionIndex[data]) // Timing // First
(wiz1D = f[data]) // Timing // First

0.235508
4.119624
0.153041

data = RandomInteger[10, 1*^6];
(jacobGen = positionIndexGeneral[data, {1, 1}, Heads -> False]) // 
  Timing // First
(wiz1D = f[data]) // Timing // First

2.294256
1.703060

Possible improvement

When we only want the expressions at level 1, the function provided by Mr.Wizard is faster. With some good metaprogramming it should be possible to get the best of both worlds.

Appendix

It would of course have been cooler to write something like

ClearAll[positionIndexGeneral]
Options[positionIndexGeneral] = {Heads -> True};
positionIndexGeneral[expr_, lev_: {1,Infinity}, OptionsPattern[]] :=
 AssociationThread @@
      {
       Unevaluated @@ #2[[All, All, 1, 1, 1]]
       ,
       Function[{x}, x[[#]] & /@ #[[All, All, 2]]]@
        Position[expr, _, lev, Heads -> OptionValue[Heads]]
       } &[#, Hold@Evaluate@#] &@
  GatherBy[Transpose[{List @@ Hold /@ #, Range[# // Length]}] &@
    Level[expr, lev, Hold, Heads -> OptionValue[Heads]], First]

but I prefer the style with Module(/Block when possible) for debugging, as well as to immediately see what happens first.

Answer 4

This is a similar approach to Mr. Wizard's I suppose, but using the function's DownValues instead of a dispatch table to store the rules.

One major difference is that this code won't work with data containing pattern objects (I guess this might be fixable with Verbatim).

The expression is traversed using MapIndexed, for each part visited the position is Sowed, with the part's value as the tag. The actual downvalues are set afterwards using the third argument of Reap (the first time I've ever used it I think).

makePositionFunction[f_Symbol, data_, level_: {-1}] := Block[{},
  ClearAll[f];
  Reap[
   MapIndexed[Sow[#2, #1] &, data, level, Heads -> True],
   _, (f[#] = #2) &];
  f[other_] := Position[data, other, level]]

Example:

data = RandomInteger[1000, {3000, 20}];

makePositionFunction[pos, data];

First @ Timing[test1 = Table[pos[i], {i, 1000}]]
(* 0. *)

First @ Timing[test2 = Table[Position[data, i], {i, 1000}]]
(* 3.532 *)

test1 == test2
(* True *)

The "PositionFunction" defaults to using plain old Position for any search patterns which have not been precomputed:

pos[n_ /; n < 20] // Length
(* 1225 *)

Answer 5

I think the complexities of things like PatternTest will obstruct any kind of data structure for searching for general pattern matches. I think some regularity, either with the pattern as in Mr Wizard's answer or with the expression as in Leonid Shifrin's remark, will be needed in order to beat using Position. Jacob Akkerboom has pointed out problems with what to do if variables in expression change after the position function is computed.

That said, here's a start along the lines of Mr. Wizard's answer, using literal matches (literally, as you will see ;-). When he mentioned NearestFunction in the question, I thought, I wonder if there is a way to measure the "pattern-distance" between expressions. I don't think it's likely in general, but you can in the case of exact matches. Recall that Nearest uses EditDistance when the data are strings.

ClearAll[position2, positionFunction]; 
positionFunction[nf_NearestFunction, expr_][x_] := 
 With[{arg = ToString[FullForm @ x]},
  With[{pos = nf[arg]}, 
   If[pos =!= {} && ToString[FullForm @ Extract[expr, Prepend[First @ pos, 1]]] === arg, 
    pos, {}]]];
Format[positionFunction[nf_NearestFunction, expr_]] := positionFunction["<>", Short[expr]];
position2[expr_, level_: 1] := positionFunction[
  Nearest[ToString[FullForm @ Extract[expr, #]] -> # & /@ Position[expr, _, level]],
  HoldForm[expr]]

Mr.Wizard's example:

pf = position2[x^2 + y^2 + q_^r_, 3]

positionFunction["<>", x^2 + y^2 + q_^r_]

pf[2]

{{1, 2}, {2, 2}}

pf[_]

{{3, 1, 2}, {3, 2, 2}}

Another example:

gf = position2[Plot3D[x^2 + y^2, {x, -2, 2}, {y, -2, 2}], Infinity]

gf[Polygon]

{{1, 2, 1, 1, 2, 1, 1, 0}, {1, 2, 1, 1, 2, 1, 2, 0},
 {1, 2, 1, 1, 2, 1, 3, 0}, {1, 2, 1, 1, 2, 1, 4, 0}}