[Scipy-svn] r7165 - in trunk/scipy: interpolate spatial
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Feb 20 11:53:11 EST 2011
Author: ptvirtan
Date: 2011-02-20 10:53:10 -0600 (Sun, 20 Feb 2011)
New Revision: 7165
Modified:
trunk/scipy/interpolate/interpnd.pyx
trunk/scipy/spatial/qhull.pxd
trunk/scipy/spatial/qhull.pyx
Log:
BUG: spatial/interpnd: never deallocate memory in a different Cython module
Different Cython modules end up as different DLL/.so files, and
apparently on some platforms freeing memory allocated in a different DLL
is not allowed. This seems to be true at least on Windows 7 + 32-bit
Python.
Modified: trunk/scipy/interpolate/interpnd.pyx
===================================================================
--- trunk/scipy/interpolate/interpnd.pyx 2011-02-20 13:30:33 UTC (rev 7164)
+++ trunk/scipy/interpolate/interpnd.pyx 2011-02-20 16:53:10 UTC (rev 7165)
@@ -31,10 +31,6 @@
# Numpy etc.
#------------------------------------------------------------------------------
-cdef extern from "stdlib.h":
- void *malloc(int size) nogil
- void free(void *ptr) nogil
-
cdef extern from "math.h":
double sqrt(double x) nogil
double fabs(double a) nogil
@@ -164,7 +160,7 @@
"""
LinearNDInterpolator(points, values)
- Piecewise linear interpolant in N dimensions.
+ Piecewise linear interpolant in N dimensions.
.. versionadded:: 0.9
@@ -205,13 +201,13 @@
cdef double c[NPY_MAXDIMS]
cdef ${CDTYPE} fill_value
cdef int i, j, k, m, ndim, isimplex, inside, start, nvalues
- cdef qhull.DelaunayInfo_t *info
+ cdef qhull.DelaunayInfo_t info
ndim = xi.shape[1]
start = 0
fill_value = self.fill_value
- info = qhull._get_delaunay_info(self.tri, 1, 0)
+ qhull._get_delaunay_info(&info, self.tri, 1, 0)
out = np.zeros((xi.shape[0], self.values.shape[1]), dtype=np.${DTYPE})
nvalues = out.shape[1]
@@ -223,7 +219,7 @@
# 1) Find the simplex
- isimplex = qhull._find_simplex(info, c,
+ isimplex = qhull._find_simplex(&info, c,
(<double*>xi.data) + i*ndim,
&start, eps)
@@ -258,7 +254,6 @@
out[i,k].imag += c[j] * values[m, k].imag
% endif
- free(<void*>info)
return out
% endfor
@@ -431,7 +426,7 @@
change = max(fabs(y[it.vertex*2 + 0] + r[0]),
fabs(y[it.vertex*2 + 1] + r[1]))
-
+
y[it.vertex*2 + 0] = -r[0]
y[it.vertex*2 + 1] = -r[1]
@@ -448,7 +443,7 @@
def estimate_gradients_2d_global(tri, y, maxiter=400, tol=1e-6):
cdef np.ndarray[np.double_t, ndim=2] data
cdef np.ndarray[np.double_t, ndim=3] grad
- cdef qhull.DelaunayInfo_t *info
+ cdef qhull.DelaunayInfo_t info
cdef int k, ret, nvalues
y = np.asanyarray(y)
@@ -476,13 +471,13 @@
data = y
grad = yi
- info = qhull._get_delaunay_info(tri, 0, 1)
+ qhull._get_delaunay_info(&info, tri, 0, 1)
nvalues = data.shape[0]
for k in xrange(nvalues):
with nogil:
ret = _estimate_gradients_2d_global(
- info,
+ &info,
<double*>data.data + info.npoints*k,
maxiter,
tol,
@@ -493,7 +488,6 @@
"the results may be inaccurate",
GradientEstimationWarning)
- free(info)
return yi.transpose(1, 0, 2).reshape(y_shape + (2,))
@@ -742,7 +736,7 @@
"""
CloughTocher2DInterpolator(points, values, tol=1e-6)
- Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
+ Piecewise cubic, C1 smooth, curvature-minimizing interpolant in 2D.
.. versionadded:: 0.9
@@ -821,13 +815,13 @@
cdef ${CDTYPE} w
cdef ${CDTYPE} fill_value
cdef int i, j, k, m, ndim, isimplex, inside, start, nvalues
- cdef qhull.DelaunayInfo_t *info
+ cdef qhull.DelaunayInfo_t info
ndim = xi.shape[1]
start = 0
fill_value = self.fill_value
- info = qhull._get_delaunay_info(self.tri, 1, 1)
+ qhull._get_delaunay_info(&info, self.tri, 1, 1)
out = np.zeros((xi.shape[0], self.values.shape[1]), dtype=np.${DTYPE})
nvalues = out.shape[1]
@@ -838,7 +832,7 @@
for i in xrange(xi.shape[0]):
# 1) Find the simplex
- isimplex = qhull._find_simplex(info, c,
+ isimplex = qhull._find_simplex(&info, c,
(<double*>xi.data) + i*ndim,
&start, eps)
@@ -870,7 +864,7 @@
df[2*j+1].imag = grad[vertices[isimplex,j],k,1].imag
% endif
- w = _clough_tocher_2d_single_${DTYPE}(info, isimplex, c,
+ w = _clough_tocher_2d_single_${DTYPE}(&info, isimplex, c,
f, df)
% if DTYPE == "double":
out[i,k] = w
@@ -879,7 +873,6 @@
out[i,k].imag = w.imag
% endif
- free(<void*>info)
return out
% endfor
Modified: trunk/scipy/spatial/qhull.pxd
===================================================================
--- trunk/scipy/spatial/qhull.pxd 2011-02-20 13:30:33 UTC (rev 7164)
+++ trunk/scipy/spatial/qhull.pxd 2011-02-20 16:53:10 UTC (rev 7165)
@@ -9,10 +9,6 @@
# Distributed under the same BSD license as Scipy.
#
-cdef extern from "stdlib.h":
- void *malloc(int size)
- void free(void *ptr)
-
cdef extern from "numpy/ndarrayobject.h":
cdef enum:
NPY_MAXDIMS
@@ -32,8 +28,9 @@
double *max_bound
double *min_bound
-cdef DelaunayInfo_t *_get_delaunay_info(obj, int compute_transform,
- int compute_vertex_to_simplex)
+cdef void _get_delaunay_info(DelaunayInfo_t *, obj,
+ int compute_transform,
+ int compute_vertex_to_simplex)
#
# N-D geometry
Modified: trunk/scipy/spatial/qhull.pyx
===================================================================
--- trunk/scipy/spatial/qhull.pyx 2011-02-20 13:30:33 UTC (rev 7164)
+++ trunk/scipy/spatial/qhull.pyx 2011-02-20 16:53:10 UTC (rev 7165)
@@ -593,22 +593,15 @@
cdef class RidgeIter2D(object):
cdef RidgeIter2D_t it
cdef object delaunay
- cdef DelaunayInfo_t *info
+ cdef DelaunayInfo_t info
def __init__(self, delaunay, ivertex):
- self.info = NULL
if delaunay.ndim != 2:
raise ValueError("RidgeIter2D supports only 2-D")
self.delaunay = delaunay
- self.info = _get_delaunay_info(delaunay, 0, 1)
- _RidgeIter2D_init(&self.it, self.info, ivertex)
+ _get_delaunay_info(&self.info, delaunay, 0, 1)
+ _RidgeIter2D_init(&self.it, &self.info, ivertex)
- def __del__(self):
- if self.info != NULL:
- free(self.info)
- self.info = NULL
- self.delaunay = None
-
def __iter__(self):
return self
@@ -1067,7 +1060,7 @@
directed search in N dimensions.
"""
- cdef DelaunayInfo_t *info
+ cdef DelaunayInfo_t info
cdef int isimplex
cdef double c[NPY_MAXDIMS]
cdef double eps
@@ -1090,23 +1083,21 @@
eps = np.finfo(np.double).eps * 10
out = np.zeros((xi.shape[0],), dtype=np.intc)
out_ = out
- info = _get_delaunay_info(self, 1, 0)
+ _get_delaunay_info(&info, self, 1, 0)
if bruteforce:
for k in xrange(x.shape[0]):
isimplex = _find_simplex_bruteforce(
- info, c,
+ &info, c,
<double*>x.data + info.ndim*k,
eps)
out_[k] = isimplex
else:
for k in xrange(x.shape[0]):
- isimplex = _find_simplex(info, c, <double*>x.data + info.ndim*k,
+ isimplex = _find_simplex(&info, c, <double*>x.data + info.ndim*k,
&start, eps)
out_[k] = isimplex
- free(info)
-
return out.reshape(xi_shape[:-1])
def plane_distance(self, xi):
@@ -1118,7 +1109,7 @@
"""
cdef np.ndarray[np.double_t, ndim=2] x
cdef np.ndarray[np.double_t, ndim=2] out_
- cdef DelaunayInfo_t *info
+ cdef DelaunayInfo_t info
cdef double z[NPY_MAXDIMS+1]
cdef int i, j, k
@@ -1130,18 +1121,16 @@
xi = xi.reshape(np.prod(xi.shape[:-1]), xi.shape[-1])
x = np.ascontiguousarray(xi.astype(np.double))
- info = _get_delaunay_info(self, 0, 0)
+ _get_delaunay_info(&info, self, 0, 0)
out = np.zeros((x.shape[0], info.nsimplex), dtype=np.double)
out_ = out
for i in xrange(x.shape[0]):
for j in xrange(info.nsimplex):
- _lift_point(info, (<double*>x.data) + info.ndim*i, z)
- out_[i,j] = _distplane(info, j, z)
+ _lift_point(&info, (<double*>x.data) + info.ndim*i, z)
+ out_[i,j] = _distplane(&info, j, z)
- free(info)
-
return out.reshape(xi_shape[:-1] + (self.nsimplex,))
def lift_points(tri, x):
@@ -1180,10 +1169,10 @@
# Delaunay triangulation interface, for low-level C
#------------------------------------------------------------------------------
-cdef DelaunayInfo_t *_get_delaunay_info(obj,
- int compute_transform,
- int compute_vertex_to_simplex):
- cdef DelaunayInfo_t *info
+cdef void _get_delaunay_info(DelaunayInfo_t *info,
+ obj,
+ int compute_transform,
+ int compute_vertex_to_simplex):
cdef np.ndarray[np.double_t, ndim=3] transform
cdef np.ndarray[np.npy_int, ndim=1] vertex_to_simplex
cdef np.ndarray[np.double_t, ndim=2] points = obj.points
@@ -1193,7 +1182,6 @@
cdef np.ndarray[np.double_t, ndim=1] min_bound = obj.min_bound
cdef np.ndarray[np.double_t, ndim=1] max_bound = obj.max_bound
- info = <DelaunayInfo_t*>malloc(sizeof(DelaunayInfo_t))
info.ndim = points.shape[1]
info.npoints = points.shape[0]
info.nsimplex = vertices.shape[0]
@@ -1215,5 +1203,3 @@
info.vertex_to_simplex = NULL
info.min_bound = <double*>min_bound.data
info.max_bound = <double*>max_bound.data
-
- return info
More information about the Scipy-svn
mailing list