Re: Help on compiling a function
- To: mathgroup at smc.vnet.net
- Subject: [mg115667] Re: Help on compiling a function
- From: Ramiro Barrantes-Reynolds <Ramiro.Barrantes at uvm.edu>
- Date: Tue, 18 Jan 2011 05:47:24 -0500 (EST)
Thanks! This is exactly what I needed!
Quoting Daniel Lichtblau <danl at wolfram.com>:
>
> ----- Original Message -----
>> From: "Ramiro" <ramiro.barrantes at gmail.com>
>> To: mathgroup at smc.vnet.net
>> Sent: Sunday, January 16, 2011 4:54:33 AM
>> Subject: [mg115614] Help on compiling a function
>> 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):
>> [...]
>
> You have products/quotients with intermediate overflow because some
> factors get too big before division brings them back down. The
> easiest way to avoid this is to use logarithms, sum, then
> exponentiate. This is workable because your factors are all explicit
> powers of positives or Gamma function evaluations, so you can use
> the power-of-log rule and LogGamma.
>
> In[285]:== example1[{97.6203, 8.4788, 21.4204,
> 46.1755}, 1, 1, {39.9342, 7.5820, 5.8656, 10.0553}]
>
> During evaluation of In[285]:== CompiledFunction::cfex: Could not
> complete external evaluation at instruction 3; proceeding with
> uncompiled evaluation. >>
>
> Out[285]== 1.062645560684518*10^-11
>
> In[286]:== example2 ==
> Compile[{{n, _Real, 1}, {a, _Real}, {b, _Real}, {t, _Real, 1}},
> Module[{log1, log2, log3, log4, log5},
> log1 == LogGamma[Total[n] + a];
> log2 == Plus @@ LogGamma[n + 1];
> log3 == LogGamma[a];
> log4 == Plus @@ (n*Log[t]);
> log5 == (Total[n] + a)*Log[Total[t] + b];
> Exp[log1 - (log2 + log3) + log4 - log5]*b^a]]
>
> In[287]:== example2[{97.6203, 8.4788, 21.4204,
> 46.1755}, 1, 1, {39.9342, 7.5820, 5.8656, 10.0553}]
>
> Out[287]== 1.062645560684574*10^-11
>
> Daniel Lichtblau
> Wolfram Research
>
--
Ramiro Barrantes-Reynolds
Ph.D. Candidate. Computational Cell and Molecular Biology
Microbiology and Molecular Genetics
University of Vermont