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