Fastest way to go from linear index to grid index

Question

I'm sure this has been asked before but I'm interested in going from a position in a vector to the index in the grid version of the vector with given strides, for example, say I have the vector:

vec = {58, 94, 19, 68, 54, 77, 1, 18, 49, 20, 90, 44, 91, 89, 15, 0, 
   60, 18, 19, 44, 87, 5, 8, 42, 51, 55, 87, 71, 83, 68, 53, 58, 27, 
   17, 8, 14, 33, 58, 86, 3, 91, 66, 3, 16, 98, 84, 72, 98, 9, 30, 90,
    99, 15, 0, 82, 76, 86, 58, 77, 58};

And say I have strides {5, 4, 3}, the position 35 in the vector would correspond to the index {3, 4, 2}:

vec[[35]]

4

ArrayReshape[vec, {5, 4, 3}][[3, 4, 2]]

4

How can I get this index fast and in a vectorized fashion because I will have potentially many positions to extract?

Answer 1

Not intended to be competitive but fun for me to write.

unrank[d : {__Integer}][n_Integer] :=
   ⌊ 1 + d*Mod[(n - 1)/Reverse@FoldList[Times, Reverse@d], 1] ⌋

unrank[{5, 4, 3}][35]

SeedRandom[0]
r = RandomInteger[{2, 99}, 20];

unrank[r][1*^30]

{3, 4, 2}

{2, 8, 6, 6, 40, 38, 5, 51, 16, 12, 15, 34, 8, 45, 5, 28, 31, 12, 9, 8}

---

This time aiming for better performance specifically for application to lists of indexes.

unrankList[dim_List][n_List] :=
  1 + FoldList[QuotientRemainder[#[[1]], #2]\[Transpose] &, {n - 1}, 
     Reverse@dim][[-1 ;; 2 ;; -1, 2]]\[Transpose]

SeedRandom[0]
r = RandomInteger[{2, 99}, 20];
x = RandomInteger[{1, 1*^30}, 100];

unrankList[r][x]; // RepeatedTiming

{0.00119, Null}

---

Here is a derivative of your own method that seems to be a bit faster on my machine. Like your code it uses machine precision so it will become incorrect with very large indexes.

unrank3[n_Integer, d_] := unrank3[{n}, d]

unrank3[n_List, dim : {__Integer}] :=
  With[{tl = N@Reverse@FoldList[Divide, 1, Reverse@Rest@dim]},
    1 + Mod[⌊Partition[tl, 1].{n - 1}⌋, dim]\[Transpose]
  ]

Answer 2

I think this is what you want:

IntegerDigits[35 - 1, MixedRadix[{5, 4, 3}], 3] + 1

In general:

gridIndex[n_Integer, shape_List] := 
 IntegerDigits[n - 1, MixedRadix[shape], Length@shape] + 1

Answer 3

If you have enough memory, then a lookup table may be fastest:

shape = {5, 4, 3};
indices = Tuples[Range /@ shape];

Lookup is fast:

indices[[35]] // RepeatedTiming
(* {3.*10^-7, {3, 4, 2}} *)

Also, it seems that doing lots of lookups simultaneously is even faster (per lookup):

indices[[{22, 45, 35, 49, 36, 9, 9, 39, 59, 14}]] // RepeatedTiming

Answer 4

Here is the obligatory compiled version; it is not as fast as Roman's lookup method but it is also less memory hungry (in particular if indexing is supposed to be done into SparseArray whose dense version does not fit into memory).

A compiled helper function:

cf = Compile[{{n, _Integer}, {mods, _Integer, 1}},
   Block[{r = n - 1, d, m},
    Table[
     m = Compile`GetElement[mods, i];
     d = Quotient[r, m];
     r = r - d m;
     d + 1,
     {i, 1, Length[mods]}]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];
getInds[idx_, shape_] :=
 cf[idx, Reverse[Most[FoldList[Times, 1, Reverse[shape]]]]]

Usage example and timing test, comparing against swish's and Roman's proposals:

RandomSeed[123];
shape = RandomInteger[{1, 10}, 4];
idx = RandomInteger[{1, Times @@ shape}, 10000];

a = gridIndex[#, shape] & /@ idx; // RepeatedTiming // First

b = getInds[idx, shape]; // RepeatedTiming // First

indices = Tuples[Range /@ shape];
c = indices[[idx]]; // RepeatedTiming // First
a == b == c

0.429

0.00059

0.0000700

True

Answer 5

Here's what I came up with:

getSubindex[index_, stride_] := {
  Mod[index, stride, 1],
  Ceiling[index/stride]
  }
getIndex[index_, strides_] := 
 Reverse@FoldPairList[getSubindex, index, Reverse@strides]

This is comparable to swish's solution speed-wise:

gridIndex[1000, {3, 5, 4, 6}] // RepeatedTiming

{0.000061, {3, 2, 3, 4}}

getIndex[1000, {3, 5, 4, 6}] // RepeatedTiming

{0.000052, {3, 2, 3, 4}}

Answer 6

So after I posted the question last night I came up with a solution that is fast and vectorized: (note that if you're working with huge numbers you'll need to remove the N for accuracy, but you'll incur a huge speed penalty)

gifs[inds_, strides : {__Integer}] :=
 Module[
  {
   accstr,
   stride = strides,
   ind = inds - 1,
   moddable,
   modres
   },
  accstr =
   N@
    Append[
     Reverse@FoldList[Times, strides[[-1 ;; 2 ;; -1]]],
     1
     ];
  moddable = If[ListQ@inds, Map[ind/# &, accstr], ind/accstr];
  modres = 1 + Mod[Floor[moddable], stride];
  If[ListQ@inds, Transpose, Identity]@modres
  ]

Obligatory performance comparison:

tests = RandomInteger[{1, 60}, 100];

res = gifs[tests, {5, 4, 3}]; // RepeatedTiming // First

0.000053

res == gridIndex[tests, {5, 4, 3}]

True

gridIndex[tests, {5, 4, 3}]; // RepeatedTiming // First

0.0064

getIndex[tests, {5, 4, 3}]; // RepeatedTiming // First

0.00032

unrankList[{5, 4, 3}][tests]; // RepeatedTiming // First

0.00039

getInds[tests, {5, 4, 3}]; // RepeatedTiming // First

0.000063

Clearly vectorization is doing what it should and getting us the performance we'd expect (which is interestingly better than a compiled implementation on my machine)

Here's a more detailed performance analysis which shows we're long-term a little bit better than Mathematica's auto-parallelization in compiled functions:

RandomSeed[123];
shape = RandomInteger[{1, 10}, 4];
sizes = {1, 5, 10, 50, 100, 500, 1000, 5000, 10000, 50000, 75000, 
   100000, 125000, 150000, 500000};
idxs = RandomInteger[{1, Times @@ shape}, #] & /@ sizes;

testC =
  MapThread[
   {#, getInds[#2, shape]; // RepeatedTiming // First} &,
   {
    sizes,
    idxs
    }
   ];

testU =
  MapThread[
   {#, gifs[#2, shape]; // RepeatedTiming // First} &,
   {
    sizes,
    idxs
    }
   ];

ListLinePlot[
 {
  testC,
  testU
  },
 PlotLegends -> {"getInds", "gifs"},
 PlotRange -> All
 ]