Need insight on code design choices. Context: scipy.integrate.odeint, Numba, and Python classes

[bhmer...@gmail.com]

Hi everyone,

scipy.integrate.odeint (http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.integrate.odeint.html) requires as its first argument, a function that computes the derivatives of the variables we want to integrate over (which I'll refer to as d_func, for "derivative function" from now on). d_func has to be written by the user, in Python code. A great way to get a boost of performance using Numba is to @jit the d_func (because d_func  is called *many* times during integration).

I have a question about how this optimization will work when d_func is complicated enough that it needs a Python class object behind it. 

Here is a "cartoon" of my code:

• there is a module called DynamicBox.py
• inside this module is a Python class, DynamicBox
• DynamicBox has a bunch of attributes
• some of these attributes are "phase variables" -- that is, I am interested in how they evolve in time
• some of these attributes are "parameters" -- that is, I use them to calculate the derivatives of the phase variables

Now, I will have a bunch of functions that will take DynamixBox phase variable or parameter attributes, in order to calculate relevant terms in the derivatives. That is:

• I will have a d_func 
• d_func itself will call lots of little helper functions to calculate relevant terms in the derivative, using DynamixBox attributes
  

So, I have to make a choice, with the following options:

1. either I can make d_func and all its helper functions methods of DynamicBox;
2. or I can make only d_func a method of DynamicBox, and all of its helper functions are in the same module as DynamicBox, but not methods of DynamicBox;
3. or only the helper functions are methods of DynamicBox, but d_func is just in the same module (DynamicBox.py), and not a method of DynamicBox;
4. or neither the helper functions, nor d_func, are methods of DynamicBox.

I do not know enough about Python to figure out which choice is best. Is it expensive to make instance attribute calls to get attributes? Is it expensive only if you are in a function that is not a method of the class? Questions like that is what I need insight on. What if Numba is in play? For instance, will Numba like it better if I am @jit-ting normal functions instead of class methods?

Can you help?

Kind regards,

Brian

[Jason Sachs]

I don't know the answer to this one, but it makes a REALLY interesting
use case that would be a great writeup of how to use Numba effectively
for a nontrivial task.

> --
> You received this message because you are subscribed to the Google Groups
> "Numba Public Discussion - Public" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to numba-users...@continuum.io.
> To post to this group, send email to numba...@continuum.io.
> To view this discussion on the web visit
> https://groups.google.com/a/continuum.io/d/msgid/numba-users/0aa079b6-b420-4208-a918-39b086c379c7%40continuum.io.
> For more options, visit https://groups.google.com/a/continuum.io/d/optout.

[Uwe Fechner]

I built a real-time kite-power flight simulator, using numba and
the RADAU5 solver from the assimulo solver suite.

The method, that computes the derivatives is called "residual".

You can look at the code here:
https://bitbucket.org/ufechner/freekitesim/src/5cc1adec08c4f9b1c7ef26813258fd1048cc6c92/KPS3_v2.py?at=master

The inner loops (calcAeroForces and loop) are functions (not methods), implemented with Numba.

This works very well, about 35 times faster then pure python.

Assimulo is better documented and supports more different solvers then scipy.integrate itself.
Have a look at:
http://www.jmodelica.org/assimulo

Best regards:

Uwe

[Frank Hellmann]

What I did in a similar, if not quite as complicated scenario (if I understand you correctly) was to use a function factory:

class ODE_Parameters(object)
    def __init__(self)
        self.x = np.ones(2)

@jit
def helper2(x)
    return(x)

def rhs_generator(ode_parameters, switch):

    x_arr = np.array(ode_paramseters.x, dtype = np.float64) # make sure we're doing double precision

    @jit
    def helper()
        return x_arr # This is a compile time constant

    @jit
    def rhs(t, y)
        return helper()

    @jit
    def rhs2(t, y)
        return helper2(x_arr)

    if switch   
        return rhs
    else
        return rhs2

This is a good approach if the parameters in your box don't change too often (or even depending on the phase evolution), as you would need to generate a new d_func every time.

[Brian Merchant]

Hi everyone. I also set up the same question on StackOverflow, along with a huge bounty (+350), so you might consider writing a nice answer there (and here -- in any case, I'll post any answer I get that's really cool over here). 

Some questions for the replies posted so far:

Uwe: did you end up testing the case where the inner loops were methods, not functions?

Frank: I am still not sure I understand why we need a function factory -- is this in order to keep Numba from constantly generating the RHS function again and again? Doesn't Numba generate only once at compile time?

[Uwe Fechner]

As far as I know, numba 0.18 cannot compile methods at all. I tried to use functions, where I pass
a reference to the object as first parameter to numba compiled functions, but this was very slow.

Currently only functions with scalars or numpy-arrays as parameters work well with numba.

There is a pull request to add simple class support, but it is not ready yet:
https://github.com/numba/numba/pull/1023