4

Quite often in numerical methods one has a lot coefficients which are static as they are fixed for the specific method. I was wondering what's the best way in Cython / C to set such arrays or variables.

In my case Runge-Kutta Integration methods are mostly the same, except for coefficients and the number of stages. Right now I'm doing something like (simplified)

# Define some struct such that it can be used for all different Runge-Kutta methods
ctypedef struct RKHelper:
    int numStages
    double* coeffs

cdef:
    RKHelper firstRKMethod
    # Later secondRKMethod, thirdRKMethod, etc.

firstRKMethod.numStages = 3
firstRKMethod.coeffs = <double*> malloc(firstRKMethod.numStages*sizeof(double))

# Arrays can be large and most entries are zero
for ii in range(firstRKMethod.numStages):
    firstRKMethod.coeffs[ii] = 0.

# Set non-zero elements
firstRKMethod.coeffs[2] = 1.3

Some points:

  • I know that malloc isn't for static arrays, but I don't know how to declare "numStages" or "RKHelper" as static in Cython, so I can't use a static array... Or I do something like "double[4]" in RKHelper, which doesn't allow to use the same struct definition for all RK methods.
  • I'm wondering if there is a better way than to do a loop. I don't wanna set the whole array manually (e.g. array = [0., 0., 0., 0., 1.3, ... lots of numbers mostly zero]).
  • As far as I can see there are no "real" static variables in Cython, are there?

Is there a nicer way of doing what I want to do?

Cheers

oli
  • 659
  • 1
  • 6
  • 18

2 Answers2

1

One way to achieve what you want to achieve is to set the coefficients of Runge Kutta scheme as global variables, that way you can use static arrays. This would be fast but definitely ugly

The ugly solution:

cdef int numStages = 3
# Using the pointer notation you can set a static array
# as well as its elements in one go
cdef double* coeffs = [0.,0.,1.3]
# You can always change the coefficients further as you wish

def RungeKutta_StaticArrayGlobal():
    # Do stuff

    # Just to check
    return numStages

A better solution would be to define a cython class with Runge Kutta coefficients as its members

The elegant solution:

cdef class RungeKutta_StaticArrayClass:
    cdef double* coeffs
    cdef int numStages
    def __cinit__(self):
        # Note that due to the static nature of self.coeffs, its elements
        # expire beyond the scope of this function    
        self.coeffs = [0.,0.,1.3]
        self.numStages = 3

    def GetnumStages(self):
        return self.numStages

    def Integrate(self):
        # Reset self.coeffs
        self.coeffs = [0.,0.,0.,0.,0.8,2.1]
        # Perform integration

Regarding your question of setting the elements, let's modify your own code with dynamically allocated arrays using calloc instead of malloc

The dynamically allocated version:

from libc.stdlib cimport calloc, free

ctypedef struct RKHelper:
    int numStages
    double* coeffs

def RungeKutta_DynamicArray():
    cdef:
        RKHelper firstRKMethod

    firstRKMethod.numStages = 3
    # Use calloc instead, it zero initialises the buffer, so you don't 
    # need to set the elements to zero within a loop
    firstRKMethod.coeffs = <double*> calloc(firstRKMethod.numStages,sizeof(double))

    # Set non-zero elements
    firstRKMethod.coeffs[2] = 1.3

    free(firstRKMethod.coeffs)

    # Just to check
    return firstRKMethod.numStages

Let's do a somewhat nonsensical benchmark, to verify that the arrays were truly static (i.e. had no runtime cost) in the first two examples

In[1]: print(RungeKutta_DynamicArray())
3
In[2]: print(RungeKutta_StaticArray())
3
In[3]: obj = RungeKutta_StaticArrayClass()
In[4]: print(obj.GetnumStages())
3

In[5]: %timeit RungeKutta_DynamicArray()
10000000 loops, best of 3: 65.2 ns per loop

In[6]: %timeit RungeKutta_StaticArray()
10000000 loops, best of 3: 25.2 ns per loop 

In[6]: %timeit RungeKutta_StaticArrayClass()
10000000 loops, best of 3: 49.6 ns per loop 

The RungeKutta_StaticArray essentially has close to a no-op cost, implying no runtime penalty for array allocation. You can choose to declare the coeffs within this function and the timing will still be the same. The RungeKutta_StaticArrayClass despite the overhead of setting a class with its members and constructors is still faster than the dynamically allocated version.

romeric
  • 2,325
  • 3
  • 19
  • 35
  • Thank you for the answer. Some things about your suggestions: StaticArray() is kinda annoying to write in the code if lot's of entries are identical (e.g. zero). I don't wanna have thousands of lines with zeros in my code. Basically I was hoping there was a way to tell Cython or the compiler to write a static array at compile time, i.e. "Cython: make a static array at compile time with this loop". A integrator class was my first implementation, but given the limitations of cdef classes with the GIL I'm trying to steer away from it for very basic blocks (like an integrator). – oli Feb 08 '16 at 09:36
0

Why not just use a numpy array? Practically it isn't actually static (see note at end), but you can allocate it at the global scope so it's created on module start-up. You can also access the raw C array underneath so there's no real efficiency cost.

import numpy as np

# at module global scope
cdef double[::1] rk_coeffs = np.zeros((50,)) # avoid having to manually fill with 0s
# illustratively fill the non-zero elements
rk_coeffs[1] = 2.0
rk_coeffs[3] = 5.0

# if you need to convert to a C array
cdef double* rk_coeffs_ptr = &rk_coeffs[0]

Note My reading of the question is that you're using "static" to mean "compiled into the module" rather than any of the numerous C-related definitions or anything to do with Python static methods/class variables.

Community
  • 1
  • 1
DavidW
  • 29,336
  • 6
  • 55
  • 86