Notebook

In [1]:

%load_ext cythonmagic
%load_ext memory_profiler
%load_ext autoreload


%autoreload 2
%pylab inline

Populating the interactive namespace from numpy and matplotlib

Pure Python¶

In [2]:

from math import sqrt
from numpy.random import randn
import numpy as np

def sdeevolvepy(x, tout, A, B, dt=.01):
    t = 0.0
    while t < tout:
        dt = min(dt, tout-t)

        ax0 = A(x)
        bx0 = B(x)
        
        x += dt * ax0 + sqrt(dt) * bx0 * randn()
        t += dt
        
    return x
    
def sdepy(x, tout, A, B, dt=.01):
    t = tout[0]
    
    xout = np.empty_like(tout)
    dt = min(dt, np.diff(tout).min())
    
    for i, (t0, tNext) in enumerate(zip(tout[:-1], tout[1:])):
        x = sdeevolvepy(x, tNext-t0, A, B, dt=dt)
        xout[i] = x
        
    return tout, xout
    

Cython¶

In [3]:

%%cython -lgsl
#cython: boundscheck=False
#cython: wraparound=False
##cython: profile=True
cimport cython
from cython.view cimport array as cvarray
cimport numpy as np
import numpy as np
from libc.math cimport sqrt, fmin



cdef extern from "gsl/gsl_rng.h":
    ctypedef struct gsl_rng_type:
        pass
    ctypedef struct gsl_rng:
        pass
    gsl_rng_type *gsl_rng_mt19937
    gsl_rng *gsl_rng_alloc(gsl_rng_type * T)
 
cdef extern from "gsl/gsl_randist.h":
    double gamma "gsl_ran_gamma"(gsl_rng * r,double,double)
    double gaussian "gsl_ran_gaussian"(gsl_rng * r,double)
 
cdef gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)





cdef class RandStream:
    cdef gsl_rng * _strm 

    def __cinit__(self):
        self._strm = gsl_rng_alloc(gsl_rng_mt19937)

    def __call__(self, double sig):
        return gaussian(self._strm, sig)
    
    cdef double normal(self, double sig):
        return gaussian(self._strm, sig)

    
cdef RandStream rv = RandStream()

cdef double A(double x):
    return -x
cdef double B(double x):
    return 1.0




cdef double sdeevolve(double x, double tout, double dt ):
    cdef double t, ax0, bx0
#     cdef np.ndarray[np.float64_t] rands = np.random.randn(int(tout/dt) + 1)
    cdef int i = 0
    cdef bint dobreak
    
    t = 0.0
    
    if tout < 1e-9 : t  = tout +1
    
    while t < tout:
            
        if t + dt > tout:
            dt = tout - t
            dobreak = True
        

        ax0 = A(x)
        bx0 = B(x)
        
#         x += dt * ax0 + sqrt(dt) * bx0 * rands[i]
        x += dt * ax0 + sqrt(dt) * bx0 * rv.normal(sqrt(dt))
        
        if dobreak:
            break

        t += dt
        i+=1
        
    return x


def sde(double x, np.ndarray[np.float64_t] tout, double dt=.01):
    
    cdef double t, tNext, t0

    
    t = tout[0]
    
    cdef np.ndarray[np.float64_t] xout = np.empty_like(tout)
    dt = min(dt, np.diff(tout).min())
    
    cdef int i

    xout[0] = x
    for i in xrange(tout.shape[0]-1):
        t0 = tout[i]
        tNext = tout[i+1]
        x= sdeevolve(x, tNext-t0, dt)
        xout[i+1] = x
        
    return tout, xout

In [4]:

x0 = 0.0
tout = np.mgrid[0:100:.1]
tout ,xout = sde(x0, tout, dt=.0001)
plot(tout, xout)

Out[4]:

[<matplotlib.lines.Line2D at 0x10d002a90>]

Using MKL¶

In [5]:

%%file sdemkl.pyx
#cython: boundscheck=False
#cython: wraparound=False
##cython: profile=True
cimport cython
from cython.view cimport array as cvarray
cimport numpy as np
import numpy as np
from libc.math cimport sqrt, fmin

ctypedef void * mklRandomStream
cdef extern from "mkl.h":
    int VSL_BRNG_MT19937 
    int VSL_RNG_METHOD_GAUSSIAN_BOXMULLER
    int mkstream "vslNewStream"(mklRandomStream *, const int  , const int)
    int mklgaussian "vdRngGaussian"(const int  ,mklRandomStream  , const int  , double [], const double  , const double  )

cdef  mklRandomStream mklstream


cdef double A(double x):
    return -x
cdef double B(double x):
    return 1.0

cdef  double a = 0



cdef class RandN:
    cdef mklRandomStream _strm 

    def __cinit__(self):
        mkstream(&self._strm, VSL_BRNG_MT19937, 111)

    def __call__(self, double sig):
        cdef double r
        mklgaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, self._strm, 1, &r, a, sig )
        return r

    cdef double n (self, double sig):
        cdef double r
        mklgaussian( VSL_RNG_METHOD_GAUSSIAN_BOXMULLER, self._strm, 1, &r, a, sig )
        return r



cdef RandN rv = RandN()

cdef double sdeevolve(double x, double tout, double dt):
    cdef double t, ax0, bx0
#     cdef np.ndarray[np.float64_t] rands = np.random.randn(int(tout/dt) + 1)
    cdef int i = 0
    cdef bint dobreak

    
    t = 0.0
    
    if tout < 1e-9 : t  = tout +1
    
    while t < tout:
            
        if t + dt > tout:
            dt = tout - t
            dobreak = True
        

        ax0 = A(x)
        bx0 = B(x)
        
#         x += dt * ax0 + sqrt(dt) * bx0 * rands[i]
        x += dt * ax0 +  bx0 * rv.n(sqrt(dt)) 
        
        if dobreak:
            break

        t += dt
        i+=1
        
    return x


def sde(double x, np.ndarray[np.float64_t] tout, double dt=.01):
    
    cdef double t, tNext, t0
    
    t = tout[0]
    
    cdef np.ndarray[np.float64_t] xout = np.empty_like(tout)
    dt = min(dt, np.diff(tout).min())
    
    cdef int i

    xout[0] = x
    for i in xrange(tout.shape[0]-1):
        t0 = tout[i]
        tNext = tout[i+1]
        x= sdeevolve(x, tNext-t0, dt)
        xout[i+1] = x
        
    return tout, xout

Overwriting sdemkl.pyx

In [6]:

%%file setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
import numpy as np                           # <---- New line

ext_modules = [Extension("sde",
    sources = ["sdemkl.pyx"],
    libraries = ['mkl_intel_lp64', 'mkl_sequential', 'mkl_core', 'pthread', 'm'],
    )]

setup(
    name = 'sde',
    cmdclass = {'build_ext': build_ext},
    include_dirs = [np.get_include(), '/opt/intel/mkl/include'],         # <---- New line
    library_dirs=['/opt/intel/mkl/lib'],
    ext_modules = ext_modules
)
# -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm

Overwriting setup.py

In [7]:

!python setup.py build_ext --inplace

/Applications/Canopy.app/appdata/canopy-1.4.1.1975.macosx-x86_64/Canopy.app/Contents/lib/python2.7/distutils/dist.py:267: UserWarning: Unknown distribution option: 'library_dirs'
  warnings.warn(msg)
running build_ext
cythoning sdemkl.pyx to sdemkl.c
building 'sde' extension
gcc -fno-strict-aliasing -fno-common -dynamic -arch x86_64 -DNDEBUG -g -O3 -arch x86_64 -I/Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include -I/opt/intel/mkl/include -I/Applications/Canopy.app/appdata/canopy-1.4.1.1975.macosx-x86_64/Canopy.app/Contents/include/python2.7 -c sdemkl.c -o build/temp.macosx-10.6-x86_64-2.7/sdemkl.o
In file included from sdemkl.c:314:
In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h:4:
In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/ndarrayobject.h:17:
In file included from /Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/ndarraytypes.h:1761:
/Users/noah/Library/Enthought/Canopy_64bit/User/lib/python2.7/site-packages/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: 
      "Using deprecated NumPy API, disable it by "          "#defining
      NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-W#warnings]
#warning "Using deprecated NumPy API, disable it by " \
 ^
1 warning generated.
gcc -bundle -undefined dynamic_lookup -g -arch x86_64 -headerpad_max_install_names -arch x86_64 build/temp.macosx-10.6-x86_64-2.7/sdemkl.o -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -o /Users/noah/Dropbox/gnl/ipynb/sde.so

In [8]:

import sde as sdemkl
tout, xout = sdemkl.sde(x0, tout, dt=.01)
plot(tout, xout)

Out[8]:

[<matplotlib.lines.Line2D at 0x10d02b290>]

Comparison¶

In [9]:

tout = np.mgrid[0:100.0:.1]
x0 = 0
dtmax= .001

%timeit sdemkl.sde(x0, tout, dt=dtmax)
%timeit sde(x0, tout, dt=dtmax)
%timeit sdepy(x0, tout, lambda x:-x, lambda x:1.0, dt=dtmax)

1000 loops, best of 3: 1.04 ms per loop
10000 loops, best of 3: 162 µs per loop
1 loops, best of 3: 428 ms per loop

GSL is far and away the winner. It is around 3000x faster than the pure python code.

Theano way to calculate SDE¶

In [10]:

import theano
import theano.tensor as T
from theano.tensor import sqrt

from theano.tensor.shared_randomstreams import RandomStreams

In [11]:

def A(x):
    return -x
def B(x):
    return 1.0

def onestep(x, dt):
    return x + A(x) * dt + sqrt(dt) * B(x) * rv

srng = RandomStreams(seed=31415)

x = T.dscalar("x")
dt = T.dscalar("dt")

rv = srng.normal()

xout, updates = theano.scan(onestep, non_sequences=[dt], n_steps=100, outputs_info=[x,])

sde = theano.function(inputs=[x, dt], outputs=xout, updates=updates, allow_input_downcast=True)

/Users/noah/epd/lib/python2.7/site-packages/theano/scan_module/scan_perform_ext.py:85: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
  from scan_perform.scan_perform import *

In [12]:

sde(0.0, .01)

Out[12]:

array([ 0.09437863,  0.24098799,  0.16236896,  0.33597058,  0.39227402,
        0.42859804,  0.4297515 ,  0.41065422,  0.42676154,  0.46194971,
        0.23951774,  0.22655675,  0.22447993,  0.2832613 ,  0.30663556,
        0.2521183 ,  0.15739451,  0.20163056,  0.44691499,  0.49963293,
        0.44309233,  0.45447732,  0.76954386,  0.99497569,  0.94931493,
        1.00052418,  0.94978512,  0.87785128,  0.94674234,  0.99830616,
        1.01836497,  1.13322821,  1.24859349,  1.10906753,  1.1890498 ,
        1.14122232,  1.09423422,  1.08703787,  1.01070417,  1.0176151 ,
        0.93632021,  0.97563049,  0.87088073,  0.697739  ,  0.73428342,
        0.71105801,  0.7295656 ,  0.77290335,  0.87487668,  0.70402478,
        0.84913249,  0.92544088,  0.86158857,  0.85057093,  0.9961133 ,
        1.07968195,  0.97135231,  0.86316631,  0.7943823 ,  0.72409596,
        0.92205431,  0.8814426 ,  1.0044605 ,  1.00126224,  0.90335298,
        0.82920857,  0.70217597,  0.7704939 ,  0.69542109,  0.56746137,
        0.60785471,  0.58045189,  0.6764967 ,  0.63270697,  0.66330246,
        0.72481229,  0.64581193,  0.7229394 ,  0.72425885,  0.82994694,
        0.78466216,  0.76627463,  0.83304857,  0.68643934,  0.58504504,
        0.42661802,  0.51153621,  0.4141255 ,  0.37583259,  0.35615542,
        0.24124443,  0.11700228,  0.0925223 ,  0.01597322, -0.04461747,
       -0.01722394,  0.0083655 ,  0.17626539,  0.09930939,  0.13939293])