Ask Your Question

Fastest way to call special function (elliptic integral) from cython code for gsl ode_solver()

asked 2014-06-27 12:55:47 +0200

argentum2f gravatar image

updated 2015-01-13 21:13:33 +0200

FrédéricC gravatar image

I am using sage.gsl.ode.ode_solver to solve an ode, and overloading the c_f() function as detailed here.

The code is working using an approximate c_f() function. For the exact one, I need to use the complete elliptic integrals. These are special functions that are not part of the libc.math library. They seem to be available through a couple packages on sage. I found I could call them by using sage.all.elliptic_ec(k) and sage.all.elliptic_kc(k)

This is causing the code to slow down by a factor of about 10^3. I know those special functions are going to be expensive, but timing the elliptic_ec() function in a separate cell returns 0 for the time, so these functions aren't that slow.

I'm wondering if the problem is mostly just that I'm calling a python function from cython and losing the speed up because of that? Is there a better way to do that? There are c libraries available with those functions - would it be better to import the files into my project, and call them from the cython code? (not sure how to do that...) Sage seems to have some gsl packages already - does it also have those special functions through gsl? Is there a different package with a faster form of those functions?

Disclaimer: I've played around with sage a bit in the past, but this is the first intensive numerical simulation I've attempted with it, so sage/python/cython are all relatively new to me. I'm more accustomed to Mathematica, Matlab/Octave and c++

Here are (I believe) the relevant bits of code:

from libc.math cimport pow
cimport sage.gsl.ode
import sage.gsl.ode
include 'gsl.pxi'

cdef class zeeman_acceleration(sage.gsl.ode.ode_system):
    ... other stuff ...
    cdef double coilFieldR(self, double rp, double zp, double coilR): # T
        cdef double B0, alpha, beta, Q, k
        ... some basic arithmetic ...
        if rp > coilR/10e4 :
            gamma = zp/rp
            EC = sage.all.elliptic_ec(k)
            EK = sage.all.elliptic_kc(k)
            Br = B0*gamma/(self.pi*sqrt(Q))*(EC*(1+alpha**2+beta**2)/(Q-4*alpha)-EK)
            return Br
        return 0
... other stuff (including c_f() function that calls coilField...
edit retag flag offensive close merge delete

1 Answer

Sort by » oldest newest most voted

answered 2014-06-27 17:01:33 +0200

kcrisman gravatar image

updated 2014-06-27 17:02:44 +0200

I'm wondering whether using mpmath directly (e.g. the elliptic integral of the first kind) might be more efficient for you. You should certainly be able to call this in Python, presumably also nicely in Cython.

Edit: mpmath is of course numeric. Our exact functionality is presumably provided by Maxima; you may be able to interact directly with the Lisp (ECL) interface, though. And it should only slow down the first time you evaluate, when Maxima is started... hmm.

edit flag offensive delete link more


I'm using - Is there any simple way of figure out what's available and how to get to it? I found that hitting tab auto completes, and after a lot of guesswork and navigating auto complete trees I did finally find the elliptic functions here: **sage.libs.mpmath.ext_main.global_context.ellipk** It has sped things up by about a factor of 5. Wouldn't it be faster if I cimport sage.mpmath or something (like I have to do with sage.gsl.ode)? I tried that and sage complains about their not being a pyx file or something. Do I have to write my own, and would it help at all?

argentum2f gravatar imageargentum2f ( 2014-06-27 21:36:35 +0200 )edit

`mpmath` isn't Cython, it's regular Python, as far as I know. You should be able to just use `from mpmath import ellipk` or something, but I am not a Cython expert so perhaps you do need that long of a tree.

kcrisman gravatar imagekcrisman ( 2014-06-30 14:58:03 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools


Asked: 2014-06-27 12:55:47 +0200

Seen: 6,170 times

Last updated: Jun 27 '14