Ask Your Question
2

Passing params to cythonized ode_system()

asked 2013-01-25 00:10:27 +0100

gitano gravatar image

Hi, I am trying to code a batch process to solve an ODE many times over a range of parameters. I am using a modified version of an example in the sage documentation:

%cython
cimport sage.gsl.ode
import sage.gsl.ode
include 'gsl.pxi'

cdef class duffing_osc(sage.gsl.ode.ode_system):
    cdef int c_f(self,double t, double *y,double *dydt):
        dydt[0]=y[1]
        dydt[1]=y[0]-params[0]*y[1]-y[0]*y[0]*y[0]+params[1]*cos(params[2]*t)
        return GSL_SUCCESS
    cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
        dfdy[0]=0
        dfdy[1]=1.0
        dfdy[2]=-3*y[0]*y[0]+1
        dfdy[3]=-params[0]*y[1]
        dfdt[0]=0
        dfdt[1]=-params[1]*params[2]*sin(params[2]*t)
        return GSL_SUCCESS

Now, the problem of course is that when trying to generate the .spyx file, param is undeclared, since I can't overload the c_f and c_j functions. Is there a way I can further modify my derived class (duffing_osc) so that when I instantiate an object I can pass in the parameters, such as d = duffing_osc(params). Then I suppose I would have to change params[i] to self.params[i] in my definitions for c_f and c_j.

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted
2

answered 2013-01-26 00:41:37 +0100

I guess you mean the example at the end of the Solving ODE numerically by GSL section.

You can add a constructor to your class which takes params as an argument and stores it in the instance created. The code to use this would be (mostly copied from the reference manual)

sage: params = [<value>, <value>, <value>]
sage: T = ode_solver()
sage: T.algorithm = "bsimp"
sage: osc = duffing_osc(params) # note the argument
sage: T.function = osc
sage: T.ode_solve(...)

The new class would look like:

cdef class duffing_osc(sage.gsl.ode.ode_system):
    def __init__(self, params):
        self.params = params

    cdef int c_f(self,double t, double *y,double *dydt):
        dydt[0]=y[1]
        dydt[1]=y[0]-self.params[0]*y[1]-y[0]*y[0]*y[0]+\
                self.params[1]*cos(self.params[2]*t)
        return GSL_SUCCESS
    cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
        dfdy[0]=0
        dfdy[1]=1.0
        dfdy[2]=-3*y[0]*y[0]+1
        dfdy[3]=-self.params[0]*y[1]
        dfdt[0]=0
        dfdt[1]=-self.params[1]*self.params[2]*sin(self.params[2]*t)
        return GSL_SUCCESS

In terms of performance, there is a problem with this code.

The elements of params are Sage types, either Python floats, RDF elements, or even multiprecision floats (backed by MPFR). Whenever these are multiplied by a double, Sage will automatically convert this value to a Sage/Python object, do the arithmetic, and convert the result back to a double. This is very slow, compared to a simple float multiplication.

To work around this, you can convert params to a vector of doubles in the constructor. Assuming params has 3 elements, then the first lines of the class would be

cdef class duffing_osc(sage.gsl.ode.ode_system):
    cdef double params[3]
    def __init__(self, params):
        for i in range(3):
            self.params[i] = params[i] # type conversion is automatic
edit flag offensive delete link more

Comments

Thanks a lot. That seems to get the job done quite nicely.

gitano gravatar imagegitano ( 2013-01-26 03:11:32 +0100 )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

Stats

Asked: 2013-01-25 00:10:27 +0100

Seen: 392 times

Last updated: Jan 26 '13