11
$\begingroup$

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?

$\endgroup$
6
$\begingroup$

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]
  ]
$\endgroup$
  • $\begingroup$ Nice refactoring. I’m gonna leave this open maybe another day but then I’ll run all the tests and accept this if no one comes up with something faster. This answer certainly covers the most bases. $\endgroup$ – b3m2a1 Jan 4 at 22:07
  • $\begingroup$ @b3m2a1 Thank you. I posted a bit prematurely and had to make a couple of code edits after the fact. Please let me know if you find any failings with the current version. $\endgroup$ – Mr.Wizard Jan 4 at 22:11
13
$\begingroup$

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
$\endgroup$
  • 1
    $\begingroup$ @C.E. You're right, it just needs dimension length specification $\endgroup$ – swish Jan 4 at 6:58
  • 1
    $\begingroup$ This is a very nice solution, +1 $\endgroup$ – C. E. Jan 4 at 6:59
  • $\begingroup$ Do you know what version is required to run this code? In v10.1 I get {1 + IntegerDigits[34, MixedRadix[5], 3], 1 + IntegerDigits[34, MixedRadix[4], 3], 1 + IntegerDigits[34, MixedRadix[3], 3]} $\endgroup$ – Mr.Wizard Jan 4 at 15:52
  • 1
    $\begingroup$ This is clean but surprisingly incredibly slow... $\endgroup$ – b3m2a1 Jan 4 at 16:56
11
$\begingroup$

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
$\endgroup$
  • $\begingroup$ I don't and would need to construct it on the fly each time, but it is a clever simple work around. $\endgroup$ – b3m2a1 Jan 4 at 16:56
9
$\begingroup$

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

$\endgroup$
  • $\begingroup$ What do you mean? It leads to the same results as the other methods... $\endgroup$ – Henrik Schumacher Jan 4 at 17:27
  • $\begingroup$ Erm. I am getting 3D indices for shape = {5, 4, 3}... $\endgroup$ – Henrik Schumacher Jan 4 at 17:31
  • 1
    $\begingroup$ Ah, that's because the second argument of cf should be Reverse[Most[FoldList[Times, 1, Reverse[shape]]]] (which has to be computed only once; that's why I did not include it into cf). $\endgroup$ – Henrik Schumacher Jan 4 at 17:33
  • $\begingroup$ Ah I see. I'll add a tiny function building on cf to make this clear since I clearly did not read enough. $\endgroup$ – b3m2a1 Jan 4 at 17:36
  • 1
    $\begingroup$ Yeah, I was a bit surprised that a memory bound lookup table should be faster than a few integer operations. But as long as the lookup table fits into relatively low level cache, memory boundedness does not seem to be a problem. Why it switches back in the end is a complete miracle to me though (I did not experience it so far). $\endgroup$ – Henrik Schumacher Jan 4 at 18:05
7
$\begingroup$

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

$\endgroup$
5
$\begingroup$

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
 ]

enter image description here

$\endgroup$
  • $\begingroup$ This doesn't agree with the output of my function, e.g. I get {2, 8, 6, 6, 40, 38, 5, 51, 16, 12, 16, 36, 6, 33, 2, 19, 1, 13, 5, 2} from gifs for the example in my answer. $\endgroup$ – Mr.Wizard Jan 4 at 17:56
  • $\begingroup$ @Mr.Wizard hmm...I don't know what's happening there. It seems to only happen for very large numbers for some reason, though. All my other test cases against other answers pan out just fine. $\endgroup$ – b3m2a1 Jan 4 at 18:03
  • 1
    $\begingroup$ @Mr.Wizard Yeah it's my use of N to speed things up that introduces some numerical instability for huge numbers. $\endgroup$ – b3m2a1 Jan 4 at 18:04
  • 1
    $\begingroup$ @Mr.Wizard it's pretty good on my machine. Comparable to C.E.'s solution with foling Mod and Ceiling. $\endgroup$ – b3m2a1 Jan 4 at 18:27
  • 1
    $\begingroup$ A refactoring of your own method, unrank3, seems to be a little bit faster or my machine. Please test it. edit I just updated it to eliminate a Transpose operation. $\endgroup$ – Mr.Wizard Jan 4 at 21:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.