一个 c 类型的函数，它在 cython 中返回一个指针数组答案

【问题标题】：a c-type function which returns a pointer array in cython一个 c 类型的函数，它在 cython 中返回一个指针数组
【发布时间】：2017-10-23 12:02:47
【问题描述】：

我一直在尝试编写一个返回 指针数组 (interv) 的 C 类型函数，我想在 python 类型函数（sampler）。我编写了以下代码，但它返回的结果是错误的。代码已编译，但doubling 函数返回到sampler 函数中的memoryview 对象的指针数组不正确。

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)

在python函数中将指针数组作为numpy数组返回的正确方法是什么？提前致谢。

【问题讨论】：

doubling 函数释放它返回的指针。你如何期望那个指针有数据？
@danny 真的吗？是这个原因吗？那么我应该在哪里释放指针的内存呢？
c malloc array pointer return in cython的可能重复

标签： python arrays numpy pointers cython

【解决方案1】：

在doubling 函数中，分配的内存被释放并返回刚刚释放的指针。在sampler 中使用时，该指针不再指向数据，因为它已被释放。

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

由于您正在使用 NumPy 数组，因此最好使用内存视图来传递指针，并使用 cython 数组来专门分配数据。

Cython 数组根据对象生命周期自动进行内存管理，而内存视图可以接受 cython 和/或 numpy 数组，而无需复制来替代指针和手动内存管理。

参见documentation for examples 和coercion of numpy arrays to cython arrays。

【讨论】：