# Applying f for a part of list (MapAt is slow)

Suppose that there are 2000 stones in a row.
You should grind 1000 of them to make beads.
But which 1000 stones?

case1) the first 1000 stones.

case2) Let B = RandomSample[Range[2000], 1000].
Then, all b-th stones, where b denotes element of B.

If this is a problem of real world, then case1) is faster then case2), because
for case2),

1. You have to read the B every time after grinding a stone.
2. If the next stone to be processed is far away, it will take a long time to go there.
3. B is even not sorted by size.

But, if it was a problem of mathematica,

2000 stones in a row = a list of 2000 elements
grinding a stone = applying function to an element

The time loads for case1) and time loads for case2) should be similar in my opinion.

I experienced that "just because some elements of a list are positionally close, there is no particular advantage to processing(=applying function) them at once(or consecutively)"

Also I believe that "there is a fast mathematica code for case2) always, that takes silmilar time compared to case1)"

Am I thinking correctly?

If so, can you make a fast code for the following problem ?

For L, square(=apply #^2&) each b-th element of L.(leave other element unchanged)

L = Range[40000, 59999];
B = RandomSample[Range[20000], 10000];


My first trial was

MapAt[#^2&, L, {#}&/@B] //Timing
2.85938


which is never successful.

After some other trials, I realized that MapAt is very slow, in all cases.
The performance was so bad that I felt MapAt should not be used in any case.

Finally I succeeded to shortening the time(using ReplacePart) but it is not neat and I don't think my code is as fast as professional programmer's one.

What skill do you use, when applying a function for a part of list? (faster!)

Can you help me?

• Much faster: cres = lx; cres[[bx]] = lx[[bx]]*lx[[bx]].
– Alan
Mar 5, 2022 at 18:55

I am sure this is not the fastest solution but it is 25 times faster.

lx = Range[40000, 59999];
bx = RandomSample[Range[20000], 10000];
ares = MapAt[#^2 &, lx, {#} & /@ bx] // Timing
bres = SubsetMap[#^2 &, lx, {#} & /@ bx] // Timing

ares[[2]] == bres[[2]]

(*True*)

{First@ares, First@bres}

{1.17001, 0.0468003}

• Wow! SubsetMap is a new feature in V12! Thank you! Mar 5, 2022 at 17:48
• For this particular variation of the SubsetMap syntax, you are executing a MapAt at the subset of the positions that you have chosen somehow, so there is no mental gymnastics involved.
– Syed
Mar 5, 2022 at 17:51
• May I ask what if you were in mathematica V11 ? Mar 5, 2022 at 17:59
• Please edit this requirement in to your post. I don't have M11 so I cannot say for sure, but if I find another way I will update the answer.
– Syed
Mar 5, 2022 at 18:17

Here are a couple more to add to Syed's & Alan's. The compiled version is much faster, if everything stays machine-size. If you're doing this repeatedly, you could compile ahead of time (to C even). You have to recompile if you change the function from #^2; and if the function to be applied is not compilable, then it obviously won't work well.

lx = Range[40000, 59999];
bx = RandomSample[Range[20000], 10000];
ares = MapAt[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
bres = SubsetMap[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
(cres = lx; Do[cres[[k]] = cres[[k]]^2, {k, bx}]); // AbsoluteTiming
dres = Compile[{{lx, _Integer, 1}, {bx, _Integer, 1}},
Block[{dres = lx},
Do[dres[[k]] = dres[[k]]^2, {k, bx}];
dres
]
][lx, bx]; // AbsoluteTiming
(eres = lx; eres[[bx]] = eres[[bx]]^2); // AbsoluteTiming

ares == bres == cres == dres == eres
(*
{0.633558, Null}
{0.013454, Null}
{0.006371, Null}
{0.000657, Null}
{0.000221, Null}  <-- Alan's is fastest
True
*)

• Note @Alan's depends on B or bx not having duplicate indices, which it doesn't in the problem given in the OP. Mar 5, 2022 at 19:34
• FunctionCompile beats them all, but there is an initial long compilation delay. I have posted the code as an extension of this answer. Mar 6, 2022 at 4:55

I just extended @Michael-E2 answer by including FunctionCompile:

lx = Range[40000, 59999];
bx = RandomSample[Range[20000], 10000];
ares = MapAt[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
bres = SubsetMap[#^2 &, lx, {#} & /@ bx]; // AbsoluteTiming
(cres = lx; Do[cres[[k]] = cres[[k]]^2, {k, bx}]); // AbsoluteTiming
dres = Compile[{{lx, _Integer, 1}, {bx, _Integer, 1}},
Block[{dres = lx}, Do[dres[[k]] = dres[[k]]^2, {k, bx}];
dres]][lx, bx]; // AbsoluteTiming
(eres = lx; eres[[bx]] = eres[[bx]]^2); // AbsoluteTiming
ffun = FunctionCompile[
Function[
{Typed[lx, TypeSpecifier["PackedArray"]["MachineInteger", 1]],
Typed[bx, TypeSpecifier["PackedArray"]["MachineInteger", 1]]},
Block[{fres = lx}, Do[fres[[k]] = fres[[k]]^2, {k, bx}]; fres]]];
fres = ffun[lx, bx]; // AbsoluteTiming

ares == bres == cres == dres == eres == fres

(*
{1.29098,Null}
{0.0542107,Null}
{0.0229712,Null}
{0.0005062,Null}
{0.0001768,Null}
{0.0000713,Null}
True
*)

• +1 because updating fres is so fast. Nice. I wonder how much the compile time should be counted. (The compile time is included in dres, but it's fairly negligible.) Mar 6, 2022 at 5:30

Another approach that is supported by the previous versions (11.3 at least) is to create a list of functions which are either (Identity which doesn't touch the data or fn which apply your function) and then MapThread this list with Construct:

Block[{temp},

(* create a list of Identity function *)
temp = ConstantArray[Identity, Length[L]];

(* replace the positions with your function *)
temp[[B]] = #^2 &;

] // MaxMemoryUsed // AbsoluteTiming

(* Out: {0.0161858, 800568} *)


Result:

result == SubsetMap[#^2 &, L, {#} & /@ B]
(* Out: True *)


Another benefit is the memory footprint (could vary):

MapAt[#^2 &, L, {#} & /@ B]; // MaxMemoryUsed // AbsoluteTiming
(* Out: {1.2336, 1912832} *)


which is around 2 and 8 times smaller respectively.
• Minimum version 11.3 is because of Construct, if you happen to have an alternative, the minimum requirement would be lowered (although performance would be affected).