# Why is there no PositionFunction in Mathematica?

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

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?

• I belong to the team of Position[]-users... :-(
– Rod
May 21, 2013 at 13:23
• I find it funny that I really only started to use Position after participating here. I think this site has warped my sensibilities. :D May 21, 2013 at 13:56
• @ruebenko That is not the case I'm thinking of. The usual misuse is something like: a = RandomInteger[{1, 1000}, {5000}]; Table[{i, Position[a, i]}, {i, 1, 1000}] // Timing // First when magnitudes faster would be Sort[{#[[1, 1]], #[[All, 2]]} & /@ GatherBy[Transpose[{a, Range@Length@a}], First]] // Timing // First May 22, 2013 at 7:56
• There does occur to me a 'reason' for why there is no PositionFunction, and that is that Nearest and Interpolation are "subject matter" functions whereas Position is more a language thing. E.g. hypothetically, Nearest could be constructed from a variety of data structures because all that matters is the mathematical distance between one unit and another, whereas Position is entirely about the specific list structure you are using. Not saying I'd agree with this, but there is this conceptual difference which may have something to do with the absence of a PositionFunction.
– amr
May 23, 2013 at 1:52
• @Mr.Wizard I have updated my answer and I feel it is now quite complete. I did look a Szabolcs GatherBy inversion trick, but I think we cannot efficiently use it in the general case. The only thing that is missing now is to use your idea when level 1 and Heads-> False is specified, but I'd say that is a minor point. Jul 16, 2014 at 12:02

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.

• This would be helpful (so +1), even without levelspec, except that it is horribly slow. Hopefully that problem can quickly be fixed for the next point release. Jul 15, 2014 at 0:16
• By the way, as an insider would you please consider answering this?: (44189) I am unable to. Jul 15, 2014 at 1:02
• @Mr.Wizard it is only horribly slow when there are a large fraction of repeated values. For small lists or large lists with unique values it is plenty fast. But yes, it will be fixed. Jul 15, 2014 at 2:22
• I have Accepted this answer, at least for now. There are certain aspects of the question that I hope to see natively addressed some day, e.g. levelspec. May 15, 2015 at 8:00

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}}

• Hm.. Are you aware that setting x=y breaks the positionFunction? i.e. f[y] -> {1,1} . Verbatim does not have HoldAll. Maybe use HoldPattern (as well as Verbatim?). I guess it would be even better if you could shield the symbols from updates. Compare, with your definitions, {Drop[Trace[f, TraceOriginal -> True], 2], Update[y]; Drop[Trace[f, TraceOriginal -> True], 2]} May 21, 2013 at 16:36
• @Jacob Good point; I should add that, or change to DownValues on a HoldAll function, etc. I only intended this to be a rough start. May 21, 2013 at 18:05
• Maybe I should have added that I respect that it is rough start. It seems like a good start though :) May 21, 2013 at 18:08
• Are you a tribal Wizard now :D? Or witch doctor, perhaps (I call dibs on that name!)
– rm -rf
May 21, 2013 at 23:16
• @rm-rf Don't miss the videos. May 22, 2013 at 5:54

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]
SetAttributes[positionIndexGeneral, HoldAll];
positionIndexGeneral[expr_, lev_: {1,Infinity}, OptionsPattern[]] :=
Module[{subExprs, positions, len, together, gathered, hGathered,
gatheredSubExprs, gatheredPos},
subExprs =
positions =
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]
positionIndexGeneral[expr_, lev_: {1,Infinity}, OptionsPattern[]] :=
{
Unevaluated @@ #2[[All, All, 1, 1, 1]]
,
Function[{x}, x[[#]] & /@ #[[All, All, 2]]]@
} &[#, Hold@Evaluate@#] &@
GatherBy[Transpose[{List @@ Hold /@ #, Range[# // Length]}] &@


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

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 *)

• +1 for the post. However, I considered MapIndexed and went with Position and Extract for one reason: held expressions (or subexpressions). To make this robust I think you'll need something additional, e.g. MapIndexed onto an explicitly held expression and then a Cases to pull out those parts. You'll also have to consider how to keep parts extracted from a held expression from evaluating. May 22, 2013 at 14:38
• @Mr.Wizard, good point, I didn't consider held expressions. I suppose one could, with sufficient code gymnastics, work with held expressions in MapIndexed - but then why bother when you've already shown how to do it neatly with Position and Extract. Great question, by the way. May 22, 2013 at 16:02
• @Mr.Wizard I find it very frustrating that there is no function like ScanIndexed. That would work perfectly with held expressions. I have spent a lot of time trying to find an alternative for held expressions to the construct with (Inner, my addition) Position and Extract (that I suppose will soon appear in my answer), to try do something with Scan and/or MapIndexed. Maybe we can discuss things soon? May 29, 2013 at 15:27

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}}

• I love how your positionFunction looks. Also, you code is so short compared to mine ;). May 23, 2013 at 21:37
• @JacobAkkerboom Thanks. I think to do this right, the code might have to be long, though. May 23, 2013 at 23:03