Commits

Anonymous committed 1cf0a4c

added samples/{logistic_map,gslctypes} for an example of railgun.cmem

Comments (0)

Files changed (10)

samples/gslctypes/gslctypes.rst

+.. _gslctypes:
+
+Using RailGun with GSL (GNU Scientific Library)
+===============================================
+
+Class :class:`gslctypes_rng.cmem` defined in the following file
+(``gslctypes_rng.py``) is a class you can give to :func:`railgun.cmem`.
+See :ref:`logistic_map` for usage.
+
+.. literalinclude:: gslctypes_rng.py
+   :language: python

samples/gslctypes/gslctypes_rng.py

+import ctypes
+from ctypes import (POINTER, c_char_p, c_size_t, c_int, c_long, c_ulong,
+                    c_double, c_void_p)
+from ctypes.util import find_library
+
+
+class _c_gsl_rng_type(ctypes.Structure):
+    _fields_ = [('name', c_char_p),
+                ('max', c_long),
+                ('min', c_size_t),
+                ('__set', c_void_p),
+                ('__get', c_void_p),
+                ('__get_double', c_void_p),
+                ]
+_c_gsl_rng_type_p = POINTER(_c_gsl_rng_type)
+
+
+class _c_gsl_rng(ctypes.Structure):
+    _fields_ = [('type', _c_gsl_rng_type_p),
+                ('state', c_void_p)]
+_c_gsl_rng_p = POINTER(_c_gsl_rng)
+
+
+class _GSLFuncLoader(object):
+
+    # see: http://code.activestate.com/recipes/576549-gsl-with-python3/
+    gslcblas = ctypes.CDLL(find_library('gslcblas'), mode=ctypes.RTLD_GLOBAL)
+    gsl = ctypes.CDLL(find_library('gsl'))
+
+    def _load_1(self, name, argtypes=None, restype=None):
+        func = getattr(self.gsl, name)
+        if argtypes is not None:
+            func.argtypes = argtypes
+        if restype is not None:
+            func.restype = restype
+        setattr(self, name, func)
+        return func
+
+    def _load(self, name, argtypes=None, restype=None):
+        if isinstance(name ,basestring):
+            return self._load_1(name, argtypes, restype)
+        else:
+            try:
+                return [self._load_1(n, argtypes, restype) for n in name]
+            except TypeError:
+                raise ValueError('name=%r should be a string or iterative '
+                                 'of string' % name)
+
+
+func = _GSLFuncLoader()
+func._load('gsl_strerror', [c_int], c_char_p)
+func._load('gsl_rng_alloc', [_c_gsl_rng_type_p], _c_gsl_rng_p)
+func._load('gsl_rng_set', [_c_gsl_rng_p, c_ulong])
+func._load('gsl_rng_free', [_c_gsl_rng_p])
+func._load('gsl_rng_types_setup',
+           restype=c_void_p)  # POINTER(_c_gsl_rng_p)
+func._load(['gsl_ran_gaussian',
+            'gsl_ran_gaussian_ziggurat',
+            'gsl_ran_gaussian_ratio_method'],
+           [_c_gsl_rng_p, c_double],
+           c_double)
+
+
+gsl_strerror = func.gsl_strerror
+
+
+def _get_gsl_rng_type_p_dict():
+    """
+    Get all ``gsl_rng_type`` as dict which has pointer to each object
+
+    This is equivalent to C code bellow which is from GSL document:
+
+    .. sourcecode:: c
+
+          const gsl_rng_type **t, **t0;
+          t0 = gsl_rng_types_setup ();
+          for (t = t0; *t != 0; t++)
+            {
+              printf ("%s\n", (*t)->name);  /* instead, store t to dict */
+            }
+
+    """
+    t = func.gsl_rng_types_setup()
+    dt = ctypes.sizeof(c_void_p)
+    dct = {}
+    while True:
+        a = c_void_p.from_address(t)
+        if a.value is None:
+            break
+        name = c_char_p.from_address(a.value).value
+        dct[name] = ctypes.cast(a, _c_gsl_rng_type_p)
+        t += dt
+    return dct
+
+
+class gsl_rng(object):
+    _gsl_rng_alloc = func.gsl_rng_alloc
+    _gsl_rng_set = func.gsl_rng_set
+    _gsl_rng_free = func.gsl_rng_free
+    _gsl_rng_type_p_dict = _get_gsl_rng_type_p_dict()
+    _ctype_ = _c_gsl_rng_p  # for railgun
+
+    def __init__(self, seed=None, name='mt19937'):
+        self._gsl_rng_name = name
+        self._gsl_rng_type_p = self._gsl_rng_type_p_dict[name]
+        self._cdata_ = self._gsl_rng_alloc(self._gsl_rng_type_p)
+        # the name '_cdata_' is for railgun
+        if seed is not None:
+            self.set(seed)
+
+    def __del__(self):
+        self._gsl_rng_free(self._cdata_)
+
+    def set(self, seed):
+        self._gsl_rng_set(self._cdata_, seed)
+
+    _gsl_ran_gaussian = {
+        '': func.gsl_ran_gaussian,
+        'ziggurat': func.gsl_ran_gaussian_ziggurat,
+        'ratio_method': func.gsl_ran_gaussian_ratio_method,
+        }
+
+    def ran_gaussian(self, sigma=1.0, method=''):
+        return self._gsl_ran_gaussian[method](self._cdata_, sigma)

samples/gslctypes/plot_gaussian.py

+import pylab
+from gslctypes_rng import gsl_rng
+
+method = 'ziggurat'
+sigma = 1.0
+rng = gsl_rng()
+pylab.hist(
+    [rng.ran_gaussian(method=method, sigma=sigma) for i in xrange(1000)],
+    bins=20, normed=True)
+pylab.show()

samples/index.rst

 
    lode/*
    lode_rk4/*
+   logistic_map/*
+   glsctypes/*
    */*

samples/logistic_map/Makefile

+liblogistic_map.so: logistic_map.c
+	$(CC) $(CFLAGS) $(CPPFLAGS) -lgsl -lgslcblas -Wall -Wextra -shared $< -o $@
+
+clean:
+	-rm *.so

samples/logistic_map/bifurcation_diagram.py

+from logistic_map import LogisticMap
+import pylab
+
+lmap = LogisticMap(1000)
+
+pylab.figure(1)
+lmap.plot_bifurcation_diagram(0, 4, 500, sigma=1e-5)
+pylab.title(r'$\sigma=10^{-5}$')
+
+pylab.figure(2)
+lmap.plot_bifurcation_diagram(0, 4, 500, sigma=1e-3)
+pylab.title(r'$\sigma=10^{-3}$')
+
+pylab.figure(3)
+lmap.plot_bifurcation_diagram(0, 4, 500, sigma=1e-2)
+pylab.title(r'$\sigma=10^{-2}$')
+
+pylab.show()

samples/logistic_map/gslctypes_rng.py

+../gslctypes/gslctypes_rng.py

samples/logistic_map/logistic_map.c

+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+typedef struct logisticmap_{
+  int num_i;
+  double *xt;
+  double mu;
+  double sigma;
+  gsl_rng *rng;
+} LogisticMap;
+
+
+int LogisticMap_gene_seq(LogisticMap *self)
+{
+  int i;
+  double eta;
+  for (i = 1; i < self->num_i; ++i){
+    eta = gsl_ran_gaussian_ziggurat(self->rng, self->sigma);
+    self->xt[i] = self->mu * self->xt[i-1] * (1 - self->xt[i-1]) + eta;
+    if (self->xt[i] < 0){
+      self->xt[i] = 0;
+    }
+  }
+  return 0;
+}
+

samples/logistic_map/logistic_map.py

+from railgun import SimObject, relpath, cmem
+import numpy
+from gslctypes_rng import gsl_rng
+
+class LogisticMap(SimObject):
+
+    _clibname_ = 'liblogistic_map.so'
+    _clibdir_ = relpath('.', __file__)
+    _cmembers_ = [
+        'num_i',
+        'double xt[i]',
+        'double mu',
+        'double sigma',
+        cmem(gsl_rng, 'rng'),
+        ]
+    _cfuncs_ = ["xt gene_seq()"]
+
+    def __init__(self, num_i, seed=None, **kwds):
+        SimObject.__init__(self, num_i=num_i, rng=gsl_rng(seed), **kwds)
+
+    def _cwrap_gene_seq(old_gene_seq):
+        def gene_seq(self, **kwds):
+            # use the previous "last state" as initial state, but setv
+            # may overwrite this by `xt_0=SOMEVAL`
+            self.xt[0] = self.xt[-1]
+            self.setv(**kwds)
+            return old_gene_seq(self)
+        return gene_seq
+
+    def bifurcation_diagram(self, mu_start, mu_stop, mu_num, **kwds):
+        mu_list = numpy.linspace(mu_start, mu_stop, mu_num)
+        self.gene_seq(mu=mu_list[0], **kwds)  # through first steps away
+        xt_list = [self.gene_seq(mu=mu).copy() for mu in mu_list]
+        return (mu_list, xt_list)
+
+    def plot_bifurcation_diagram(self, mu_start, mu_stop, mu_num, **kwds):
+        import pylab
+        bd = self.bifurcation_diagram(mu_start, mu_stop, mu_num, **kwds)
+        ones = numpy.ones(self.num_i)
+
+        for (mu, x) in zip(*bd):
+            pylab.plot(ones * mu, x, 'ko', markersize=0.5)
+        pylab.ylim(0, 1)

samples/logistic_map/logistic_map.rst

+.. _logistic_map:
+
+Logistic map with noise additive (usage of :func:`railgun.cmem`)
+================================================================
+
+You will need the file ``gslctypes_rng.py`` (see :ref:`gslctypes`).
+
+Python code
+-----------
+``logistic_map.py``:
+
+.. literalinclude:: logistic_map.py
+   :language: python
+
+``bifurcation_diagram.py``:
+
+.. literalinclude:: bifurcation_diagram.py
+   :language: python
+
+C code
+------
+``logistic_map.c``:
+
+.. literalinclude:: logistic_map.c
+   :language: c
+
+Output
+------
+.. plot:: source/samples/logistic_map/bifurcation_diagram.py