6
$\begingroup$

These days I am trying to rewrite the following recursive function, in order to gain some speed. It is crucial because I end up using 4-5 similar functions, and I need to decrease execution time.

This is my original function

Clear[function]
difference = 0;
remaining = 80;


p6h = 0.0171;
p4h = 0.0171;
p2h = 0.00857;
p1h = 0.001906;

p6a = 0.02446;
p4a = 0.0244;
p2a = 0.012234;
p1a = 0.00271;

pns = 1 - p6h - p4h - p2h - p1h - p1a - p2a - p4a - p6a;


function[t_, d_] := 
  function[t, d] = 
   If[t > 80, 0, 
    If[d == difference && t == remaining, 1, 
     p6h*function[t + 1, d - 6] + p4h*function[t + 1, d - 4] + 
      p2h*function[t + 1, d - 2] + p1h*function[t + 1, d - 1] + 
      p6a*function[t + 1, d + 6] + p4a*function[t + 1, d + 4] + 
      p2a*function[t + 1, d + 2] + p1a*function[t + 1, d + 1] + 
      pns*function[t + 1, d]]];

Once function is defined I seek several summations similar to these (below my timings)

function[0, 0] // AbsoluteTiming
Sum[function[0, i], {i, 1, 560}] // AbsoluteTiming
Sum[function[0, -i], {i, 1, 560}] // AbsoluteTiming
Sum[function[0, i], {i, -560, 560}] // AbsoluteTiming

(*
{0.589095, 0.0399939}    
{0.679503, 0.293243}    
{0.710139, 0.666763}    
{0.00101605, 1.} *)

I searched stackexchange and found several topics discussing about dynamic programming and recursive function. I even tried to replicate some methods provided by @LeonidShifrin but either I fail, or I end up with slower evaluation times.

Based on his answer here: i wrote this

  functionB[t_, d_] :=
 Block[{functionB},
  functionB[t1_, d1_] := 
   functionB[t1, d1] = 
    If[t1 > 80, 0, 
     If[d1 == difference && t1 == remaining, 1, 
      p6h*functionB[t1 + 1, d1 - 6] + p4h*functionB[t1 + 1, d1 - 4] + 
       p2h*functionB[t1 + 1, d1 - 2] + p1h*functionB[t1 + 1, d1 - 1] +
        p6a*functionB[t1 + 1, d1 + 6] + 
       p4a*functionB[t1 + 1, d1 + 4] + p2a*functionB[t1 + 1, d1 + 2] +
        p1a*functionB[t1 + 1, d1 + 1] + pns*functionB[t1 + 1, d1]]];
  functionB[t, d]
  ]

but it is dramatically slow.

I tried this approach too

Clear[functionB2, "functionB2`*"]

Begin["functionB2`"];


functionB2[t1_, d1_] := 
 functionB2[t1, d1] = 
  Expand[If[t1 > 80, 0, 
    If[d1 == difference && t1 == remaining, 1, 
     p6h*functionB2[t1 + 1, d1 - 6] + p4h*functionB2[t1 + 1, d1 - 4] +
       p2h*functionB2[t1 + 1, d1 - 2] + 
      p1h*functionB2[t1 + 1, d1 - 1] + 
      p6a*functionB2[t1 + 1, d1 + 6] + 
      p4a*functionB2[t1 + 1, d1 + 4] + 
      p2a*functionB2[t1 + 1, d1 + 2] + 
      p1a*functionB2[t1 + 1, d1 + 1] + pns*functionB2[t1 + 1, d1]]]]



End[];

but times are slightly worse than original approach.

(*
{0.641918, 0.0399939}
{0.738379, 0.293243}
{0.747237, 0.666763}
{0.00101013, 1.}
*)

Based on a previous question here and here I tried to compile it, but unfortunately I found out that it is not as easy as I thought (or even possible).

Finally, I though about setting all values calculated into a table, like this

> s = Table[function[0, i], {i, -560, 560}];

in order to add table's elements accordingly, only to find out that it is actually the same thing with memoization.

Could anyone help me rewrite the original function in order to gain speed?


update 18/3/2016

based on this usage of Block suggested here I again failed to improve speeds. I actually made it terribly slow once again

Clear[function]
difference = 0;
remaining = 80;
p6h = 0.0171;
p4h = 0.0171;
p2h = 0.00857;
p1h = 0.001906;
p6a = 0.02446;
p4a = 0.0244;
p2a = 0.012234;
p1a = 0.00271;
pns = 1 - p6h - p4h - p2h - p1h - p1a - p2a - p4a - p6a;

ClearAll[defWrap];
SetAttributes[defWrap, HoldFirst];
defWrap[fcall_] := Block[{function},
   function[t_, d_] := 
    function[t, d] = 
     If[t > 80, 0, 
      If[d == difference && t == remaining, 1, 
       p6h*function[t + 1, d - 6] + p4h*function[t + 1, d - 4] + 
        p2h*function[t + 1, d - 2] + p1h*function[t + 1, d - 1] + 
        p6a*function[t + 1, d + 6] + p4a*function[t + 1, d + 4] + 
        p2a*function[t + 1, d + 2] + p1a*function[t + 1, d + 1] + 
        pns*function[t + 1, d]]];
   fcall];

ClearAll[function];
function[tt_, kk_] := defWrap[function[tt, kk]];

Sum[function[0, i], {i, 1, 10}] // AbsoluteTiming
(* {6.28821, 0.197245} *)
$\endgroup$

Your Answer

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

Browse other questions tagged or ask your own question.