ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 13 Mar 2011 18:15:47 -0500convert mpmath function to cythonhttp://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.
Here is the function:
from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
return num/denom
def myfisher_mp(a1,a2,b1,b2):
p = hyp_mp(a1,a2,b1,b2)
outp = p
for i in range(0,a1):
temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
if temp > p:
break
else:
outp += temp
for i in range(0,a2):
temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
if temp > p:
break
else:
outp += temp
return outp
A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.Fri, 04 Mar 2011 08:39:58 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/Answer by DSM for <p>I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.</p>
<p>Here is the function:</p>
<pre><code>from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
return num/denom
def myfisher_mp(a1,a2,b1,b2):
p = hyp_mp(a1,a2,b1,b2)
outp = p
for i in range(0,a1):
temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
if temp > p:
break
else:
outp += temp
for i in range(0,a2):
temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
if temp > p:
break
else:
outp += temp
return outp
</code></pre>
<p>A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12170#post-id-12170Assuming you mean convert this function to Cython [ah, you do, it's only the title that's wrong-- edited], I'm not sure whether it will help.
First, I get
sage: time z = myfisher_mp(1286, 9548, 133437, 148905)
CPU times: user 1.56 s, sys: 0.00 s, total: 1.56 s
Wall time: 1.60 s
sage: z
4.191535038879969055586166e-1316
Is that right? Almost all of the computations in the loops over i in the function don't matter to the final value because they're too small.
Second, you only get real benefits from Cython when the work can be pushed into C. If the values are so small, though, C floats won't work, as they'll underflow. And it's not that most of the time is being spent in the Python loops, which are linear anyway: the vast majority of the time is being spent in the gamma function itself, which is already pretty fast.
sage: timeit('z=mpmath.gamma(100)')
625 loops, best of 3: 5.67 µs per loop
sage: timeit('z=mpmath.gamma(100000)')
625 loops, best of 3: 50.9 µs per loop
Do you have a link to the definition of this function? I'm sure we can find a way to compute it more efficiently, especially given how many terms are currently noncontributing.Fri, 04 Mar 2011 12:04:47 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12170#post-id-12170Comment by mhampton for <p>Assuming you mean convert this function to Cython [ah, you do, it's only the title that's wrong-- edited], I'm not sure whether it will help.</p>
<p>First, I get</p>
<pre><code>sage: time z = myfisher_mp(1286, 9548, 133437, 148905)
CPU times: user 1.56 s, sys: 0.00 s, total: 1.56 s
Wall time: 1.60 s
sage: z
4.191535038879969055586166e-1316
</code></pre>
<p>Is that right? Almost all of the computations in the loops over i in the function don't matter to the final value because they're too small.</p>
<p>Second, you only get real benefits from Cython when the work can be pushed into C. If the values are so small, though, C floats won't work, as they'll underflow. And it's not that most of the time is being spent in the Python loops, which are linear anyway: the vast majority of the time is being spent in the gamma function itself, which is already pretty fast.</p>
<pre><code>sage: timeit('z=mpmath.gamma(100)')
625 loops, best of 3: 5.67 µs per loop
sage: timeit('z=mpmath.gamma(100000)')
625 loops, best of 3: 50.9 µs per loop
</code></pre>
<p>Do you have a link to the definition of this function? I'm sure we can find a way to compute it more efficiently, especially given how many terms are currently noncontributing.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22019#post-id-22019Its Fisher's exact test. I was using the R implementation, but the smallest value it will give is 10^(-16), and I'd like more precision.Sat, 05 Mar 2011 03:55:07 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22019#post-id-22019Answer by fredrik for <p>I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.</p>
<p>Here is the function:</p>
<pre><code>from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
return num/denom
def myfisher_mp(a1,a2,b1,b2):
p = hyp_mp(a1,a2,b1,b2)
outp = p
for i in range(0,a1):
temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
if temp > p:
break
else:
outp += temp
for i in range(0,a2):
temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
if temp > p:
break
else:
outp += temp
return outp
</code></pre>
<p>A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12158#post-id-12158Use a recurrence instead of computing each term from scratch. The quotient between successive values of hyp_mp should be a rational function of i.
To avoid manual labor, I would usually use Mathematica to rewrite a hypergeometric quotient to rational form. I'm not sure if Sage currently has direct support for this.Fri, 04 Mar 2011 13:17:16 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12158#post-id-12158Comment by mhampton for <p>Use a recurrence instead of computing each term from scratch. The quotient between successive values of hyp_mp should be a rational function of i.</p>
<p>To avoid manual labor, I would usually use Mathematica to rewrite a hypergeometric quotient to rational form. I'm not sure if Sage currently has direct support for this.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22018#post-id-22018Thanks. Yes, using a recurrence would be much better, but I'm trying to get away with not thinking too much. If I had more time for this project it would be cool to implement Fisher's exact test for the general case (larger contingency tables, not just 2 by 2) but that is much harder.Sat, 05 Mar 2011 03:56:24 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22018#post-id-22018Answer by fredrik for <p>I think I could figure this out eventually, but I'm hoping it will be a very easy question for someone out there. I would like to convert the following function to cython to be as fast as possible. I am not sure exactly what I need to import from mpmath to do that.</p>
<p>Here is the function:</p>
<pre><code>from mpmath import *
mp.dps = 25; mp.pretty = True
def hyp_mp(a1,a2,b1,b2):
num = gamma(a1+a2+1)*gamma(b1+b2+1)*gamma(a1+b1+1)*gamma(a2+b2+1)
denom = gamma(a1+1)*gamma(a2+1)*gamma(b1+1)*gamma(b2+1)*gamma(a1+a2+b1+b2+1)
return num/denom
def myfisher_mp(a1,a2,b1,b2):
p = hyp_mp(a1,a2,b1,b2)
outp = p
for i in range(0,a1):
temp = hyp_mp(i,a2+a1-i,b1-i+a1,b2-a1+i)
if temp > p:
break
else:
outp += temp
for i in range(0,a2):
temp = hyp_mp(a1+a2-i,i,b1-a2+i,b2+a2-i)
if temp > p:
break
else:
outp += temp
return outp
</code></pre>
<p>A good test case would be myfisher_mp(1286, 9548, 133437, 148905), which takes about 1 second on my desktop. The inputs a1,a2,b1,b2 can be assumed to be positive ints.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12174#post-id-12174Actually, there is a simple method to speed this up quite significantly. Since the only needed gamma values are at integers (up to a few several thousand), you can just store all of them in a list or a dict.
In fact mpmath provides a decorator "memoize" that effectively does this. So you could just add
gamma = memoize(gamma)
to the program. With this change, the test case myfisher_mp(1286, 9548, 133437, 148905) runs about 2x faster the first time, and 7x faster the second time. But it should be even faster if you precompute the needed gamma values and replace the function calls with list lookups. Note that to precompute a list up to say n = 100000, you could just use a loop instead of repeatedly calling gamma, and this would be very fast.Sat, 05 Mar 2011 12:02:51 -0600http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?answer=12174#post-id-12174Comment by mhampton for <p>Actually, there is a simple method to speed this up quite significantly. Since the only needed gamma values are at integers (up to a few several thousand), you can just store all of them in a list or a dict.</p>
<p>In fact mpmath provides a decorator "memoize" that effectively does this. So you could just add</p>
<p>gamma = memoize(gamma)</p>
<p>to the program. With this change, the test case myfisher_mp(1286, 9548, 133437, 148905) runs about 2x faster the first time, and 7x faster the second time. But it should be even faster if you precompute the needed gamma values and replace the function calls with list lookups. Note that to precompute a list up to say n = 100000, you could just use a loop instead of repeatedly calling gamma, and this would be very fast.</p>
http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22002#post-id-22002Thanks, that's very helpful.Sun, 13 Mar 2011 18:15:47 -0500http://ask.sagemath.org/question/7983/convert-mpmath-function-to-cython/?comment=22002#post-id-22002