Interfacing GSL with Python using Cython / comparison with weave

The second example I wrote while testing how easy it is to use a C library from Python using Cython is the Gnu Scientific Library, gsl. The goal here is to compare how easy it is to do it using Cython or weave.

The first example is a bessel function (from gsl/gsl_sf_bessel.h):

double gsl_sf_bessel_J0(const double x);

The Cython interface is (gsl_test.pyx):

cdef extern from "gsl/gsl_sf_bessel.h":
    double gsl_sf_bessel_J0(double x)
        
def gsl_sf_bessel_jo(double x):
     return gsl_sf_bessel_J0 (x)     

The weave interface is defined in a pure python file. There is no need to have a full compilation round (setup.py ,etc) :

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

import numpy
import sys

if sys.platform == "linux2" :
    include_gsl_dir = "/usr/local/include/"
    lib_gsl_dir = "/usr/local/lib/"
elif sys.platform == "win32":    
    include_gsl_dir = r"c:\msys\1.0\local\include"
    lib_gsl_dir = r"c:\msys\1.0\local\lib"    
    
ext = Extension("gsl", ["gsl.pyx"],
    include_dirs=[numpy.get_include(), 
                  include_gsl_dir],
    library_dirs=[lib_gsl_dir],
    libraries=["gsl"]
)

setup(ext_modules=[ext],
    cmdclass = {'build_ext': build_ext})

The next step is to build you Cython module :

python setup.py build_ext --inplace

There is nothing to do with the weave extension that is built on the fly.

Finally all what you need is a function to compare the results from the two methods :

"""
Test the GSL cython interface
"""

from nose.tools import assert_almost_equal

import numpy

import gsl_test
import weave_test

def test_cython_gsl_bessel():    
    x = 5.0
    assert_almost_equal(gsl_test.gsl_sf_bessel_jo(x), -1.775967713143382920e-01, 15)
    
def test_weave_gsl_bessel():
    x = 5.0
    assert_almost_equal(weave_test.weave_gsl_bessel(x), -1.775967713143382920e-01, 15)

The two tests passes ok on my Windows system.

Going to the next step with more complex calls where we use numpy arrays with some stats function :

# importing the needed function definition
cdef extern from "gsl/gsl_sort.h":
    void gsl_sort (double * data, size_t stride, size_t n)
    
cdef extern from "gsl/gsl_statistics_double.h":
    double gsl_stats_median_from_sorted_data (double sorted_data[], size_t stride, size_t n)
    double gsl_stats_quantile_from_sorted_data (double sorted_data[], size_t stride, size_t n, double f)
      
def gsl_stat_demo(np.ndarray[np.float_t, ndim=1] inreal):
       """
       Computing some statistics using GSL on a one dimensional numpy array
       """
       cdef double median, upperq, lowerq   
     
       gsl_sort (<double *>inreal.data, 1, 5)
     
       median = gsl_stats_median_from_sorted_data (<double *>inreal.data, 1, 5)     
       upperq = gsl_stats_quantile_from_sorted_data (<double *>inreal.data, 1, 5, 0.75);
       lowerq = gsl_stats_quantile_from_sorted_data (<double *>inreal.data, 1, 5, 0.25);
     
       return (median, upperq, lowerq)

The same code in using weave is :

def weave_gsl_stat(ary):
       
    gsl_stat_code = """
    double median, upperq, lowerq;
    py::tuple results(3);
    
    gsl_sort (ary, 1, 5);
             
    median = gsl_stats_median_from_sorted_data (ary, 1, 5);     
    upperq = gsl_stats_quantile_from_sorted_data (ary, 1, 5, 0.75);
    lowerq = gsl_stats_quantile_from_sorted_data (ary, 1, 5, 0.25);
    
    results[0] = median;
    results[1] = upperq;
    results[2] = lowerq;
    return_val = results;
   """
    
    return inline(gsl_stat_code, ['ary'],           
           headers=['<gsl/gsl_statistics.h>', '<gsl/gsl_sort.h>'],
           libraries=['gsl'],
           include_dirs=[include_gsl_dir],
           library_dirs=[lib_gsl_dir],
           compiler='gcc')      

Writing two tests to check that we have the right results will confirm that the two version run ok :

def test_cython_gsl_stats():
    a = numpy.array([17.2, 18.1, 16.5, 18.3, 12.6])
    (median, upperq, lowerq) = gsl.gsl_stat_demo(a)
    assert_almost_equal(median, 17.2, 1)
    assert_almost_equal(upperq, 18.1, 1)
    assert_almost_equal(lowerq, 16.5, 1)
        
def test_weave_gsl_stats():
    a = numpy.array([17.2, 18.1, 16.5, 18.3, 12.6])
    (median, upperq, lowerq) = weave_test.weave_gsl_stat(a)
    assert_almost_equal(median, 17.2, 1)
    assert_almost_equal(upperq, 18.1, 1)
    assert_almost_equal(lowerq, 16.5, 1)

In conclusion, even if interfacing gsl with Python using Cython was quiet easy, it seems that using weave removes one layer of “complexity” (no setup.py file). This is convenient but requires that the environment where you run the file has a compiler, the needed headers and ibraries.

Advertisement