0
$\begingroup$

I wish to use a preassigned variable to index a large array many times. For example:

a = {i,j+1,k-1};
b = {i+1,j,k-2};

I know I can do this several ways but for conciseness in formula I want to use the double bracket assignment such as:

Ex[[ Sequence@@a ]] + Ey[[ Sequence@@b]]

In a large number of calculations however the sequence operation adds about 50% time overhead. Is there a way without the time cost overhead? thanks


The form of these two expressions are constructed elsewhere in the program

ind = {{i, j + 1, k}, {i - 1, j, k - 1}};

and need to be passed to a function, e.g.

f[in_] := Module[{a = in[[1]], b = in[[2]]},
    Table[Ex[[Sequence @@ a]] + Ey[[Sequence @@ b]], {m, 500}, {i, 2, 8}, 
      {j, 2, 8}, {k, 2, 8}]];

If I did not have the requirement to construct the index form somewhere else I could simply put:

g[]:= Table[
  Ex[[i, j + 1, k]] + Ey[[i - 1, j, k - 1]], {m, 500}, {i, 2, 8}, {j, 
 2,8}, {k, 2, 8}]

Ex = RandomReal[{0, 10}, {10, 10, 10, 10}];
Ey = RandomReal[{0, 10}, {10, 10, 10, 10}];

AbsoluteTiming[f[ind];]

{0.651337, Null}

AbsoluteTiming[g[];]

{0.0892097, Null}

$\endgroup$
  • $\begingroup$ You say "assignment" but I see only an addition. Is this a typo? $\endgroup$ – Henrik Schumacher Dec 5 '18 at 22:37
5
$\begingroup$

You may use Extract:

m = 100;
n = 10000;
Ex = RandomReal[{-1, 1}, {m, m, m}];
Ey = RandomReal[{-1, 1}, {m, m, m}];
a = RandomInteger[{1, m}, {n, 3}];
b = RandomInteger[{1, m}, {n, 3}];

r1 = Table[Ex[[Sequence @@ a[[i]]]] + Ey[[Sequence @@ b[[i]]]], {i, 1, Length[a]}]; // AbsoluteTiming // First
r2 = Extract[Ex, a] + Extract[Ey, b]; // AbsoluteTiming // First
r1 == r2

0.103876

0.000698

True


Edit

In the meantime, I got to the conclusion that the reason for f being slower than g is that Sequence is not compilable. So let us perform the compilation by hand:

f2[in_] := Quiet@Block[{X, XX, Y, YY},
    With[{XX = Part[X, Sequence @@ in[[1]]], YY = Part[Y, Sequence @@ in[[2]]]},
     Compile[{{X, _Real, 4}, {Y, _Real, 4}},
       Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
       CompilationTarget -> "WVM"
       ];
     cf[Ex, Ey]
     ]
    ];

Now we have:

aa = f[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
bb = g[]; // RepeatedTiming // First
cc = f2[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == bb == cc

0.692

0.090

0.0775

True

A further speed-up can be obtained by replacing Part by CompileGetElement`:

f3[in_] := Block[{X, XX, Y, YY},
   With[{XX = Compile`GetElement[X, Sequence @@ in[[1]]], YY = Compile`GetElement[Y, Sequence @@ in[[2]]]},
    Compile[{{X, _Real, 4}, {Y, _Real, 4}},
      Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
      CompilationTarget -> "WVM",
      RuntimeOptions -> "Speed"
      ];
    cf[Ex, Ey]
    ]
   ];

Runtime test:

dd = f3[{{i, j + 1, k}, {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == dd

0.028

True

Finally, let me say that handing over the code for the indices by the symbolic code {{i, j + 1, k}, {i - 1, j, k - 1}} is very fragile. A more robust way would be to submit functions that generate the symbolic code:

f4[in_] := Quiet@Block[{X, XX, Y, YY, i, j, k},
   With[{XX = Compile`GetElement[X, Sequence @@ in[[1]][i, j, k]], YY = Compile`GetElement[Y, Sequence @@ in[[2]][i, j, k]]},
    Compile[{{X, _Real, 4}, {Y, _Real, 4}},
      Table[XX + YY, {m, 1, 500}, {i, 2, 8}, {j, 2, 8}, {k, 2, 8}],
      CompilationTarget -> "WVM",
      RuntimeOptions -> "Speed"
      ];
    cf[Ex, Ey]
    ]
   ];

Test:

ee = f4[{{i, j, k} \[Function] {i, j + 1, k}, {i, j, k} \[Function] {i - 1, j, k - 1}}]; // RepeatedTiming // First
aa == ee

0.028

True

$\endgroup$
  • $\begingroup$ I have not explained myself well. The indices to be used are created $\endgroup$ – user8281 Dec 5 '18 at 23:12
  • $\begingroup$ I have not explained myself well: $\endgroup$ – user8281 Dec 5 '18 at 23:12
  • $\begingroup$ Then please clarify in your post. It may help to post the code that you actually used. $\endgroup$ – Henrik Schumacher Dec 5 '18 at 23:16
  • $\begingroup$ Thank you for your detailed explanation. It is very helpful $\endgroup$ – user8281 Dec 9 '18 at 19:15
  • $\begingroup$ You're welcome. $\endgroup$ – Henrik Schumacher Dec 9 '18 at 19:19

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.