Re: Help on compiling a function
- To: mathgroup at smc.vnet.net
- Subject: [mg115647] Re: Help on compiling a function
- From: "Sjoerd C. de Vries" <sjoerd.c.devries at gmail.com>
- Date: Mon, 17 Jan 2011 05:39:51 -0500 (EST)
- References: <iguiqe$8sv$1@smc.vnet.net>
Hi Ramiro,
One trick would be to work with the Logs of the numbers. There happens
to be a Log of the Gamma function on board, LogGamma, that is able to
deal with large numbers. Below the uncompiled function example0, your
version (example1), my first approach (example2), and the same with
compilationtarget set to C (example3). You can speed up by a factor of
7.
example0[n_, a_, b_, t_] :=
Gamma[Total[n] + a]/(Times @@ (Gamma[n + 1])*Gamma[a])*b^a*
Times @@ (t^n)/(Total[t] + b)^(Total[n] + a);
example1 =
Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
Gamma[Total[n] + a]/(Times @@ (Gamma[n + 1])*Gamma[a])*b^a*
Times @@ (t^n)/(Total[t] + b)^(Total[n] + a)];
example2 =
Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
Exp[LogGamma[
Total[n] + a] - (Plus @@ (LogGamma[n + 1]) + LogGamma[a]) +
Log[a b] +
Plus @@ (n Log[(t)]) - (Total[n] + a) Log[(Total[t] + b)]]]
example3 =
Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
Exp[LogGamma[
Total[n] + a] - (Plus @@ (LogGamma[n + 1]) + LogGamma[a]) +
Log[a b] +
Plus @@ (n Log[(t)]) - (Total[n] + a) Log[(Total[t] + b)]],
CompilationTarget -> "C"]
In[119]:=
x0 = Table[
example0[{97.6203, 8.4788, 21.4204, 46.1755 + i}, 1,
1, {39.9342, 7.5820, 5.8656, 10.0553}], {i, 0, 1, 0.00001}] //
Timing;
x1 = Table[
example1[{97.6203, 8.4788, 21.4204, 46.1755 + i}, 1,
1, {39.9342, 7.5820, 5.8656, 10.0553}], {i, 0, 1, 0.00001}] //
Timing;
x2 = Table[
example2[{97.6203, 8.4788, 21.4204, 46.1755 + i}, 1,
1, {39.9342, 7.5820, 5.8656, 10.0553}], {i, 0, 1, 0.00001}] //
Timing;
x3 = Table[
example2[{97.6203, 8.4788, 21.4204, 46.1755 + i}, 1,
1, {39.9342, 7.5820, 5.8656, 10.0553}], {i, 0, 1, 0.00001}] //
Timing;
Speed up factors (native->yours, native->Compile with Log, native->
Compile with Log + C):
In[123]:=
{First[x0]/First[x1], First[x0]/First[x2], First[x0]/First[x3]}
Out[123]= {0.949367, 7.17735, 7.27961}
The results are extremely similar but -due to numerical errors- not
identical. The numerical errors are small; on average they relative
difference is:
In[124]:= Mean[2 (x1[[2]] - x2[[2]])/(x1[[2]] + x2[[2]])]
Out[124]= -4.79378*10^-14
Cheers -- Sjoerd
On Jan 16, 11:55 am, Ramiro <ramiro.barran... at gmail.com> wrote:
> Hi everyone,
>
> I had written about this but wanted to revisit the problem in light of
> the latest Mathematica options for compiling (C and parallelization),
> which I was hoping to try to use on this function, which is the main
> arithmetic and time-consuming part of my simulation.
>
> I have the following code:
>
> example1 =
> Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
> Gamma[Total[n] + a]/(Times @@ (Gamma[n + 1])*Gamma[a])*b^a*
> Times @@ (t^n)/(Total[t] + b)^(Total[n] + a)];
>
> but under input such as the following it breaks:
>
> example1[{97.6203, 8.4788, 21.4204, 46.1755}, 1, 1, {39.9342, 7.5820,
> 5.8656, 10.0553}]
> CompiledFunction::cfse: Compiled expression
> 1.33128164105722332870399207`12.920368310128136*^315 should be a
> machine-size real number. >>
> CompiledFunction::cfex: Could not complete external evaluation at
> instruction 4; proceeding with uncompiled evaluation. >>
>
> the problem is just that in the calculation Gamma ends up being really
> big, larger than $MaxMachineNumber so it complains.
>
> Anybody see a way around this? This is the workhorse function of my
> simulation, so any gain in speed helps a lot
>
> Thanks in advance,
> Ramiro
>
> Longer explanation (if interested):
>
> example =
> Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
> Module[{totaln, totalt, gammanp1times, tpowntimes, times1, div1,
> times2, times3, pow1, gammatotnpa}, totaln = Total[n];
> totalt = Total[t];
> gammanp1times = Apply[Times, Gamma[n + 1]];
> tpowntimes = Apply[Times, t^n];
> times1 = Times[gammanp1times, Gamma[a]];
> gammatotnpa = Gamma[totaln + a];
> div1 = Divide[gammatotnpa, times1];
> times2 = Times[div1, b^a];
> times3 = Times[times2, tpowntimes];
> pow1 = Power[totalt + b, totaln + a];
> Divide[times3, pow1]]
> ]
> example[{97.6203, 8.4788, 21.4204, 46.1755}, 1, 1, {39.9342, 7.5820,
> 5.8656, 10.0553}]
>
> the result of gammatotnpa is greater than 10^315 which is greater than
> $MaxMachineNumber on my machine. Any suggestions?