Saturday, December 10, 2011

cython odeint example

In response to the request by Ryan I now post a full example of using cython and odeint together. This is a specially written example because the original work contains too much code from my work that I do not want to post on the web. The equations used in the example have been chosen rather randomly.

In this example I have created a sample Model class which is a callable class (you can use an instance of the class as the function parameter for odeint). The Model class is defined in model_cython.pyx and contains all the Cython code. In ode_script.py the model function is integrated with odeint and thus shows how to call the Cython class from a normal python script. Setup.py is needed to compile model_cython.pyx (see this post).

model_cython.pyx:
import sys

#define some constants a c floats but still callable from python scripts
#not really used here but left in here as an example
cpdef float R_GAS,K_B
R_GAS=8.31      #Gas constant J/molK
K_B=1.38066e-23 #Boltzmann constant J/K

#import some often used math functions from the c math library
cdef extern from "math.h":
    double sin(double)

cdef extern from "math.h":
    double fabs(double)

cdef extern from "math.h":
    double exp(double)

#define a function callable from c and python
cpdef float Velocity(float t,float amplitude=1.0,float period=4.):
    cdef f1
    f1=fabs(sin(t*period))*amplitude
    return(f1)

class Model:
    """
    This class can be used to integrate a differential equation
    """
    def __init__(self,float a,float p):
        """
        -a is the amplitude of the sine function that describes the 
            speed oscillation as a function time
        -p controls the period of the sine function
        """
        self.a=a
        self.p=p
        
    
    def __call__(self,float y,float t):
        """
        y is the current state of the initial value problem:
        y: position
        t: time
        """
        cdef float pos
        pos=y
        #Workhardening
        v=Velocity(t,amplitude=self.a,period=self.p)*exp(-pos)
        return(v)

ode_script.py:
#!/bin/env/python

from numpy import r_
import pylab,scipy.integrate,sys
sys.path.append(".")
from model_cython import *


#ampltiude and period
amp=1.0
period=4.0
#set initial position
pos0=0.0
#create instance of Model object with speed amplitude of 1
myModel=Model(amp,period)
#create time vector (4 seconds in 500 steps)
tvec=r_[0.0:4.0:500j]
#integrate ov
y=(scipy.integrate.odeint(myModel,pos0,tvec))
    

#plot results
pylab.figure(1)
pylab.clf()
pylab.subplot(211)
pylab.plot(tvec,y)
pylab.xlabel("Time (s)")
pylab.ylabel("Position (m)")
pylab.grid(b=True)
pylab.subplot(212)
v_list=[]
for t in tvec:
    v_list.append(Velocity(t,amplitude=amp,period=period))
pylab.plot(tvec,v_list)
pylab.xlabel("Time (s)")
pylab.ylabel("Undamped velocity (m/s)")

pylab.grid(b=True)
pylab.show()

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

ext_modules = [Extension("model_cython",["model_cython.pyx"],\
    libraries=["m"])]

setup(
    name= 'Generic model class',
    cmdclass = {'build_ext': build_ext},
    ext_modules = ext_modules
)

Hope this helps.

Sunday, November 13, 2011

Ubuntu 11 samba mounting problems during booting

After upgrading to Ubuntu 11 I started to receive errors during booting mentioning CIFS-VFS. A quick search on the internet did not reveal a simple solution. Now I had a more thorough look. The CIFS-VFS error apparently has something to do with the automatic mounting of a samba share on a server in my network. I now know why the error pops up twice. From posts on the web I found that the most likely reason that I get this error is because the network is not up before the mounting is tried. Which makes sense when you get the error

mount error(101): Network is unreachable

(look in /var/log/boot.log)


Ubuntu runs the boot scripts in /etc/init. Here I found mountall.conf and mountall-net.conf. Here I found that the order of the scripts is determined by a line in the scripts in the form

start on runlevel [2345] stop on runlevel [!2345]

When I looked in /etc/init/mountall-net I found the line:

start on net-device-up

which suggests that the script waits for the network to be up before the script is run. The other script mountall.conf has the line "start on startup" in it which I think means that this script is called more or less at the beginning when the network indeed may not yet be up. When the network is up (or supposed to be up) mountall-net is called. This explains why I get the mount error twice during booting. While searching the web for "start on net-device-up" I somewhere found the variant


start on (net-device-up IFACE=eth0):

And this appears to make the problem go away (so now I only have one mnt error left during booting (from mountall.conf) but my samba share does get mounted automatically during booting).

Friday, October 14, 2011

working with openmp and gcc vectorization output

When working with openmp for parallelisation you can set the number of threads in a bash shell (cygwin) by

export OMP_NUM_THREADS=4

this can be good to do if you do not want to use all "cores" detected automatically on a cpu with hyperthreading.

To see which loops have been automatically vectorized by gcc you have to add the compiler switch

-ftree-vectorizer-verbose=5

Other verbose levels can also be used but for me level 5 seems to give the most relevant information.

Sunday, September 18, 2011

march=native

On the web I found an easy way to find what the gcc compiler option -march=native translates into (can't remember where I found this, it was somewhere on a gentoo forum):

echo "" | gcc -march=native -v -E - 2>&1 | grep cc1

On my Q6600 with gcc 4.4.3 this results in:

/usr/lib/gcc/i486-linux-gnu/4.4.3/cc1 -E -quiet -v - -D_FORTIFY_SOURCE=2 -march=core2 -mcx16 -msahf --param l1-cache-size=32 --param l1-cache-line-size=64 --param l2-cache-size=4096 -mtune=core2 -fstack-protector

Surprisingly on my laptop with an i7-2630QM and (cygwin) gcc 4.5.3 I get:

/usr/lib/gcc/i686-pc-cygwin/4.5.3/cc1.exe -E -quiet -v -D__CYGWIN32__ -D__CYGWIN__ -Dunix -D__unix__ -D__unix -idirafter /usr/lib /gcc/i686-pc-cygwin/4.5.3/../../../../include/w32api -idirafter /usr/lib/gcc/i686-pc-cygwin/4.5.3/../../../../i686-pc-cygwin/lib/. ./../include/w32api - -march=core2 -mcx16 -msahf -mpclmul -mpopcnt -mavx --param l1-cache-size=32 --param l1-cache-line-size=64 -- param l2-cache-size=6144 -mtune=generic

I was expecting a -march=corei7-avx but I guess that is only from gcc 4.6?

Sunday, July 24, 2011

Using Cython

The easiest way to convert old python code for use with cython is to rename your python file into something with the extension .pyx. This .pyx file can then be compiled with the setup.py (for distutils)script.

Here is an example setup.py which also shows how to include the math libray (-lm for gcc)

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

ext_modules = [Extension("BvL_cython",["BvL_cython.pyx"],\
    libraries=["m"])]

setup(
    name= 'BvL model class',
    cmdclass = {'build_ext': build_ext},
    ext_modules = ext_modules
)

The reason that I was interested in Cython was the long calculation times I encountered while doing a multi-variable optimization with a function evaluation that involved solving a differential equation with scipy.integrate.odeint. By simply replacing the class that contained the differential equation with a Cython version the calculation time dropped by a factor 5. Not bad for half a Sunday afternoons work.

First steps with Cython and Python(x,y)

Started today experimenting with Cython. With Cython you can write normal python code and compile it to C code. When you start to insert type information the speed up you get compared to normal python can be quite dramatic. The compiled c-code can be called from python (seemingly as if it was normal python code your calling). It is all nicely introduced here: Cython tutorial

On windows I am using Python(x,y). When you install it you have to tell the installer to include Cython. As c-compiler the Python (x,y) package includes mingw. Unfortunately Cython by default seems to expect a Visual C compiler. So the only way (so far) I found to use Cython is to run python from a normal command prompt when using the setup.py method (distutils) as described in the Cython Tutorial as follows:

python setup.py build_ext --inplace --compiler=mingw32

otherwise you get this error message:

error: Unable to find vcvarsall.bat