summaryrefslogtreecommitdiff
path: root/numpy
diff options
context:
space:
mode:
Diffstat (limited to 'numpy')
-rw-r--r--numpy/core/include/numpy/_neighborhood_iterator_imp.h10
-rw-r--r--numpy/core/src/multiarray/iterators.c114
-rw-r--r--numpy/core/src/multiarray/multiarray_tests.c.src317
-rw-r--r--numpy/core/tests/test_multiarray.py3
4 files changed, 424 insertions, 20 deletions
diff --git a/numpy/core/include/numpy/_neighborhood_iterator_imp.h b/numpy/core/include/numpy/_neighborhood_iterator_imp.h
index 37a03ebd1..b8679bc88 100644
--- a/numpy/core/include/numpy/_neighborhood_iterator_imp.h
+++ b/numpy/core/include/numpy/_neighborhood_iterator_imp.h
@@ -367,7 +367,8 @@ PyArrayNeighborhoodIter_Next(PyArrayNeighborhoodIterObject* iter)
iter->dataptr = iter->translate((PyArrayIterObject*)iter, iter->coordinates);
break;
case NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING:
- _PyArrayNeighborhoodIter_SetPtrMirror(iter);
+ //_PyArrayNeighborhoodIter_SetPtrMirror(iter);
+ iter->dataptr = iter->translate((PyArrayIterObject*)iter, iter->coordinates);
break;
case NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING:
_PyArrayNeighborhoodIter_SetPtrCircular(iter);
@@ -435,14 +436,14 @@ PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject* iter)
for (i = 0; i < iter->nd; ++i) {
iter->coordinates[i] = iter->bounds[i][0];
}
- iter->dataptr = _neigh_iter_get_ptr_const(iter);
+ iter->dataptr = iter->translate((PyArrayIterObject*)iter, iter->coordinates);
#if 0
switch (iter->mode) {
case NPY_NEIGHBORHOOD_ITER_ZERO_PADDING:
case NPY_NEIGHBORHOOD_ITER_ONE_PADDING:
case NPY_NEIGHBORHOOD_ITER_CONSTANT_PADDING:
- _PyArrayNeighborhoodIter_SetPtrConstant(iter);
+ iter->dataptr = _neigh_iter_get_ptr_const(iter);
break;
case NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING:
_PyArrayNeighborhoodIter_SetPtrMirror(iter);
@@ -450,6 +451,9 @@ PyArrayNeighborhoodIter_Reset(PyArrayNeighborhoodIterObject* iter)
case NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING:
_PyArrayNeighborhoodIter_SetPtrCircular(iter);
break;
+ default:
+ printf("oups\n");
+ exit(-1);
}
#endif
diff --git a/numpy/core/src/multiarray/iterators.c b/numpy/core/src/multiarray/iterators.c
index 7cdb199ab..4e42064cf 100644
--- a/numpy/core/src/multiarray/iterators.c
+++ b/numpy/core/src/multiarray/iterators.c
@@ -1813,10 +1813,120 @@ get_ptr_constant(PyArrayIterObject* _iter, npy_intp *coordinates)
_INF_SET_PTR(i)
}
+ // printf("%s: coordinates is %ld, ptr is %f\n", __func__, _coordinates[0], *((double*)p->translate(p, _coordinates)));
return p->translate(p, _coordinates);
}
#undef _INF_SET_PTR
+#define _NPY_IS_EVEN(x) ((x) % 2 == 0)
+
+/* For an array x of dimension n, and given index i, returns j, 0 <= j < n
+ * such as x[i] = x[j], with x assumed to be mirrored. For example, for x =
+ * {1, 2, 3} (n = 3)
+ *
+ * index -5 -4 -3 -2 -1 0 1 2 3 4 5 6
+ * value 2 3 3 2 1 1 2 3 3 2 1 1
+ *
+ * _npy_pos_index_mirror(4, 3) will return 1, because x[4] = x[1]*/
+static inline npy_intp
+__npy_pos_remainder(npy_intp i, npy_intp n)
+{
+ npy_intp k, l, j;
+
+ /* Mirror i such as it is guaranteed to be positive */
+ if (i < 0) {
+ i = - i - 1;
+ }
+
+ /* compute k and l such as i = k * n + l, 0 <= l < k */
+ k = i / n;
+ l = i - k * n;
+
+ if (_NPY_IS_EVEN(k)) {
+ j = l;
+ } else {
+ j = n - 1 - l;
+ }
+ return j;
+}
+#undef _NPY_IS_EVEN
+
+#define _INF_SET_PTR_MIRROR(c) \
+ printf("bounds: %ld - %ld\n", p->bounds[c][0], p->bounds[c][1]); \
+ bd = coordinates[c] + p->coordinates[c]; \
+ bd -= p->bounds[c][0]; \
+ printf("bd: %ld - max %ld\n", bd, p->bounds[c][1] - p->bounds[c][0]); \
+ truepos = __npy_pos_remainder(bd, niter->dimensions[c] - p->bounds[c][0]); \
+ printf("truepos: %ld \n", truepos); \
+ _coordinates[c] = (truepos + p->bounds[c][0]); \
+ printf("_cordinates: %ld \n", _coordinates[c]); \
+
+/* set the dataptr from its current coordinates */
+static char*
+get_ptr_mirror(PyArrayIterObject* _iter, npy_intp *coordinates)
+{
+ int i;
+ npy_intp bd, _coordinates[NPY_MAXDIMS];
+ npy_intp truepos;
+ PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter;
+ PyArrayIterObject *p = niter->_internal_iter;
+
+ //printf("%s\n", __func__);
+ for(i = 0; i < niter->nd; ++i) {
+ // _INF_SET_PTR_MIRROR(i)
+ // printf("bounds: %ld - %ld\n", p->bounds[i][0], p->bounds[i][1]);
+ bd = coordinates[i] + p->coordinates[i];
+ bd -= p->bounds[i][0];
+ // printf("bd: %ld - max %ld\n", bd, p->bounds[i][1] - p->bounds[i][0]);
+ truepos = __npy_pos_remainder(bd, p->bounds[i][1] + 1 - p->bounds[i][0]);
+ // printf("truepos: %ld \n", truepos);
+ _coordinates[i] = (truepos + p->bounds[i][0]);
+ // printf("_cordinates: %ld \n", _coordinates[i]);
+ }
+
+ // printf("%s: coordinates is %ld | %ld -> %ld\n", __func__, coordinates[0], p->coordinates[0], _coordinates[0]);
+ return p->translate(p, _coordinates);
+}
+#undef _INF_SET_PTR_MIRROR
+
+#if 0
+/* compute l such as i = k * n + l, 0 <= l < |k| */
+static inline npy_intp
+_npy_euclidean_division(npy_intp i, npy_intp n)
+{
+ npy_intp l;
+
+ l = i % n;
+ if (l < 0) {
+ l += n;
+ }
+ return l;
+}
+#endif
+
+#define _INF_SET_PTR_CIRCULAR(c) \
+ bd = coordinates[c] + niter->_internal_iter->coordinates[c]; \
+ truepos = _npy_euclidean_division(bd, niter->dimensions[c]); \
+ offset = (truepos - niter->_internal_iter->coordinates[c]) * niter->strides[c]; \
+ ret += offset;
+
+static char*
+get_coordinates_circular(PyArrayIterObject* _iter, npy_intp *coordinates)
+{
+ int i;
+ npy_intp offset, bd, truepos;
+ char *ret;
+ PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter;
+
+ ret = niter->_internal_iter->dataptr;
+
+ for(i = 0; i < niter->nd; ++i) {
+ _INF_SET_PTR_CIRCULAR(i)
+ }
+
+ return ret;
+}
+#undef _INF_SET_PTR_CIRCULAR
/*
* fill and x->ao should have equivalent types
@@ -1875,12 +1985,12 @@ PyArray_NeighborhoodIterNew(PyArrayIterObject *x, intp *bounds,
ret->mode = mode;
ret->translate = &get_ptr_constant;
break;
-#if 0
case NPY_NEIGHBORHOOD_ITER_MIRROR_PADDING:
ret->mode = mode;
ret->constant = NULL;
- ret->translate = get_coordinates_mirror;
+ ret->translate = &get_ptr_mirror;
break;
+#if 0
case NPY_NEIGHBORHOOD_ITER_CIRCULAR_PADDING:
ret->mode = mode;
ret->constant = NULL;
diff --git a/numpy/core/src/multiarray/multiarray_tests.c.src b/numpy/core/src/multiarray/multiarray_tests.c.src
index ba74d5de0..6b0887f3d 100644
--- a/numpy/core/src/multiarray/multiarray_tests.c.src
+++ b/numpy/core/src/multiarray/multiarray_tests.c.src
@@ -23,8 +23,22 @@ static int copy_@type@(PyArrayIterObject *itx, PyArrayNeighborhoodIterObject *ni
/* For each point in itx, copy the current neighborhood into an array which
* is appended at the output list */
for(i = 0; i < itx->size; ++i) {
+ npy_intp coord[10];
PyArrayNeighborhoodIter_Reset(niterx);
+#if 0
+ printf("+++++++++++++++++++++++++++++++\n");
+ coord[0] = -2;
+ printf("%f\n", *((double*)niterx->translate(niterx, coord)));
+ coord[0] = -1;
+ printf("%f\n", *((double*)niterx->translate(niterx, coord)));
+ coord[0] = 0;
+ printf("%f\n", *((double*)niterx->translate(niterx, coord)));
+ coord[0] = 1;
+ printf("%f\n", *((double*)niterx->translate(niterx, coord)));
+ coord[0] = 2;
+ printf("%f\n", *((double*)niterx->translate(niterx, coord)));
+#endif
for(j = 0; j < itx->ao->nd; ++j) {
odims[j] = bounds[2 * j + 1] - bounds[2 * j] + 1;
}
@@ -199,22 +213,217 @@ clean_ax:
return NULL;
}
+/*
+ * Update to next item of the iterator
+ *
+ * Note: this simply increment the coordinates vector, last dimension
+ * incremented first , i.e, for dimension 3
+ * ...
+ * -1, -1, -1
+ * -1, -1, 0
+ * -1, -1, 1
+ * ....
+ * -1, 0, -1
+ * -1, 0, 0
+ * ....
+ * 0, -1, -1
+ * 0, -1, 0
+ * ....
+ */
+#define _UPDATE_COORD_ITER(c) \
+ wb = iter->coordinates[c] < iter->bounds[c][1]; \
+ if (wb) { \
+ iter->coordinates[c] += 1; \
+ return 0; \
+ } \
+ else { \
+ iter->coordinates[c] = iter->bounds[c][0]; \
+ }
+
+static inline int
+neigh_iter_incr_coord(PyArrayNeighborhoodIterObject* iter)
+{
+ int i, wb;
+
+ for (i = iter->nd - 1; i >= 0; --i) {
+ _UPDATE_COORD_ITER(i)
+ }
+
+ return 0;
+}
+#undef _UPDATE_COORD_ITER
+
+#if 0
+#define _INF_SET_PTR(c) \
+ bd = niter->coordinates[c] + niter->_internal_iter->coordinates[c]; \
+ printf("%d|%d - %d|%d,\n", niter->coordinates[c], \
+ niter->_internal_iter->coordinates[c], \
+ niter->_internal_iter->bounds[c][0], \
+ niter->_internal_iter->bounds[c][1]); \
+ if (bd < niter->_internal_iter->bounds[c][0] || \
+ bd > niter->_internal_iter->bounds[c][1]) { \
+ printf("OOB\n"); \
+ return niter->constant; \
+ } \
+ coordinates[c] = bd; \
+ printf("coord is %d\n", coordinates[c]);
+#endif
+
+#define _INF_SET_PTR(c) \
+ bd = niter->coordinates[c] + p->coordinates[c]; \
+ if (bd < p->bounds[c][0] || \
+ bd > p->bounds[c][1]) { \
+ return niter->constant; \
+ } \
+ _coordinates[c] = bd;
+
+static inline char*
+neigh_iter_get_ptr_const(PyArrayNeighborhoodIterObject* niter)
+{
+ int i;
+ npy_intp offset, bd;
+ char *ret;
+ PyArrayIterObject *p = niter->_internal_iter;
+ npy_intp _coordinates[2*NPY_MAXDIMS];
+
+ for(i = 0; i < niter->nd; ++i) {
+ _INF_SET_PTR(i)
+ }
+
+ return niter->_internal_iter->translate(niter->_internal_iter, _coordinates);
+}
+
+static inline int
+neigh_iter_next_fptr(PyArrayNeighborhoodIterObject *iter)
+{
+ neigh_iter_incr_coord(iter);
+ iter->dataptr = neigh_iter_get_ptr_const(iter);
+ return 0;
+}
+
+static inline int
+neigh_iter_reset(PyArrayNeighborhoodIterObject *iter)
+{
+ npy_intp i;
+
+ for(i = 0; i < iter->nd; ++i) {
+ iter->coordinates[i] = iter->bounds[i][0];
+ }
+ iter->dataptr = neigh_iter_get_ptr_const(iter);
+ return 0;
+}
+
+//#include "cycle.h"
+
+//ticks T0, T1;
+double RES;
+double _MAX;
+
+#define _NPY_IS_EVEN(x) ((x) % 2 == 0)
+
+/* For an array x of dimension n, and given index i, returns j, 0 <= j < n
+ * such as x[i] = x[j], with x assumed to be mirrored. For example, for x =
+ * {1, 2, 3} (n = 3)
+ *
+ * index -5 -4 -3 -2 -1 0 1 2 3 4 5 6
+ * value 2 3 3 2 1 1 2 3 3 2 1 1
+ *
+ * _npy_pos_index_mirror(4, 3) will return 1, because x[4] = x[1]*/
+static inline npy_intp
+__npy_pos_remainder(npy_intp i, npy_intp n)
+{
+ npy_intp k, l, j;
+
+ /* Mirror i such as it is guaranteed to be positive */
+ if (i < 0) {
+ i = - i - 1;
+ }
+
+ /* compute k and l such as i = k * n + l, 0 <= l < k */
+ k = i / n;
+ l = i - k * n;
+
+ if (_NPY_IS_EVEN(k)) {
+ j = l;
+ } else {
+ j = n - 1 - l;
+ }
+ return j;
+}
+#undef _NPY_IS_EVEN
+
+#define _INF_SET_PTR_MIRROR(c) \
+ bd = coordinates[c] + p->coordinates[c]; \
+ bd -= p->bounds[c][0]; \
+ truepos = __npy_pos_remainder(bd, niter->dimensions[c] - p->bounds[c][0]); \
+ _coordinates[c] = (truepos - p->coordinates[c] + p->bounds[c][0]);
+
+/* set the dataptr from its current coordinates */
+static char*
+get_ptr_mirror(PyArrayIterObject* _iter, npy_intp *coordinates)
+{
+ int i;
+ npy_intp bd, _coordinates[NPY_MAXDIMS];
+ npy_intp truepos;
+ PyArrayNeighborhoodIterObject *niter = (PyArrayNeighborhoodIterObject*)_iter;
+ PyArrayIterObject *p = niter->_internal_iter;
+
+ printf("%s\n", __func__);
+ for(i = 0; i < niter->nd; ++i) {
+ _INF_SET_PTR_MIRROR(i)
+ }
+
+ printf("%s:%s: coordinates is %ld, ptr is %f\n", __FILE__, __func__, _coordinates[0], *((double*)p->translate(p, _coordinates)));
+ return p->translate(p, _coordinates);
+}
+#undef _INF_SET_PTR_MIRROR
+
static int
copy_double_double(PyArrayNeighborhoodIterObject *itx,
PyArrayNeighborhoodIterObject *niterx,
npy_intp *bounds,
PyObject **out)
{
- npy_intp i, j;
- double *ptr;
+ npy_intp i, j, k;
+ register double *ptr;
npy_intp odims[NPY_MAXDIMS];
PyArrayObject *aout;
/* For each point in itx, copy the current neighborhood into an array which
* is appended at the output list */
+ // itx->translate = &get_ptr_mirror;
+ PyArrayNeighborhoodIter_Reset(itx);
+ _MAX = 1e100;
for(i = 0; i < itx->size; ++i) {
- PyArrayNeighborhoodIter_Reset(niterx);
+ //for(i = 0; i < 1; ++i) {
+ npy_intp yo[10];
+ yo[0] = i;
+ printf("====================\n");
+ printf("(%f) \n", *((double*)itx->translate(itx, yo)));
+ // printf("(%f) \n", *((double*)itx->dataptr));
+ PyArrayNeighborhoodIter_Reset(niterx);
+ printf("------ begin --------------\n");
+ //T0 = getticks();
+ for(j = 0; j < niterx->size; ++j) {
+ printf("(%f) \n", *((double*)niterx->translate(niterx, niterx->coordinates)));
+ // printf("(%f) \n", *((double*)niterx->dataptr));
+ PyArrayNeighborhoodIter_Next(niterx);
+ }
+ printf("------ end -------------\n");
+
+ // yo[0] = 1;
+ // printf("(%f) \n", *((double*)itx->translate(itx, yo)));
+ // yo[0] = 2;
+ // printf("(%f) \n", *((double*)itx->translate(itx, yo)));
+ // yo[0] = 3;
+ // printf("(%f) \n", *((double*)itx->translate(itx, yo)));
+ // yo[0] = 4;
+ // printf("(%f) \n", *((double*)itx->translate(itx, yo)));
+
+ // yo[0] = -1;
+ // printf("(%f) \n", *((double*)niterx->translate(niterx, yo)));
+#if 0
for(j = 0; j < itx->ao->nd; ++j) {
odims[j] = bounds[2 * j + 1] - bounds[2 * j] + 1;
}
@@ -225,17 +434,63 @@ copy_double_double(PyArrayNeighborhoodIterObject *itx,
ptr = (double*)aout->data;
+ PyArrayNeighborhoodIter_Reset(niterx);
+ printf("===================\n");
+ T0 = getticks();
for(j = 0; j < niterx->size; ++j) {
+ npy_intp yo[10];
+ yo[0] = 0;
+ printf("-------------------\n");
+ niterx->_internal_iter->translate(niterx->_internal_iter, yo);
+ //ptr = (double*)niterx->dataptr;
+ //printf("%ld | %ld -> %p (%f) \n", niterx->coordinates[0],
+ //niterx->_internal_iter->coordinates[0],
+ //niterx->_internal_iter->translate(niterx->_internal_iter, yo),
+ //*((double*)niterx->translate(niterx, yo)));
+ //printf("%f\n", *((double*)niterx->dataptr));
*ptr = *((double*)niterx->dataptr);
+ //*ptr = *((double*)neigh_iter_get_ptr_const(niterx));
+
+ // coord[0] = -2;
+ // ptr = niterx->_internal_iter->translate(niterx->_internal_iter, coord);
+ // printf("%f \n", *ptr);
+
+ // neigh_iter_next_fptr(niterx);
PyArrayNeighborhoodIter_Next(niterx);
+ // neigh_iter_incr_coord(niterx);
ptr += 1;
}
+ T1 = getticks();
+ RES = elapsed(T1, T0);
+ if (RES < _MAX) {
+ _MAX = RES;
+ }
Py_INCREF(aout);
PyList_Append(*out, (PyObject*)aout);
Py_DECREF(aout);
+#endif
PyArrayNeighborhoodIter_Next(itx);
}
+ // printf("%f - %f\n", _MAX, _MAX / niterx->size);
+
+#if 0
+ {
+ npy_intp coord[10];
+ coord[0] = 0;
+ coord[1] = 0;
+ *ptr = *((double*)niterx->translate(niterx, coord));
+ printf("%f\n", *ptr);
+ coord[0] = 1;
+ coord[1] = 0;
+ *ptr = *((double*)niterx->translate(niterx, coord));
+ printf("%f\n", *ptr);
+ coord[0] = 2;
+ coord[1] = 0;
+ *ptr = *((double*)niterx->translate(niterx, coord));
+ printf("%f\n", *ptr);
+ }
+#endif
return 0;
}
@@ -243,14 +498,18 @@ copy_double_double(PyArrayNeighborhoodIterObject *itx,
static PyObject*
test_neighborhood_iterator_oob(PyObject* NPY_UNUSED(self), PyObject* args)
{
- PyObject *x, *out;
+ PyObject *x, *out, *b1, *b2;
PyArrayObject *ax;
PyArrayIterObject *itx;
- int i, typenum, mode, st;
+ int i, typenum, mode1, mode2, st;
npy_intp bounds[NPY_MAXDIMS*2];
PyArrayNeighborhoodIterObject *niterx1, *niterx2;
- if (!PyArg_ParseTuple(args, "O", &x)) {
+ if (!PyArg_ParseTuple(args, "OOiOi", &x, &b1, &mode1, &b2, &mode2)) {
+ return NULL;
+ }
+
+ if (!PySequence_Check(b1) || !PySequence_Check(b2)) {
return NULL;
}
@@ -260,6 +519,16 @@ test_neighborhood_iterator_oob(PyObject* NPY_UNUSED(self), PyObject* args)
if (ax == NULL) {
return NULL;
}
+ if (PySequence_Size(b1) != 2 * ax->nd) {
+ PyErr_SetString(PyExc_ValueError,
+ "bounds sequence 1 size not compatible with x input");
+ goto clean_ax;
+ }
+ if (PySequence_Size(b2) != 2 * ax->nd) {
+ PyErr_SetString(PyExc_ValueError,
+ "bounds sequence 2 size not compatible with x input");
+ goto clean_ax;
+ }
out = PyList_New(0);
if (out == NULL) {
@@ -272,27 +541,47 @@ test_neighborhood_iterator_oob(PyObject* NPY_UNUSED(self), PyObject* args)
}
/* Compute boundaries for the neighborhood iterator */
- for(i = 0; i < ax->nd; ++i) {
- bounds[2*i] = -2;
- bounds[2*i+1] = 1 + ax->dimensions[i];
+ for(i = 0; i < 2 * ax->nd; ++i) {
+ PyObject* bound;
+ bound = PySequence_GetItem(b1, i);
+ if (bounds == NULL) {
+ goto clean_itx;
+ }
+ if (!PyInt_Check(bound)) {
+ PyErr_SetString(PyExc_ValueError, "bound not long");
+ Py_DECREF(bound);
+ goto clean_itx;
+ }
+ bounds[i] = PyInt_AsLong(bound);
+ Py_DECREF(bound);
}
/* Create the neighborhood iterator */
niterx1 = (PyArrayNeighborhoodIterObject*)PyArray_NeighborhoodIterNew(
(PyArrayIterObject*)itx, bounds,
- NPY_NEIGHBORHOOD_ITER_ZERO_PADDING, NULL);
+ mode1, NULL);
if (niterx1 == NULL) {
goto clean_out;
}
- for(i = 0; i < 2 * niterx1->nd; ++i) {
- bounds[2*i] = 0;
- bounds[2*i+1] = 0;
+ for(i = 0; i < 2 * ax->nd; ++i) {
+ PyObject* bound;
+ bound = PySequence_GetItem(b2, i);
+ if (bounds == NULL) {
+ goto clean_itx;
+ }
+ if (!PyInt_Check(bound)) {
+ PyErr_SetString(PyExc_ValueError, "bound not long");
+ Py_DECREF(bound);
+ goto clean_itx;
+ }
+ bounds[i] = PyInt_AsLong(bound);
+ Py_DECREF(bound);
}
niterx2 = (PyArrayNeighborhoodIterObject*)PyArray_NeighborhoodIterNew(
(PyArrayIterObject*)niterx1, bounds,
- NPY_NEIGHBORHOOD_ITER_ZERO_PADDING, NULL);
+ mode2, NULL);
if (niterx1 == NULL) {
goto clean_niterx1;
}
diff --git a/numpy/core/tests/test_multiarray.py b/numpy/core/tests/test_multiarray.py
index ebce40c5d..15422d341 100644
--- a/numpy/core/tests/test_multiarray.py
+++ b/numpy/core/tests/test_multiarray.py
@@ -1244,7 +1244,8 @@ class TestNeighborhoodIter2(TestCase):
np.array([3], dtype=dt),
np.array([0], dtype=dt),
np.array([0], dtype=dt)]
- l = test_neighborhood_iterator_oob(x)
+ l = test_neighborhood_iterator_oob(x, [-2, 2], [0])
+ print l, r
assert_array_equal(l, r)
if __name__ == "__main__":
View Lorry Controller Status