python/sddekit.pyx

   1 # copyright 2016 Apache 2 sddekit authors
   2
   3 from libc.stdlib cimport malloc, free
   4 from libc.stdint cimport uint32_t
   5 from cpython cimport PyObject, Py_INCREF, Py_DECREF
   6 cimport numpy as np
   7 import numpy as np
   8 np.import_array()
   9
  10 # gc-friendly c array wrapper
  11 # based on https://gist.github.com/GaelVaroquaux/1249305 BSD-licensed.
  12 cdef class CArrayWrapper:
  13     "Python wrapper for C array."
  14
  15     cdef void* ptr
  16     cdef int size
  17     cdef int typenum
  18     cdef int free_on_gc
  19
  20     cdef init(self, int size, void *ptr, int typenum, int free_on_gc=0):
  21         "Initialize wrapper."
  22         self.ptr = ptr
  23         self.size = size
  24         self.typenum = typenum
  25         self.free_on_gc = free_on_gc
  26
  27     def __array__(self):
  28         "Obtain NumPy array from self."
  29         cdef np.npy_intp shape[1]
  30         shape[0] = <np.npy_intp> self.size
  31         return np.PyArray_SimpleNewFromData(1, shape, self.typenum, self.ptr)
  32
  33     def __dealloc__(self):
  34         "Deallocate C array on Python garbage collect."
  35         if self.free_on_gc:
  36             free(<void*> self.ptr)
  37
  38 cdef np.ndarray np_of_ptr(int size, void *ptr, int typenum, int free_on_gc):
  39     "Wrap a given C array, returning a NumPy array."
  40     cdef np.ndarray ndarray
  41     wrapper = CArrayWrapper()
  42     wrapper.init(size, ptr, typenum, free_on_gc=free_on_gc)
  43     ndarray = np.array(wrapper, copy=False)
  44     ndarray.base = <PyObject*> wrapper
  45     Py_INCREF(wrapper)
  46     return ndarray
  47
  48 # hist
  49 cdef extern from "sddekit.h":
  50
  51     ctypedef enum sd_stat:
  52         SD_OK
  53         SD_ERR
  54         SD_CONT
  55         SD_STOP
  56
  57     ctypedef struct sd_hfill:
  58         void *ptr
  59         sd_stat (*apply)(sd_hfill*, uint32_t n, double * t, uint32_t *indices, double * buf)
  60         void (*free)(sd_hfill*)
  61
  62     ctypedef struct sd_hist:
  63         uint32_t(*get_maxvi)(sd_hist*)
  64         uint32_t(*get_vi2i)(sd_hist*, uint32_t vi)
  65         uint32_t(*get_nu)(sd_hist *h)
  66         void (*free)(sd_hist *h)
  67         sd_stat (*fill)(sd_hist *h, sd_hfill *filler)
  68         void (*get)(sd_hist *h, double t, double *aff)
  69         void (*set)(sd_hist *h, double t, double *eff)
  70         uint32_t (*nbytes)(sd_hist *h)
  71         double (*get_buf_lin)(sd_hist *h, uint32_t index)
  72         uint32_t (*get_nd)(sd_hist *h)
  73         double (*get_t)(sd_hist *h)
  74         double (*get_dt)(sd_hist *h)
  75         uint32_t (*get_lim)(sd_hist *h, uint32_t i)
  76         uint32_t (*get_len)(sd_hist *h, uint32_t i)
  77         uint32_t (*get_pos)(sd_hist *h, uint32_t i)
  78         uint32_t (*get_uvi)(sd_hist *h, uint32_t i)
  79         double (*get_maxd)(sd_hist *h, uint32_t i)
  80         uint32_t (*get_vi)(sd_hist *h, uint32_t i)
  81         double (*get_vd)(sd_hist *h, uint32_t i)
  82
  83     sd_hist * sd_hist_new_default(uint32_t nd, uint32_t *vi, double *vd, double t0, double dt)
  84
  85
  86 # hfill
  87 ctypedef struct hfill_py:
  88     sd_hfill hf
  89     void *fn
  90
  91 cdef sd_stat hfill_py_apply(sd_hfill *hfill, uint32_t n, double * t, uint32_t *i, double * buf):
  92     cdef hfill_py *hp = <hfill_py*> hfill.ptr
  93     cdef np.ndarray[np.float64_t, ndim=1] np_buf
  94     t_ = np_of_ptr(n, t, np.NPY_DOUBLE, 0)
  95     i_ = np_of_ptr(n, t, np.NPY_UINT64, 0)
  96     try:
  97          np_buf = (<object> hp.fn)(t_, i_).astype(np.float64)
  98     except Exception as exc:
  99         print exc # TODO
 100         return SD_ERR
 101     for j in xrange(n):
 102         buf[j] = np_buf[j]
 103     return SD_OK
 104
 105 cdef void hfill_py_free(sd_hfill *h):
 106     free(h)
 107
 108 cdef sd_hfill * hfill_py_new(object fn):
 109     cdef hfill_py *h = <hfill_py *> malloc(sizeof(hfill_py))
 110     if h==NULL:
 111         raise Exception('alloc hfill_py failed.')
 112     h.fn = <void*> fn
 113     h.hf.ptr = h
 114     h.hf.apply = &hfill_py_apply
 115     h.hf.free = &hfill_py_free
 116     return &(h.hf)
 117
 118 cdef class Hist(object):
 119     cdef sd_hist *h
 120     def __init__(self,
 121             np.ndarray[uint32_t, ndim=1, mode="c"] vi not None,
 122             np.ndarray[double, ndim=1, mode="c"] vd not None, double t0, double dt):
 123         self.h = sd_hist_new_default(vi.size, <uint32_t*> vi.data, <double*> vd.data, np.float64(t0), np.float64(dt))
 124     def __del__(self):
 125         if self.h != NULL:
 126             self.h.free(self.h)
 127     def fill(self, filler):
 128         cdef sd_hfill *hf = hfill_py_new(filler)
 129         cdef sd_stat stat = self.h.fill(self.h, hf)
 130         hf.free(hf)
 131         if stat == SD_ERR:
 132             raise Exception('fill failed.')
 133     def get(self, double t):
 134         cdef np.ndarray[double, ndim=1] aff = np.empty((self.nd(),), np.float64)
 135         self.h.get(self.h, t, <double*> aff.data)
 136         return aff
 137     def set(self, double t, np.ndarray[double, ndim=1] eff):
 138         # void sd_hist_set(sd_hist *  h, double t, double *  eff)
 139         self.h.set(self.h, t, <double*> eff.data)
 140     def maxvi(self):
 141         return self.h.get_maxvi(self.h)
 142     def nu(self):
 143         return self.h.get_nu(self.h)
 144     def nd(self):
 145         return self.h.get_nd(self.h)
 146     def t(self):
 147         return self.h.get_t(self.h)
 148     def dt(self):
 149         return self.h.get_dt(self.h)
 150     def vi2i(self, i):
 151         return self.h.get_vi2i(self.h, i)
 152     def buf_lin(self, i):
 153         return self.h.get_buf_lin(self.h, i)
 154     def lim(self, i):
 155         return self.h.get_lim(self.h, i)
 156     def len(self, i):
 157         return self.h.get_len(self.h, i)
 158     def pos(self, i):
 159         return self.h.get_pos(self.h, i)
 160     def uvi(self, i):
 161         return self.h.get_uvi(self.h, i)
 162     def maxd(self, i):
 163         return self.h.get_maxd(self.h, i)
 164     def vi(self, i):
 165         return self.h.get_vi(self.h, i)
 166     def vd(self, i):
 167         return self.h.get_vd(self.h, i)
 168
 169 # rng
 170 cdef extern from "sddekit.h":
 171     ctypedef struct sd_rng:
 172         void *ptr
 173         void (*seed)(sd_rng*, uint32_t seed)
 174         double (*norm)(sd_rng*)
 175         void (*fill_norm)(sd_rng*, uint32_t n, double *x)
 176         uint32_t (*nbytes)(sd_rng*)
 177         void (*free)(sd_rng*)
 178     sd_rng * sd_rng_new_default()
 179
 180 cdef class Rng(object):
 181     cdef sd_rng *r
 182
 183 cdef object rng_of_ptr(sd_rng *r):
 184     rng = Rng()
 185     rng.r = r
 186     return rng
 187
 188 # sys
 189 cdef extern from "sddekit.h":
 190     ctypedef struct sd_sys:
 191         uint32_t (*ndim)(sd_sys*)
 192         uint32_t (*ndc)(sd_sys*)
 193         uint32_t (*nobs)(sd_sys*)
 194         uint32_t (*nrpar)(sd_sys*)
 195         uint32_t (*nipar)(sd_sys*)
 196         sd_stat (*apply)(sd_sys*, sd_sys_in*, sd_sys_out*)
 197         void (*free)(sd_sys*)
 198         void *ptr
 199     ctypedef struct sd_sys_in:
 200         uint32_t nx, nc, id
 201         double t, *x, *i
 202         sd_hist *hist
 203         sd_rng *rng
 204     ctypedef struct sd_sys_out:
 205         double *f, *g, *o
 206
 207 cdef uint32_t sys_py_ndim(sd_sys*s):
 208     cdef object sys = <object> s.ptr
 209     cdef uint32_t i = sys.ndim
 210     return i
 211
 212 cdef uint32_t sys_py_ndc(sd_sys*s):
 213     cdef object sys = <object> s.ptr
 214     cdef uint32_t i = sys.ndc
 215     return i
 216
 217 cdef uint32_t sys_py_nobs(sd_sys*s):
 218     cdef object sys = <object> s.ptr
 219     cdef uint32_t i = sys.nobs
 220     return i
 221
 222 cdef uint32_t sys_py_nrpar(sd_sys*s):
 223     cdef object sys = <object> s.ptr
 224     cdef uint32_t i = sys.nrpar
 225     return i
 226
 227 cdef uint32_t sys_py_nipar(sd_sys*s):
 228     cdef object sys = <object> s.ptr
 229     cdef uint32_t i = sys.nipar
 230     return i
 231
 232 cdef sd_stat sys_py_apply(sd_sys *s, sd_sys_in *i, sd_sys_out *o):
 233     cdef object sys = <object> s.ptr
 234     # TODO conversions
 235     f, g, o = sys.apply(i.hist, i.rng, i.id, i.t, i.x, i.i)
 236     # TODO conversions
 237     o.f = f
 238     o.g = g
 239     o.o = o
 240     return SD_OK
 241
 242 cdef void sys_py_free(sd_sys *s):
 243     Py_DECREF(<object> s.ptr)
 244     free(s)
 245
 246 cdef class Sys(object):
 247     "Base class for system defintions, may be defined in C."
 248     cdef sd_sys *s
 249
 250     @property
 251     def ndim(self):
 252         return self.s.ndim(self.s)
 253
 254     @property
 255     def ndc(self):
 256         return self.s.ndc(self.s)
 257
 258     @property
 259     def nobs(self):
 260         return self.s.nobs(self.s)
 261
 262     @property
 263     def nrpar(self):
 264         return self.s.nrpar(self.s)
 265
 266     @property
 267     def nipar(self):
 268         return self.s.nipar(self.s)
 269
 270     def apply(self, hist, rng, id, t, x, i):
 271         cdef:
 272             sd_sys_in in_
 273             sd_sys_out out_
 274         raise NotImplementedError
 275
 276     def __del__(self):
 277         self.s.free(self.s)
 278
 279 cdef class PySys(Sys):
 280     "Base class for systems defined in Python."
 281
 282     def __init__(self):
 283         "Construct system defined by Python subclass."
 284         if self.__class__.__name__ == 'PySys':
 285             raise TypeError('PySys must be subclassed.')
 286         self.s = <sd_sys *> malloc(sizeof(sd_sys))
 287         Py_INCREF(self)
 288         self.s.ptr = <void*> self
 289         self.s.ndim = &sys_py_ndim
 290         self.s.ndc = &sys_py_ndc
 291         self.s.nobs = &sys_py_nobs
 292         self.s.nrpar = &sys_py_nrpar
 293         self.s.nipar = &sys_py_nipar
 294         self.s.apply = &sys_py_apply
 295         self.s.free = &sys_py_free
 296
 297     @property
 298     def ndim(self):
 299         raise NotImplementedError
 300
 301     @property
 302     def ndc(self):
 303         raise NotImplementedError
 304
 305     @property
 306     def nobs(self):
 307         raise NotImplementedError
 308
 309     @property
 310     def nrpar(self):
 311         raise NotImplementedError
 312
 313     @property
 314     def nipar(self):
 315         raise NotImplementedError
 316
 317     def apply(self, hist, rng, id, t, x, i):
 318         raise NotImplementedError
 319
 320 cdef object sys_of_ptr(sd_sys *s):
 321     sys = Sys()
 322     sys.s = s
 323     return sys
 324
 325
 326 # vim: sw=4 et