CODEWITHSUNDEEP
pythonarraysnumpypointerscython
Dalek
Dalek

Reputation: 4318

a c-type function which returns a pointer array in cython

I have been trying to write a C-type function which returns a pointer array (interv) and I'd like to call this function inside a python-type function(sampler). I wrote the following code but the results that it returns are wrong. The code gets compiled, but the pointer array that doubling function returns to a memoryview object in the sampler function is not correct. 

from cpython cimport array
import cython
import numpy as np
import ctypes
cimport numpy as np
cimport cython
from libc.stdlib cimport malloc, free

from libcpp.vector cimport vector
cdef extern from "gsl/gsl_rng.h":#nogil:
     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 gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937)

cdef extern from "gsl/gsl_randist.h" nogil:
     double unif "gsl_rng_uniform"(gsl_rng * r)
     double unif_interval "gsl_ran_flat"(gsl_rng * r,double,double)    ## syntax; (seed, lower, upper)
     double exponential "gsl_ran_exponential"(gsl_rng * r,double) ## syntax; (seed, mean) ... mean is 1/rate


cdef double* doubling(double x0, double y, double w, int p):
     cdef double* interv = <double *>malloc(2 * cython.sizeof(double))
     if interv is NULL:
        raise MemoryError()
     cdef double u
     cdef int K    
     cdef bint now_left
     cdef double g_interv[2]
     u = unif_interval(r,0,1)
     interv[0] = x0 - w*u
     interv[1] = interv[0] + w
     if p>0:
        K = p
     g_interv[0]= f(interv[0])
     g_interv[1]= f(interv[1])
     while ((g_interv[0] > y) or (g_interv[1] > y)):
           u = unif_interval(r,0,1)             
           now_left = (u < 0.5)          
           if (now_left):
              interv[0] -= (interv[1] - interv[0])
              g_interv[0]=f(interv[0])
           else:
              interv[1] += (interv[1] - interv[0])
              g_interv[1]=f(interv[1])
           if p>0:
              K-=1
              if (K<=0):
                  break
     try:
         return interv
     finally:
         if interv is not NULL:
            free(interv)

def sampler(int n_sample,
                  int p = 0,
                  double x0=0.0, 
                  double w=0.1):
     cdef vector[double] samples
     cdef double vertical
     cdef Py_ssize_t i
     cdef np.ndarray[ndim=1, dtype=np.float64_t] interv
     cdef np.float64_t[:] view
     for 0<= i <n_sample:
              vertical = f(x0) - exponential(r, 1)      
              view=<np.float64_t[:2]>doubling(x0, vertical, w, p)
              interv= np.asarray(view)
              samples.push_back(interv[0])
     return samples

cdef double f(double x):
     cdef double a=5.
     return log(a)+(a-1.)*log(x)

What is the right way to return a pointer array as a numpy array in a python function? Thanks in advance.

Upvotes: 0

Views: 501

Answers (1)

danny
danny

Reputation: 5270

In the doubling function, the allocated memory is freed and the just-freed pointer returned. When used in sampler, that pointer no longer points to data as it has been freed.

 try:
     return interv
 finally:
     if interv is not NULL:
        free(interv)

Since you are working with NumPy arrays, probably best to use memory views for passing pointers around and cython arrays for allocating data exclusively.

Cython arrays automate memory management based on object life time while memory views can accept cython and/or numpy arrays without copying as a replacement for pointers and manual memory management.

See documentation for examples and also on coercion of numpy arrays to cython arrays.

Upvotes: 3

Related Questions