Faster position list construction from a "take instructions" list

Question

Take some list, e.g., {1,2,3,4,5}, and a "take instructions" list, say {1,1,2,1,1}.

The latter is read as "take the first element from the list, drop it from the list, take the first element of the new list, drop it from the list..." while retaining the taken elements in the order taken., resulting in taking from the list the positions {1, 2, 4, 3, 5} in order.

I'm using

Module[{p = Range@Length@#},
   Reap[Scan[(Sow[p[[#]]]; p = Drop[p, {#}]) &, #]]] &@{1, 1, 2, 1, 1}

to do the conversion from the take instructions to the positions, but I have the sneaking suspicion there's a more efficient method that eludes me for the moment. The above returns the list of positions in the target, so target[[{list of positions}]] would give me the ultimate desired result. Something that avoids the intermediate step of getting positions would be even better (of course, one could just put the target in place of the range, but In a few tests I did, unpacked targets seemed to be slower than getting positions and using those on the target)...

Any ideas?

Edit: For timing, simply pad the above example "take instructions" with ones on the left, say 20K of them, since any solution w/b quick on examples as small as above. Also, the "take instructions" lists will always be well-formed (that is, never something like {1,2,1,1,2}) so no need to do any validation.

Update: Thanks to all for the very illuminating answers, it was tough picking the accept, but since "faster" was the criteria, 2012rcampion took the cake for overall average best performance for my needs. I really enjoyed Daniel's and Halirutan's detailed answers also.

I still have the lurking feeling there's a sneaky and super-fast method and will continue to ponder it.

Answer 1

My answer is based on a modification of a binary heap. Basically the construction looks something like this. We start with a binary tree:

Notice that if we label the nodes breadth-first, the labels have an interesting property. Each parent node $n$ has two children, $2n$ and $2n+1$. This also works in reverse: the parent of node $n$ is node $\left\lfloor n/2 \right\rfloor$. This means that we can store the nodes in a straight array, without having to store pointers.

We can renumber the nodes to correspond to their left-to-right position. Think of each node as a list, with the left subtree representing the left sublist, the right subtree the right sublist, and the node itself an element between the two lists. We end up with a numbering scheme like this:

This allows for fast binary-search lookup of each element (order $O(\log n)$). However, each removal would involve rebasing the heap. With this data structure this is an $O(n)$ operation, so we need an alternative.

What if, at each node, we store the number of elements $l$ in the left subtree, plus one? Then when looking at the $i$-th element in the list, we have three possibilites:

• $i<l$: element $i$ is in the left subtree, recurse left
• $i=l$: element $i$ is the current node, stop
• $i>l$: element $i$ is in the right subtree, recurse right

When I remove an element from the heap, I simply mark it as deleted, and then subtract one from all of the left-subtree-counts that contain it. Since this operation only affects the parent nodes, the operation is $O(\log n)$ (and our overall algorithm is $O(n\log n)$).

All that said, here is my algorithm:

robert = Compile[{{derangement, _Integer, 1}},
  With[{
    n = Length[derangement]
    },
   Module[{
     position = ConstantArray[1, n],
     totalCount = ConstantArray[1, n],
     leftCount = ConstantArray[1, n],
     result = ConstantArray[0, n]
     },
    Do[
     totalCount[[Quotient[i, 2]]] += totalCount[[i]];
     If[Mod[i, 2] == 0,
      leftCount[[Quotient[i, 2]]] = 1 + totalCount[[i]];
      ];,
     {i, n, 2, -1}];
    position[[1]] = leftCount[[1]];
    Do[
     position[[i]] = 
      position[[Quotient[i, 2]]] + 
       If[EvenQ[i], -totalCount[[i]] + leftCount[[i]], 
        1 + leftCount[[i]]] - 1,
     {i, 2, n}];
    Do[Module[{
        p = 1,
        i = derangement[[j]]
        },
       While[i != leftCount[[p]] || position[[p]] == 0,
        If[i > leftCount[[p]],
          i -= leftCount[[p]];
          p = 2 p + 1;,
          p = 2 p;
          ];
        ];
       result[[j]] = position[[p]];
       position[[p]] = 0;
       leftCount[[p]] -= 1;
       While[p > 1,
        If[Mod[p, 2] == 0,
          p = Quotient[p, 2];
          leftCount[[p]] -= 1;,
          p = Quotient[p, 2];
          ];
        ];
       ];
     , {j, n}];
    result
    ]
   ]
  ]

Here are my benchmarking results using AbsoluteTiming. I used this function to generate the take instructions/derangement lists:

randomDerange[n_] := Reverse@Table[RandomInteger[{1, i}], {i, n}]

This should give all possible derangements with equal probability. I may run mostly-ordered and -reverse-ordered cases later.

All the functions give the same results:

x = randomDerange[15];
Through[functions[x]] // Column

{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}
{6,9,2,8,7,13,3,5,15,4,10,1,11,12,14}

Although Daniel's solution is quite slow, it is clearly $O(n \log n)$ and crosses the others somewhere around a million elements.

Interestingly enough, all but two of the other solutions follow almost exactly the same $O(n^2)$ curve. I guess I shouldn't be too surprised, since the approaches are all fairly similar and Mathematica's builtins are well-optimized.

As for the two slow algorithms, halirutan's is slow because he's manually shifting around the array elements, and the compiler isn't optimizing them as well. Mike's is slow because he's using FoldList. (Changing it to Fold moves it down with the rest of the functional solutions.)

Incidentally, that makes Mike's $O(n^2)$ (!) in space complexity. All the others have very similar $O(n)$ memory usage, except for Daniel's: due to the way he stores the tree, his memory usage is around 5x higher (but still linear).

Answer 2

Preface

Below, you will find two different solutions. For understanding the problem itself, the first, iterative solution is better suited since it gives insight in how the solution can be found without directly executing the instructions given as input.

Iterative Solution

Detailed explanation

To explain the idea behind this approach let us work with a concrete example instruction

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

When simply executing these instructions on an input, we have to constantly drop elements from anywhere, even inside the list. My intent was to find an approach then does not do this but instead calculates the solution without doing any drop operations.

The important question on the way to the solution is

Can we (in any stage of the execution process) tell at which position the currently first element will be in the output?

The answer is yes. If you think a minute about it, you will see that the current first element of our list (no matter how many elements we have already dropped) will be used at the next 1 in the instruction set. Let's say we have already worked the first 1's in our example then our instruction set looks like this

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

and our input is probably

{., ., c, d, e, f, g}

Now it becomes clear, as long as we don't see a 1, we will not touch the c. Let's say we want to do this step now instead of using the 2 and dropping the d.

We can do that, but afterwards the leading 2 in the instruction set is no longer a 2 because we have altered the list. As it turns out after some sleepless minutes, this can be fixed by decreasing all indices that are left and including the 1 by exactly one. So what we get is the following altered instructions and input

{., ., 1, 2, 0, 2, 1}
{., ., ., d, e, f, g}

The important part is that we need to put the extracted c at the correct position in the result list which is the 5th position. This can easily be seen: When you work by dropping from the former input

{., ., 2, 3, 1, 2, 1}
{., ., c, d, e, f, g}

we drop one element at the time. Therefore, 2 would lead to inserting the d at the 3rd position, 3 would insert f at the 4th position in the result and finally, the following 1 would insert the (still leading) c at position 5.

Here is a more formal description of the algorithm. It calculates the final positions that can be directly used with Part on the input list:

1. in are the instructions, out is an empty list with the same length as in, c = 1 is the current step
2. Find pos, the first 1 in in beginning from left. Return out if this is not possible anymore.
3. Set out[[pos]]=c
4. Decrease all positions of in that are small-equal pos by 1
5. Goto (2.)

Here is a complete run for the sample instructions

Improvements

As can easily seen when working manually through an example, all left parts of in that solely contain numbers smaller 1 don't need to be processed any more. In the blackboard image I have inserted a small blue line at those positions. The implementation below contains this improvement (look out for start) but for the explanation I have left it out since it would have probably been confusing.

Implementation

fc = Compile[{{instr, _Integer, 1}},
  Module[
   {in = instr, out = instr, l = Length[instr], 
    start = 1, j = 1, i = 1, p = 1},
   While[i <= l,
    While[j <= l,
     If[in[[j]] == 1,
      out[[j]] = p++;
      in[[start ;; j]]--;
      Break[];
      ];
     j++;
     ];
    If[j > l, Break[]];
    While[in[[start]] < 1 && start < l, start++];
    j = start;
    i = start;
    ];
   out
   ]
];

This method will be superior on many inputs, but slow when there are many long ways in the input. This means, a worst case scenario is the input {5,4,3,2,1}. Let us compare 3 data-sets with rashers original solution:

worstCase = Reverse[Range[10^4]];
First[AbsoluteTiming[rasher[worstCase]]]
First[AbsoluteTiming[fc[worstCase]]]
(* rasher 0.04, halirutan 2.61 *)

averageCase = Flatten[ConstantArray[{5, 4, 3, 2, 1}, 10^4]];
First[AbsoluteTiming[rasher[averageCase]]]
First[AbsoluteTiming[fc[averageCase]]]
(* rasher 0.57, halirutan 0.04 *)

bestCase = ConstantArray[1, 10^5];
First[AbsoluteTiming[rasher[bestCase]]]
First[AbsoluteTiming[fc[bestCase]]]
(* rasher 2.51, halirutan 0.06 *)

Functional Style Solution

There is another variant, that is more functional because it passes the shrinking arrays around instead of assigning them. The advantage is probably more the code-style, although it performs as fast as rashers solution.

This function is calculating the output instead of the final positions, for the obvious reason that the calculation of the positions is as hard as directly calculating the output:

take[l_, p_] := 
  Flatten[Reverse[Reap@Fold[(Sow[Part[##]]; Delete[##]) &, l, p]]]

take[{1, 2, 3, 4, 5}, {1, 1, 2, 1, 1}]
(* {1, 2, 4, 3, 5} *)

Answer 3

I have a tree-based method that has the right asymptotics but a very high coefficient. The upshot being, it will not compete with other methods until we get past 10^6 or so in list size. With considerable work that tree structure could be flattened so that Compile might be brought into play.

The basic tree layout is {left subtree, node, right subtree} where a node is a 4-tuple of the form {value, size of left subtree, size of right subtree, present}. That last is a boolean flag indicating whether the node is still around or has been removed.

{leftsubtree, node, rightsubtree} = Range[3];
emptyTree = {};

leftChild[tree_] := tree[[leftsubtree]]
SetAttributes[setLeftChild, HoldFirst]
setLeftChild[tree_, left_] := tree[[leftsubtree]] = left

rightChild[tree_] := tree[[rightsubtree]]
SetAttributes[setRightChild, HoldFirst]
setRightChild[tree_, right_] := tree[[rightsubtree]] = right

value[tree_] := tree[[node, 1]]

numLeft[emptyTree] := 0
numLeft[tree_] := tree[[node, 2]]
SetAttributes[decNumLeft, HoldFirst]
decNumLeft[node_] := node[[2]]--

numRight[emptyTree] := 0
numRight[tree_] := tree[[node, 3]]
SetAttributes[decNumRight, HoldFirst]
decNumRight[node_] := node[[3]]--

isPresent[tree_] := tree[[node, 4]]
SetAttributes[setIsPresent, HoldFirst]
setIsPresent[node_, flag_] := node[[4]] = flag

makeNode[val_] := {val, 1, 0, True};
makeNode[val_, numl_, numr_] := {val, numl, numr, True}

makeTree[val_] := {emptyTree, makeNode[val], emptyTree}

makeTree[val_, left_, right_] := {left, 
  makeNode[val, numLeft[left] + numRight[left] + 1, 
   numLeft[right] + numRight[right]], right}

newTree[tlist_List] := Block[
  {len, len2, left, right},
  len = Length[tlist];
  If[len == 0, Return[emptyTree]];
  If[len == 1, Return[makeTree[tlist[[1]]]]];
  len2 = Quotient[len, 2];
  left = newTree[Take[tlist, len2]];
  right = newTree[Drop[tlist, len2 + 1]];
  makeTree[tlist[[len2 + 1]], left, right]
  ]

SetAttributes[removeElement, HoldFirst];
removeElement[tree_, elem_Integer] := Module[{newnode, val},
  newnode = tree[[node]];
  If[isPresent[tree] && numLeft[tree] == elem, val = value[tree];
   setIsPresent[newnode, False];
   decNumLeft[newnode];
   tree[[node]] = newnode;
   Return[val];];
  If[numLeft[tree] < elem, decNumRight[newnode];
   tree[[node]] = newnode;
   removeElement[tree[[rightsubtree]], elem - numLeft[tree]]
   ,(*else*)
   decNumLeft[newnode];
   tree[[node]] = newnode;
   removeElement[tree[[leftsubtree]], elem]]]

Here is an example. For the "take" list we have a max jump size of 100. To be safe we set the final 101 elements equal to 1.

n = 10^4;
maxJump = 100;
originals = Range[n];
takelist = 
  Join[RandomInteger[{1, maxJump}, n - maxJump - 1], 
   ConstantArray[1, maxJump + 1]];

AbsoluteTiming[tree = newTree[originals];]
AbsoluteTiming[
 result = Table[removeElement[tree, takelist[[j]]], {j, n}];]

(* Out[66]= {0.110308, Null}

Out[67]= {2.858037, Null} *)

This is, suffice it to say, hugely slower than the other solutions presented thus far. But the complexity is O(n log n) and it does in fact cross over with the originally presented method somewhere between one and two million.

Edit 1

Since we do not need an on line (that is, dynamically changing) structure, it is faster to use a "flat" tree. We have basically the same tree as I used above but flattened into a tensor form.

Here is a reference implementation. I simplified the code somewhat by removing things I wasn't using in the first place.

makeNode[val_] := {val, 1,(*0,*)1};
makeNode[val_, numl_] := {val, numl, 1}

makeTree[val_] := {makeNode[val]}

makeTree[val_, left_, right_] := Module[{node},
  node = makeNode[val, Length[left] + 1];
  Join[left, {node}, right]]

newTree[tlist_List] := 
 Block[{len, len2, left, right}, len = Length[tlist];
  If[len == 0, Return[{}]];
  If[len == 1, Return[makeTree[tlist[[1]]]]];
  len2 = Quotient[len, 2];
  left = newTree[Take[tlist, len2]];
  right = newTree[Drop[tlist, len2 + 1]];
  makeTree[tlist[[len2 + 1]], left, right]]

SetAttributes[removeElement, HoldFirst];
removeElement[tree_, elem_Integer] := Catch[Module[
   {node, mid = Quotient[Length[tree], 2] + 1, start = 0, 
    end = Length[tree] + 1, j = 0, sum = 0},
   While[True && j < 20, j++;
    If[mid == start || mid == end, Throw[-1]];
    node = tree[[mid]];
    If[node[[3]] == 1 && sum + node[[2]] == elem,
     tree[[mid, 3]] = 0;
     tree[[mid, 2]]--;
     Throw[node[[1]]];];
    If[sum + node[[2]] < elem,
     sum += node[[2]];
     {start, mid} = {mid, mid + Quotient[(end - 1 - mid), 2] + 1};
     ,(*else*)
     tree[[mid, 2]]--;
     {end, mid} = {mid, start + Quotient[(mid - start - 1), 2] + 1};];
    ]
   ]]

To be continued.

Edit 2

We can now put this through Compile as below.

newTreeC = Compile[{{tlist, _Integer, 1}},
   Catch[Module[{len, len2, left, right},
     len = Length[tlist];
     If[len == 1, Throw[{{tlist[[1]], 1, 1}}]];
     If[len == 2, Throw[{{tlist[[1]], 1, 1}, {tlist[[2]], 2, 1}}]];
     len2 = Quotient[len, 2];
     left = newTreeC[tlist[[1 ;; len2]]];
     right = newTreeC[tlist[[len2 + 2 ;; -1]]];
     Join[left, {{tlist[[len2 + 1]], Length[left] + 1, 1}}, 
      right]]], {{newTreeC[_], _Integer, 2}}, 
   RuntimeOptions -> "Speed"];

takeByListC = Compile[{{otree, _Integer, 2}, {tlist, _Integer, 1}},
   Module[
    {tree = otree, node, mid, start, end, sum, elem},
    Table[elem = tlist[[k]];
     mid = Quotient[Length[otree], 2] + 1;
     start = 0;
     end = Length[otree] + 1;
     sum = 0;
     Catch[While[True,
       If[mid == start || mid == end, Throw[-1]];
       node = tree[[mid]];
       If[node[[3]] == 1 && sum + node[[2]] == elem,
        tree[[mid, 3]] = 0;
        tree[[mid, 2]]--;
        Throw[node[[1]]],
        If[sum + node[[2]] < elem,
          sum += node[[2]];
          {start, mid} = {mid, mid + Quotient[(end - 1 - mid), 2] + 1};
          ,(*else*)
          tree[[mid, 2]]--;
          {end, mid} = {mid, 
            start + Quotient[(mid - start - 1), 2] + 1}];
        ]]]
     , {k, Length[tlist]}]], RuntimeOptions -> "Speed"];

This test, on a 6 year old laptop, will at least give an idea of performance. I think I may finally have all the bugs out but caveat emptor.

n = 10^6;
maxJump = 100;
originals = Range[n];
takelist = 
  Join[RandomInteger[{1, maxJump}, n - maxJump - 1], 
   ConstantArray[1, maxJump + 1]];

AbsoluteTiming[tree = newTreeC[originals];]
AbsoluteTiming[result = takeByListC[tree, takelist];]

(* Out[45]= {5.290000, Null}

Out[46]= {18.920000, Null} *)

Might be faster with CompilationTarget -> "C". If I had a machine that knew what a C compiler looked like.

Edit 3

[I guess I'm spending too much time on this. Nice problem it has turned out to be.]

Here is an improved compiled version. Still around twice as slow as the code by @2012rcampion.

newTreeC = Compile[{{len, _Integer}}, Catch[Block[
     {len2, left, right},
     If[len == 1, Throw[{{1, 1}}]];
     If[len == 2, Throw[{{1, 1}, {2, 1}}]];
     len2 = Quotient[len, 2];
     left = newTreeC[len2];
     right = newTreeC[len - len2 - 1];
     Join[left, {{Length[left] + 1, 1}}, right]]],
   {{newTreeC[_], _Integer, 2}},
   RuntimeOptions -> "Speed", CompilationTarget -> "C"];

takeByListC = 
  Compile[{{otree, _Integer, 2}, {tlist, _Integer, 1}}, 
   Block[{tree = otree, node, mid, start, end, sum, elem, 
     saveend = Length[otree] + 1, 
     savemid = Floor[Length[otree]/2.] + 1},
    Table[
     elem = tlist[[k]];
     mid = savemid;
     start = 0;
     end = saveend;
     sum = 0;
     Catch[While[True,
       node = tree[[mid]];
       If[node[[2]] == 1 && sum + node[[1]] == elem,
        tree[[mid, 2]] = 0;
        tree[[mid, 1]]--;
        Throw[mid]
        ];
       If[sum + node[[1]] < elem,
        sum += node[[1]];
        start = mid;
        mid = mid + Floor[(end - 1 - mid)/2.] + 1;
        ,(*else*)
        tree[[mid, 1]]--;
        end = mid;
        mid = start + Floor[(mid - start - 1)/2.] + 1;
        ];
       ]]
     , {k, Length[tlist]}]],
   RuntimeOptions -> "Speed", CompilationTarget -> "C"];

Here is code similar to the method of @halirutan, but moving result values instead of index list values. It is slightly longer than need be so as to avoid using RotateRight. That usage would make is somewhat slower; I guess it involves calls to the run time library.

iterativeTakeByList = Compile[{{tlist, _Integer, 1}}, Block[
    {len = Length[tlist], result, elem},
    result = Range[len];
    Do[
     elem = result[[j + tlist[[j]] - 1]];
     Do[result[[j + k + 1]] = result[[j + k]], 
       {k, tlist[[j]] - 2, 0, -1}];
      result[[j]] = elem;
     , {j, len}
     ];
    result
    ],
RuntimeOptions -> "Speed", CompilationTarget -> "C"];

Here is a test on a list of a million elements. With a maximum jump of 1000 it is still competitive with the best of the tree-based codes.

n = 10^6;
maxJump = 1000;
takelist = 
  Join[RandomInteger[{1, maxJump}, n - maxJump - 1], 
   ConstantArray[1, maxJump + 1]];
AbsoluteTiming[result = iterativeTakeByList[takelist];]

(* Out[179]= {0.762602, Null} *)

Answer 4

Module[{p = Range[Length@#]}, 
   Reap@Fold[(Sow[#[[#2]]]; Drop[#, {#2}]) &, p, #]] &@{1, 1, 2, 1, 1}
(* {{}, {{1, 2, 4, 3, 5}}} *)

Answer 5

If I've understood correctly you appear to be wanting a nested delete in which you keep the deleted items in their preserved order of deletion. This solution is probably not any easier than yours but here goes:

Rest@FoldList[{Join[#1[[1]], {#1[[2, #2]]}], 
    Delete[#1[[2]], #2]} &, {{}, Range[5]}, {1, 1, 2, 1, 1}]

(*  {{{1}, {2, 3, 4, 5}}, {{1, 2}, {3, 4, 5}}, {{1, 2, 4}, {3, 5}}, {{1, 
   2, 4, 3}, {5}}, {{1, 2, 4, 3, 5}, {}}}  *)

In the above the you start with two lists {{},Range[5]} and your "take instructions" list. The first empty list gets populated with the deleted elements of the second list -- based on the take instructions list.

Answer 6

l = {1, 2, 3, 4, 5};
t = {1, 1, 2, 1, 1};
l2 = t;
Fold[(Delete[l2[[#2]] = #1[[t[[#2]]]]; #1, t[[#2]]]) &, l, Range[5]];
l2
(*{1, 2, 4, 3, 5}*)

Answer 7

I failed to think of anything with competitive performance, however in simply working through the problem I produced a function that is somewhat different from the rest. Despite being slower than the original perhaps it will inspire something better.

f1 = 
  Module[{a = Range @ Length @ #},
    Scan[(a[[# ;;]] = RotateLeft @ a[[# ;;]]) &, #];
    a
  ] &;

Answer 8

Using FoldPairList:

list = {1, 2, 3, 4, 5};
instr = {1, 1, 2, 1, 1};
FoldPairList[{#1[[#2]], Delete[#1, #2]} &, list, instr]

{1, 2, 4, 3, 5}