Compile a program with Fold function

Question

I have the following small program:

(*Parameters*)
x0 = 4.;
g = 10.;
k = 0.5;
(*Define function*)    
fun[nTot_?IntegerQ, dt_?NumberQ] :=
 (fac = 1. - dt*k/g;
  list = Reap[Fold[Sow[fac*#] &, x0, Range[nTot]]][[2, 1]];
  list2 = ({(# - 1.)*dt, list[[#]]} & /@ Range[Length[list]]);)
(*Parameters*)
nTot = 10^7; dt = 10.^-3;
(*Computefunction*)
fun[nTot, dt]

which takes on my machine about 13 seconds. I want to speed it up and thought about Mathematica's Compile capabilities. Is there any way to Compile the function fun to C? For example here it is done for a Do loop, but I cannot apply this to my code. Any suggestions?

[Edit]

Thank you again for your answer. To illustrate what I mean see the following change of my initial program:

f[nTot_?IntegerQ, nMem_?IntegerQ, 
  dt_?NumberQ] :=
 (fun[t_] = -163.6 E^(-2. t) + 40. E^(-1. t);
  fac1 = 1. - dt*(41.8/0.1 + 1.);
  fac2 = dt^2/0.1;
  listExp = 
   Developer`ToPackedArray[
    fun[#*dt] & /@ Range[0, nMem - 1] // Reverse];
  list = Prepend[
    Reap[Fold[Append[Rest@#, Sow[fac1*Last[#] - fac2*listExp.#]] &, 
       ConstantArray[1., nMem], Range[nTot]]][[2, 1]], 1.];
  list2 = ({(# - 1.)*dt, list[[#]]} & /@ Range[Length[list]]);)
(*Define parameters*)
nTot = 10^4; nMem = 10^5; dt = 10.^-4;
(*Compute function*)
f[nTot, nMem, dt]

which takes 2.3 secs on my machine. And now the (equivalent) compiled version which takes 8.4 secs on my machine:

fun[t_] = -163.6 E^(-2.* t) + 40. E^(-1.* t);
cfun = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {nMem, \
_Integer}}, 
   Block[{}, x = ConstantArray[x0, nMem]; 
    fac1 = 1. - dt*(41.8/0.1 + 1.);
    fac2 = dt^2/0.1; 
    listExp = Reverse[fun[#*dt] & /@ Range[0, nMem - 1]];
    Table[{i*dt, 
      Last[x = Append[Rest@x, fac1*Last[x] - fac2*(listExp.x)]]}, {i, 
      1, nTot}]], CompilationTarget -> "C"];
clist2 = Prepend[cfun[1., 10^4, 10.^-4, 10^5], {0., 1.}];

Why is the compiled version now so slow? Is there a way to speed it up? Or did I do some simple mistakes? Help is highly appreciated again!

[Edit2]:

Ok, a small modification in the compiled version speeds it up to around 3 secs on my machine. But still it is slower compared to the solution with Fold:

fac1 = 1. - dt*(41.8/0.1 + 1.);
fac2 = dt^2/0.1;
fun[t_] = -163.6 E^(-2.* t) + 40. E^(-1.* t);
cfun = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {nMem, \
_Integer}}, Block[{fac2 = fac2, fac1 = fac1, x, listExp},
    x = Table[x0, {i, 1, 10^5}];
    listExp = Reverse[fun[#*dt] & /@ Range[0, nMem - 1]];
    Table[{i*dt, 
      Last[x = Append[Rest@x, fac1*Last[x] - fac2*listExp.x]]}, {i, 1,
       nTot}]], {{listExp, _Real, 1}}, 
   CompilationOptions -> {"InlineExternalDefinitions" -> True, 
     "InlineCompiledFunctions" -> True, 
     "ExpressionOptimization" -> True}, 
   RuntimeAttributes -> {Listable}, Parallelization -> True];
clist2 = Prepend[cfun[1, 10^4, 10^-4, 10^5], {0, 1}];

[Edit3]

Getting rid of Append and using Part instead is a little bit helpful:

dt = 10.^-4;
nMem = 10^5;
nTot = 10^5;
fac1 = 1. - dt*(41.8/0.1 + 1.);
fac2 = dt^2/0.1;
x0 = 1.;
x = Array[x0 &, nTot + nMem];
fun[t_] = -163.6 E^(-2.* t) + 40. E^(-1.* t);
listExp = 
  Developer`ToPackedArray[Reverse[fun[#*dt] & /@ Range[0, nMem - 1]]];
list = Table[{i*dt, 
    x[[i + nMem]] = 
     fac1*x[[i + nMem - 1]] - 
      fac2*listExp.x[[i ;; nMem + i - 1]]}, {i, 1, nTot}];

and a compiled version:

fac1 = 1. - dt*(41.8/0.1 + 1.);
fac2 = dt^2/0.1;
fun[t_] = -163.6 E^(-2.* t) + 40. E^(-1.* t);
cfun = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {nMem, \
_Integer}}, Block[{fac2 = fac2, fac1 = fac1, x, listExp},
    x = Array[x0 &, nTot + nMem];
    listExp = Reverse[fun[#*dt] & /@ Range[0, nMem - 1]];
    Table[{i*dt, 
      x[[i + nMem]] = 
       fac1*x[[i + nMem - 1]] - 
        fac2*listExp.x[[i ;; nMem + i - 1]]}, {i, 1, 
      nTot}]], {{listExp, _Real, 1}}, 
   CompilationOptions -> {"InlineExternalDefinitions" -> True, 
     "InlineCompiledFunctions" -> True, 
     "ExpressionOptimization" -> True}, 
   RuntimeAttributes -> {Listable}, Parallelization -> True];
clist2 = cfun[1, 10^5, 10^-4, 10^5];

[Edit4]

Managed to get rid of the MainEvaluate code by using pure function definition:

Clear["Global`*"]
fac1 = 1. - dt*(41.8/0.1 + 1.);
fac2 = dt^2/0.1;
fun = -163.6 E^(-2.* #) + 40. E^(-1.* #) &;
cfun = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {nMem, \
_Integer}}, Block[{fac2 = fac2, fac1 = fac1, x, listExp},
    x = Array[x0 &, nTot + nMem];
    listExp = Reverse[fun[#*dt] & /@ Range[0, nMem - 1]];
    Table[{i*dt, 
      x[[i + nMem]] = 
       fac1*x[[i + nMem - 1]] - 
        fac2*listExp.x[[i ;; nMem + i - 1]]}, {i, 1, 
      nTot}]], {{listExp, _Real, 1}}, 
   CompilationOptions -> {"InlineExternalDefinitions" -> True, 
     "InlineCompiledFunctions" -> True, 
     "ExpressionOptimization" -> True}];
clist3 = cfun[1., 10^4, 10.^-4.5, 4*10^5];

Any further recommendations how to speed it up?

You want to avoid MainEvaluate. See the output of Needs["CompiledFunctionTools`"]; CompilePrint@cfun. — Michael E2, Jan 15 '18 at 14:58
Developer`ToPackedArray is basically an expensive no-op inside Compile, because the only arrays Compile deals with are packed arrays. — Michael E2, Jan 15 '18 at 15:00
Can you elaborate a little bit more on how to avoid MainEvaluate in my specific case? I inserted these commands but what does the output exactly tells me? — NeverMind, Jan 15 '18 at 15:06
Ok, so DeveloperToPackedArrayis basically a useless function withinCompile` if I get you right. But why is there an effect in my program? — NeverMind, Jan 15 '18 at 15:07
Yep, ToPackedArray is useless inside Compile. MainEvaluate is a call out of the compiler (actually WVM) back to the main kernel to evaluate something the compiler cannot do. For more explanation, explore some of the Q&A in the link. — Michael E2, Jan 15 '18 at 15:17
@Micheal E2: I read about pre-defining constants and use them as a=a in the With[] environment. But how can I do this with a vector like x or listExp? I don't see another way to avoid MainEvaluate. @Henrik Schumacher: Your first suggestion doesn't seem to have a huge effect. Your second tip is worth a try, but I think I will end up with Part then, and I don't know whether it is really faster to keep the very long list. I actually can predict the length of the output list since it shall be pre-defined by nTot. — NeverMind, Jan 15 '18 at 15:37
External constants can be injected with With or passed as arguments. For variables like x, localize them inside Block or Module (e.g. Block[{x, listExp}, <code>], as Henrik did in his answer). — Michael E2, Jan 15 '18 at 15:49
Alright, there is one MainEvaluate left, related to R8 = MainEvaluate[ Hold[fun][ R6]]. How to treat this function fun[t]. Cannot be localized within Block, can it? — NeverMind, Jan 15 '18 at 16:21
@Henrik Schumacher: In my [Edit3] I got rid of Append and used indexing of a pre-defined list instead. I used Array[] here to create this list since according to Compile\CompilerFunctions[]` it is a compilable function. Still, I think one can improve the performance. If you any further idea let me please know. — NeverMind, Jan 15 '18 at 18:53
@Michael E2: Managed now to get rid of the MainEvaluate by using a pure function definition. See my [Edit4] in the OP. However, it still is very slow unfortunately. So far your tips and the ones from Henrik Schumacher have been very helpful, so I wonder whether there is more I can do to improve the overall performance of the code. — NeverMind, Jan 15 '18 at 23:47

Henrik Schumacher · Answer 1 · 2018-01-15T10:19:02.683

7

Your use of Reap and Sow along with Fold could be better formulated with FoldList, or even better with NestList. Moreover, the second sweep over Range[Length[list]] can be avoided at all. "Unfolding" your code and compiling it, this could look like this:

cfun = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {g, _Real}, {k, _Real}},
   Block[{fac = 1. - dt*k/g, x = x0},
    Table[x *= fac; {(i - 1.)*dt, x}, {i, 1, nTot}]
    ],
   CompilationTarget -> "C"
   ];

clist2 = cfun[x0, nTot, dt, g, k]; // AbsoluteTiming // First

0.49941

Note that constructs like Fold and Nest are bound to evaluate sequentially while your result can actually obtained in parallel. A simple vectorized version of your code is this:

vfun[x0_?NumberQ, nTot_Integer, dt_?NumberQ, g_?NumberQ, k_?NumberQ] :=
   With[{ran = N@Range[1, nTot]},
   Transpose[{(ran - 1.) dt, x0 ((1. - dt*k/g)^ran)}]
   ];

On my Haswell quad core, this performs like this

vlist2 = vfun[x0, nTot, dt, g, k]; // AbsoluteTiming // First

0.397826

Note that the most expensive operation in vfun is Transpose; skipping will half the cost.

Comparing the results:

Max[Abs[clist2 - vlist2]]

2.79776*10^-14

edited Jan 15 '18 at 10:19

answered Jan 15 '18 at 10:11

Henrik Schumacher

106,770
7
179
309

Thank you very much for your recommendations! For some reasons, vfun is becoming slow if nTot=2*10^7 (25 secs) compared to cfun with the same nTot (2 secs) on my machine. What could be reason for this? – NeverMind Jan 15 '18 at 10:33
When the exponents become too large, Mathematica has to use extended precision as the results get much smaller than $MachineEpsilon. Extended precision calculations are much more expensive since they have to be emulated in software and since bigger data types have to be used. The compiled version does not care; it will simply return 0. in that case. – Henrik Schumacher Jan 15 '18 at 10:39
Of course, makes totally sense. Can one avoid this behaviour? If yes, how? – NeverMind Jan 15 '18 at 10:40
1

Ok, SetSystemOptions["CatchMachineUnderflow" -> False] should do the job. – NeverMind Jan 15 '18 at 10:48
Good point. Yeah, a thought about it in the meantime: I think Mathematica cannot switch to extended precision because (1. - dt*k/g) has machine precision. So, it's the error handling that slows things down. Very good! – Henrik Schumacher Jan 15 '18 at 10:50
Another observation is that skipping Transpose does not really half the cost. At least not after setting the "CatchMachineUnderflow" to False. For nTot=10^9 timing is almost identical.
One further question about the vectorized code: Here, the next element only depends on the previous one. And it is easy to see that vectorization is possible.

But what about the case, where I have memory and the next element depends on a (infinite) list of previous elements. In this case vectorization is probably not possible anymore, right? Hope you can understand what I mean.
– NeverMind Jan 15 '18 at 10:56
Hm. I cannot confirm that. However, these things may depend quite much on the specific architecture, e.g. on Flops/s vs. Byte/s (transposition is involves copying and thus is memory bound). My machine is quite good at Flops/s but the memory is rather old (due to the Haswell architecture I use). – Henrik Schumacher Jan 15 '18 at 11:01
In principle: whenever the next element cannot be computed without knowning the previous element, parallelization cannot be used. (That's not entirely true as there are such algorithms as Parareal). But in your case, we are able to predict the value of an arbitrary element, without knowing any other. – Henrik Schumacher Jan 15 '18 at 11:05
Thank you again for your comment. Please see my edit in the opening post of this topic. – NeverMind Jan 15 '18 at 12:18
I added a minimal working example in my [Edit] in the opening post. For some reason the compiled version seems to be slower. Is there a way to speed it up? [Edit] It seems to be that slow because in the compiled version I cannot use a PackedArray for the listExp. Or is there a way to do it? – NeverMind Jan 15 '18 at 14:29
I may have a look at that tonight (in 5 hours or so...) – Henrik Schumacher Jan 15 '18 at 15:11
Slightly faster vfun: With[{fac = 1. - dt*k/g, r = Range[0, nTot - 1]}, Transpose@{dt*r, Exp[Log[x0*fac] + Log[fac] r]}] – Michael E2 Jan 16 '18 at 04:05

Henrik Schumacher · Accepted Answer · 2018-01-17T12:43:07.037

This new task with dotting listExp against a moving window in x would have deserved a new question post.

Your problem as a very specific structure that I am going to exploit: The list listExp is a sum of two lists listExp1 and listExp2 that are both of the form Table[a base^i, {i,1,nMem - 1}]. These list are dotted against a window moving along x resulting in numbers sum1 and sum2. The main obserservation for speeding up your code is that after a moving step, sum1 can be computed as sum1 = (sum1 + a1 xout)/base1 + b1 xin where xin is the element of x entering the window and xout is the one leaving the window. Hence, we don't have to compute the expensive dot product any more! Along with some other tweaks (replacing Part in read operations by Compile`GetElement, compiling into binaries via "C", and switching off some secrity checks with RuntimeOptions -> "Speed"), the code looks like this:

cfunDotLess = Compile[{{x0, _Real}, {nTot, _Integer}, {dt, _Real}, {nMem, _Integer}}, 
   Block[{fac2, fac1, x, listExp1, listExp2, base1, base2, sum1, sum2,
      xnew, a1, a2, b1, b2, xout, xin, j}, 
    x = Table[x0, {i, 1, nMem}];
    fac1 = 1. - dt (41.8/0.1 + 1.);
    fac2 = dt^2/0.1;
    base1 = E^(-2. dt);
    base2 = E^(-1. dt);
    a1 = -163.6 fac2;
    a2 = 40. fac2;
    sum1 = 0.;
    sum2 = 0.;
    b1 = (-a1);
    b2 = (-a2);
    Do[
     sum1 += b1 x0;
     sum2 += b2 x0;
     b1 *= base1;
     b2 *= base2;, {i, 0, nMem - 1}];
    xnew = x0;
    j = 0;
    Table[
     If[j >= nMem, j = 1, j++];
     xout = Compile`GetElement[x, j];
     xin = xnew;
     sum1 = (sum1 + a1 xout)/base1 + b1 xin;
     sum2 = (sum2 + a2 xout)/base2 + b2 xin;
     x[[j]] = xnew = Plus[fac1 xin, sum1, sum2];
     {i dt, xnew}, {i, 1, nTot}]
    ],
   CompilationTarget -> "C",
   RuntimeOptions -> "Speed"
   ];

The following shows that clistDotLess is about 8000 times faster than cfun from Edit 4 at a cost of an acceptable loss of accuracy:

clistDotLess = cfunDotLess[1., 10^5, 10.^-4.5, 4.*10^5]; // AbsoluteTiming // First
clist = cfun[1., 10^5, 10.^-4.5, 4.*10^5]; // AbsoluteTiming // First
Max[Abs[clistDotLess - clist]]

0.003614

37.028

4.63854*10^-9

Edit: I removed some superfluous indexing into x and shortened x to be an array of length nMem by using indices j modulo nMem.

First of all, thank you very much for your effort! I wasn't aware that I should open a new topic for this one, but I will certainly do next time. I realized that in my [Edit4] I somehow lost the Reverse in the definition of listExp. I re-edited my [Edit4] to provide the correct definition as intended. I think this won't affect your answer too much, but feel free to edit it so that it is fitting to the opening post again. I will thouroughly look at your answer tomorrow, but if it is 8000 times faster than cfun it is exactly what I was looking and hoping for. Great! — NeverMind, Jan 16 '18 at 20:36
What I do not understand is how you obtain the expression for the new sum after a moving step. Shouldn't it be sum=base*sum+a*xin-a*base^(nMem+1)*xout? — NeverMind, Jan 17 '18 at 12:16
Maybe. Remember that you reverted the ordering in listExp. Division by base might become multiplication and so on. Quite likely, I also made some mistakes. I am sure you will adapt this to you problem... — Henrik Schumacher, Jan 17 '18 at 12:31
Sure, I think this will just take some time. But more importantly I got the idea behind it. — NeverMind, Jan 17 '18 at 12:33
Adapted it now to my problem. Now, it is 'only' 5000 faster than the previous solution with cfun but that is fast enough anyways for my purposes. Again, thank you so much for your help. If I could give you more than one up vote I would certainly do. — NeverMind, Jan 17 '18 at 13:45

Compile a program with Fold function

2 Answers2