1 #define NPY_NO_DEPRECATED_API NPY_API_VERSION
2 #include <Python.h>
3 #include <structmember.h>
4 #include <string.h>
5
6 #define _MULTIARRAYMODULE
7 #include "numpy/arrayobject.h"
8 #include "numpy/npy_3kcompat.h"
9 #include "numpy/npy_math.h"
10 #include "npy_config.h"
11 #include "templ_common.h" /* for npy_mul_with_overflow_intp */
12 #include "lowlevel_strided_loops.h" /* for npy_bswap8 */
13 #include "alloc.h"
14 #include "ctors.h"
15 #include "common.h"
16
17
18 /*
19 * Returns -1 if the array is monotonic decreasing,
20 * +1 if the array is monotonic increasing,
21 * and 0 if the array is not monotonic.
22 */
23 static int
check_array_monotonic(const double * a,npy_intp lena)24 check_array_monotonic(const double *a, npy_intp lena)
25 {
26 npy_intp i;
27 double next;
28 double last;
29
30 if (lena == 0) {
31 /* all bin edges hold the same value */
32 return 1;
33 }
34 last = a[0];
35
36 /* Skip repeated values at the beginning of the array */
37 for (i = 1; (i < lena) && (a[i] == last); i++);
38
39 if (i == lena) {
40 /* all bin edges hold the same value */
41 return 1;
42 }
43
44 next = a[i];
45 if (last < next) {
46 /* Possibly monotonic increasing */
47 for (i += 1; i < lena; i++) {
48 last = next;
49 next = a[i];
50 if (last > next) {
51 return 0;
52 }
53 }
54 return 1;
55 }
56 else {
57 /* last > next, possibly monotonic decreasing */
58 for (i += 1; i < lena; i++) {
59 last = next;
60 next = a[i];
61 if (last < next) {
62 return 0;
63 }
64 }
65 return -1;
66 }
67 }
68
69 /* Find the minimum and maximum of an integer array */
70 static void
minmax(const npy_intp * data,npy_intp data_len,npy_intp * mn,npy_intp * mx)71 minmax(const npy_intp *data, npy_intp data_len, npy_intp *mn, npy_intp *mx)
72 {
73 npy_intp min = *data;
74 npy_intp max = *data;
75
76 while (--data_len) {
77 const npy_intp val = *(++data);
78 if (val < min) {
79 min = val;
80 }
81 else if (val > max) {
82 max = val;
83 }
84 }
85
86 *mn = min;
87 *mx = max;
88 }
89
90 /*
91 * arr_bincount is registered as bincount.
92 *
93 * bincount accepts one, two or three arguments. The first is an array of
94 * non-negative integers The second, if present, is an array of weights,
95 * which must be promotable to double. Call these arguments list and
96 * weight. Both must be one-dimensional with len(weight) == len(list). If
97 * weight is not present then bincount(list)[i] is the number of occurrences
98 * of i in list. If weight is present then bincount(self,list, weight)[i]
99 * is the sum of all weight[j] where list [j] == i. Self is not used.
100 * The third argument, if present, is a minimum length desired for the
101 * output array.
102 */
103 NPY_NO_EXPORT PyObject *
arr_bincount(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwds)104 arr_bincount(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
105 {
106 PyObject *list = NULL, *weight = Py_None, *mlength = NULL;
107 PyArrayObject *lst = NULL, *ans = NULL, *wts = NULL;
108 npy_intp *numbers, *ians, len, mx, mn, ans_size;
109 npy_intp minlength = 0;
110 npy_intp i;
111 double *weights , *dans;
112 static char *kwlist[] = {"list", "weights", "minlength", NULL};
113
114 if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|OO:bincount",
115 kwlist, &list, &weight, &mlength)) {
116 goto fail;
117 }
118
119 lst = (PyArrayObject *)PyArray_ContiguousFromAny(list, NPY_INTP, 1, 1);
120 if (lst == NULL) {
121 goto fail;
122 }
123 len = PyArray_SIZE(lst);
124
125 /*
126 * This if/else if can be removed by changing the argspec to O|On above,
127 * once we retire the deprecation
128 */
129 if (mlength == Py_None) {
130 /* NumPy 1.14, 2017-06-01 */
131 if (DEPRECATE("0 should be passed as minlength instead of None; "
132 "this will error in future.") < 0) {
133 goto fail;
134 }
135 }
136 else if (mlength != NULL) {
137 minlength = PyArray_PyIntAsIntp(mlength);
138 if (error_converting(minlength)) {
139 goto fail;
140 }
141 }
142
143 if (minlength < 0) {
144 PyErr_SetString(PyExc_ValueError,
145 "'minlength' must not be negative");
146 goto fail;
147 }
148
149 /* handle empty list */
150 if (len == 0) {
151 ans = (PyArrayObject *)PyArray_ZEROS(1, &minlength, NPY_INTP, 0);
152 if (ans == NULL){
153 goto fail;
154 }
155 Py_DECREF(lst);
156 return (PyObject *)ans;
157 }
158
159 numbers = (npy_intp *)PyArray_DATA(lst);
160 minmax(numbers, len, &mn, &mx);
161 if (mn < 0) {
162 PyErr_SetString(PyExc_ValueError,
163 "'list' argument must have no negative elements");
164 goto fail;
165 }
166 ans_size = mx + 1;
167 if (mlength != Py_None) {
168 if (ans_size < minlength) {
169 ans_size = minlength;
170 }
171 }
172 if (weight == Py_None) {
173 ans = (PyArrayObject *)PyArray_ZEROS(1, &ans_size, NPY_INTP, 0);
174 if (ans == NULL) {
175 goto fail;
176 }
177 ians = (npy_intp *)PyArray_DATA(ans);
178 NPY_BEGIN_ALLOW_THREADS;
179 for (i = 0; i < len; i++)
180 ians[numbers[i]] += 1;
181 NPY_END_ALLOW_THREADS;
182 Py_DECREF(lst);
183 }
184 else {
185 wts = (PyArrayObject *)PyArray_ContiguousFromAny(
186 weight, NPY_DOUBLE, 1, 1);
187 if (wts == NULL) {
188 goto fail;
189 }
190 weights = (double *)PyArray_DATA(wts);
191 if (PyArray_SIZE(wts) != len) {
192 PyErr_SetString(PyExc_ValueError,
193 "The weights and list don't have the same length.");
194 goto fail;
195 }
196 ans = (PyArrayObject *)PyArray_ZEROS(1, &ans_size, NPY_DOUBLE, 0);
197 if (ans == NULL) {
198 goto fail;
199 }
200 dans = (double *)PyArray_DATA(ans);
201 NPY_BEGIN_ALLOW_THREADS;
202 for (i = 0; i < len; i++) {
203 dans[numbers[i]] += weights[i];
204 }
205 NPY_END_ALLOW_THREADS;
206 Py_DECREF(lst);
207 Py_DECREF(wts);
208 }
209 return (PyObject *)ans;
210
211 fail:
212 Py_XDECREF(lst);
213 Py_XDECREF(wts);
214 Py_XDECREF(ans);
215 return NULL;
216 }
217
218 /* Internal function to expose check_array_monotonic to python */
219 NPY_NO_EXPORT PyObject *
arr__monotonicity(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwds)220 arr__monotonicity(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
221 {
222 static char *kwlist[] = {"x", NULL};
223 PyObject *obj_x = NULL;
224 PyArrayObject *arr_x = NULL;
225 long monotonic;
226 npy_intp len_x;
227 NPY_BEGIN_THREADS_DEF;
228
229 if (!PyArg_ParseTupleAndKeywords(args, kwds, "O|_monotonicity", kwlist,
230 &obj_x)) {
231 return NULL;
232 }
233
234 /*
235 * TODO:
236 * `x` could be strided, needs change to check_array_monotonic
237 * `x` is forced to double for this check
238 */
239 arr_x = (PyArrayObject *)PyArray_FROMANY(
240 obj_x, NPY_DOUBLE, 1, 1, NPY_ARRAY_CARRAY_RO);
241 if (arr_x == NULL) {
242 return NULL;
243 }
244
245 len_x = PyArray_SIZE(arr_x);
246 NPY_BEGIN_THREADS_THRESHOLDED(len_x)
247 monotonic = check_array_monotonic(
248 (const double *)PyArray_DATA(arr_x), len_x);
249 NPY_END_THREADS
250 Py_DECREF(arr_x);
251
252 return PyLong_FromLong(monotonic);
253 }
254
255 /*
256 * Returns input array with values inserted sequentially into places
257 * indicated by the mask
258 */
259 NPY_NO_EXPORT PyObject *
arr_insert(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwdict)260 arr_insert(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
261 {
262 char *src, *dest;
263 npy_bool *mask_data;
264 PyArray_Descr *dtype;
265 PyArray_CopySwapFunc *copyswap;
266 PyObject *array0, *mask0, *values0;
267 PyArrayObject *array, *mask, *values;
268 npy_intp i, j, chunk, nm, ni, nv;
269
270 static char *kwlist[] = {"input", "mask", "vals", NULL};
271 NPY_BEGIN_THREADS_DEF;
272 values = mask = NULL;
273
274 if (!PyArg_ParseTupleAndKeywords(args, kwdict, "O!OO:place", kwlist,
275 &PyArray_Type, &array0, &mask0, &values0)) {
276 return NULL;
277 }
278
279 array = (PyArrayObject *)PyArray_FromArray((PyArrayObject *)array0, NULL,
280 NPY_ARRAY_CARRAY | NPY_ARRAY_WRITEBACKIFCOPY);
281 if (array == NULL) {
282 goto fail;
283 }
284
285 ni = PyArray_SIZE(array);
286 dest = PyArray_DATA(array);
287 chunk = PyArray_DESCR(array)->elsize;
288 mask = (PyArrayObject *)PyArray_FROM_OTF(mask0, NPY_BOOL,
289 NPY_ARRAY_CARRAY | NPY_ARRAY_FORCECAST);
290 if (mask == NULL) {
291 goto fail;
292 }
293
294 nm = PyArray_SIZE(mask);
295 if (nm != ni) {
296 PyErr_SetString(PyExc_ValueError,
297 "place: mask and data must be "
298 "the same size");
299 goto fail;
300 }
301
302 mask_data = PyArray_DATA(mask);
303 dtype = PyArray_DESCR(array);
304 Py_INCREF(dtype);
305
306 values = (PyArrayObject *)PyArray_FromAny(values0, dtype,
307 0, 0, NPY_ARRAY_CARRAY, NULL);
308 if (values == NULL) {
309 goto fail;
310 }
311
312 nv = PyArray_SIZE(values); /* zero if null array */
313 if (nv <= 0) {
314 npy_bool allFalse = 1;
315 i = 0;
316
317 while (allFalse && i < ni) {
318 if (mask_data[i]) {
319 allFalse = 0;
320 } else {
321 i++;
322 }
323 }
324 if (!allFalse) {
325 PyErr_SetString(PyExc_ValueError,
326 "Cannot insert from an empty array!");
327 goto fail;
328 } else {
329 Py_XDECREF(values);
330 Py_XDECREF(mask);
331 PyArray_ResolveWritebackIfCopy(array);
332 Py_XDECREF(array);
333 Py_RETURN_NONE;
334 }
335 }
336
337 src = PyArray_DATA(values);
338 j = 0;
339
340 copyswap = PyArray_DESCR(array)->f->copyswap;
341 NPY_BEGIN_THREADS_DESCR(PyArray_DESCR(array));
342 for (i = 0; i < ni; i++) {
343 if (mask_data[i]) {
344 if (j >= nv) {
345 j = 0;
346 }
347
348 copyswap(dest + i*chunk, src + j*chunk, 0, array);
349 j++;
350 }
351 }
352 NPY_END_THREADS;
353
354 Py_XDECREF(values);
355 Py_XDECREF(mask);
356 PyArray_ResolveWritebackIfCopy(array);
357 Py_DECREF(array);
358 Py_RETURN_NONE;
359
360 fail:
361 Py_XDECREF(mask);
362 PyArray_ResolveWritebackIfCopy(array);
363 Py_XDECREF(array);
364 Py_XDECREF(values);
365 return NULL;
366 }
367
368 #define LIKELY_IN_CACHE_SIZE 8
369
370 #ifdef __INTEL_COMPILER
371 #pragma intel optimization_level 0
372 #endif
373 static NPY_INLINE npy_intp
_linear_search(const npy_double key,const npy_double * arr,const npy_intp len,const npy_intp i0)374 _linear_search(const npy_double key, const npy_double *arr, const npy_intp len, const npy_intp i0)
375 {
376 npy_intp i;
377
378 for (i = i0; i < len && key >= arr[i]; i++);
379 return i - 1;
380 }
381
382 /** @brief find index of a sorted array such that arr[i] <= key < arr[i + 1].
383 *
384 * If an starting index guess is in-range, the array values around this
385 * index are first checked. This allows for repeated calls for well-ordered
386 * keys (a very common case) to use the previous index as a very good guess.
387 *
388 * If the guess value is not useful, bisection of the array is used to
389 * find the index. If there is no such index, the return values are:
390 * key < arr[0] -- -1
391 * key == arr[len - 1] -- len - 1
392 * key > arr[len - 1] -- len
393 * The array is assumed contiguous and sorted in ascending order.
394 *
395 * @param key key value.
396 * @param arr contiguous sorted array to be searched.
397 * @param len length of the array.
398 * @param guess initial guess of index
399 * @return index
400 */
401 static npy_intp
binary_search_with_guess(const npy_double key,const npy_double * arr,npy_intp len,npy_intp guess)402 binary_search_with_guess(const npy_double key, const npy_double *arr,
403 npy_intp len, npy_intp guess)
404 {
405 npy_intp imin = 0;
406 npy_intp imax = len;
407
408 /* Handle keys outside of the arr range first */
409 if (key > arr[len - 1]) {
410 return len;
411 }
412 else if (key < arr[0]) {
413 return -1;
414 }
415
416 /*
417 * If len <= 4 use linear search.
418 * From above we know key >= arr[0] when we start.
419 */
420 if (len <= 4) {
421 return _linear_search(key, arr, len, 1);
422 }
423
424 if (guess > len - 3) {
425 guess = len - 3;
426 }
427 if (guess < 1) {
428 guess = 1;
429 }
430
431 /* check most likely values: guess - 1, guess, guess + 1 */
432 if (key < arr[guess]) {
433 if (key < arr[guess - 1]) {
434 imax = guess - 1;
435 /* last attempt to restrict search to items in cache */
436 if (guess > LIKELY_IN_CACHE_SIZE &&
437 key >= arr[guess - LIKELY_IN_CACHE_SIZE]) {
438 imin = guess - LIKELY_IN_CACHE_SIZE;
439 }
440 }
441 else {
442 /* key >= arr[guess - 1] */
443 return guess - 1;
444 }
445 }
446 else {
447 /* key >= arr[guess] */
448 if (key < arr[guess + 1]) {
449 return guess;
450 }
451 else {
452 /* key >= arr[guess + 1] */
453 if (key < arr[guess + 2]) {
454 return guess + 1;
455 }
456 else {
457 /* key >= arr[guess + 2] */
458 imin = guess + 2;
459 /* last attempt to restrict search to items in cache */
460 if (guess < len - LIKELY_IN_CACHE_SIZE - 1 &&
461 key < arr[guess + LIKELY_IN_CACHE_SIZE]) {
462 imax = guess + LIKELY_IN_CACHE_SIZE;
463 }
464 }
465 }
466 }
467
468 /* finally, find index by bisection */
469 while (imin < imax) {
470 const npy_intp imid = imin + ((imax - imin) >> 1);
471 if (key >= arr[imid]) {
472 imin = imid + 1;
473 }
474 else {
475 imax = imid;
476 }
477 }
478 return imin - 1;
479 }
480
481 #undef LIKELY_IN_CACHE_SIZE
482
483 NPY_NO_EXPORT PyObject *
arr_interp(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwdict)484 arr_interp(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
485 {
486
487 PyObject *fp, *xp, *x;
488 PyObject *left = NULL, *right = NULL;
489 PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
490 npy_intp i, lenx, lenxp;
491 npy_double lval, rval;
492 const npy_double *dy, *dx, *dz;
493 npy_double *dres, *slopes = NULL;
494
495 static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL};
496
497 NPY_BEGIN_THREADS_DEF;
498
499 if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO:interp", kwlist,
500 &x, &xp, &fp, &left, &right)) {
501 return NULL;
502 }
503
504 afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_DOUBLE, 1, 1);
505 if (afp == NULL) {
506 return NULL;
507 }
508 axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1);
509 if (axp == NULL) {
510 goto fail;
511 }
512 ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 0, 0);
513 if (ax == NULL) {
514 goto fail;
515 }
516 lenxp = PyArray_SIZE(axp);
517 if (lenxp == 0) {
518 PyErr_SetString(PyExc_ValueError,
519 "array of sample points is empty");
520 goto fail;
521 }
522 if (PyArray_SIZE(afp) != lenxp) {
523 PyErr_SetString(PyExc_ValueError,
524 "fp and xp are not of the same length.");
525 goto fail;
526 }
527
528 af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax),
529 PyArray_DIMS(ax), NPY_DOUBLE);
530 if (af == NULL) {
531 goto fail;
532 }
533 lenx = PyArray_SIZE(ax);
534
535 dy = (const npy_double *)PyArray_DATA(afp);
536 dx = (const npy_double *)PyArray_DATA(axp);
537 dz = (const npy_double *)PyArray_DATA(ax);
538 dres = (npy_double *)PyArray_DATA(af);
539 /* Get left and right fill values. */
540 if ((left == NULL) || (left == Py_None)) {
541 lval = dy[0];
542 }
543 else {
544 lval = PyFloat_AsDouble(left);
545 if (error_converting(lval)) {
546 goto fail;
547 }
548 }
549 if ((right == NULL) || (right == Py_None)) {
550 rval = dy[lenxp - 1];
551 }
552 else {
553 rval = PyFloat_AsDouble(right);
554 if (error_converting(rval)) {
555 goto fail;
556 }
557 }
558
559 /* binary_search_with_guess needs at least a 3 item long array */
560 if (lenxp == 1) {
561 const npy_double xp_val = dx[0];
562 const npy_double fp_val = dy[0];
563
564 NPY_BEGIN_THREADS_THRESHOLDED(lenx);
565 for (i = 0; i < lenx; ++i) {
566 const npy_double x_val = dz[i];
567 dres[i] = (x_val < xp_val) ? lval :
568 ((x_val > xp_val) ? rval : fp_val);
569 }
570 NPY_END_THREADS;
571 }
572 else {
573 npy_intp j = 0;
574
575 /* only pre-calculate slopes if there are relatively few of them. */
576 if (lenxp <= lenx) {
577 slopes = PyArray_malloc((lenxp - 1) * sizeof(npy_double));
578 if (slopes == NULL) {
579 PyErr_NoMemory();
580 goto fail;
581 }
582 }
583
584 NPY_BEGIN_THREADS;
585
586 if (slopes != NULL) {
587 for (i = 0; i < lenxp - 1; ++i) {
588 slopes[i] = (dy[i+1] - dy[i]) / (dx[i+1] - dx[i]);
589 }
590 }
591
592 for (i = 0; i < lenx; ++i) {
593 const npy_double x_val = dz[i];
594
595 if (npy_isnan(x_val)) {
596 dres[i] = x_val;
597 continue;
598 }
599
600 j = binary_search_with_guess(x_val, dx, lenxp, j);
601 if (j == -1) {
602 dres[i] = lval;
603 }
604 else if (j == lenxp) {
605 dres[i] = rval;
606 }
607 else if (j == lenxp - 1) {
608 dres[i] = dy[j];
609 }
610 else if (dx[j] == x_val) {
611 /* Avoid potential non-finite interpolation */
612 dres[i] = dy[j];
613 }
614 else {
615 const npy_double slope =
616 (slopes != NULL) ? slopes[j] :
617 (dy[j+1] - dy[j]) / (dx[j+1] - dx[j]);
618
619 /* If we get nan in one direction, try the other */
620 dres[i] = slope*(x_val - dx[j]) + dy[j];
621 if (NPY_UNLIKELY(npy_isnan(dres[i]))) {
622 dres[i] = slope*(x_val - dx[j+1]) + dy[j+1];
623 if (NPY_UNLIKELY(npy_isnan(dres[i])) && dy[j] == dy[j+1]) {
624 dres[i] = dy[j];
625 }
626 }
627 }
628 }
629
630 NPY_END_THREADS;
631 }
632
633 PyArray_free(slopes);
634 Py_DECREF(afp);
635 Py_DECREF(axp);
636 Py_DECREF(ax);
637 return PyArray_Return(af);
638
639 fail:
640 Py_XDECREF(afp);
641 Py_XDECREF(axp);
642 Py_XDECREF(ax);
643 Py_XDECREF(af);
644 return NULL;
645 }
646
647 /* As for arr_interp but for complex fp values */
648 NPY_NO_EXPORT PyObject *
arr_interp_complex(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwdict)649 arr_interp_complex(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwdict)
650 {
651
652 PyObject *fp, *xp, *x;
653 PyObject *left = NULL, *right = NULL;
654 PyArrayObject *afp = NULL, *axp = NULL, *ax = NULL, *af = NULL;
655 npy_intp i, lenx, lenxp;
656
657 const npy_double *dx, *dz;
658 const npy_cdouble *dy;
659 npy_cdouble lval, rval;
660 npy_cdouble *dres, *slopes = NULL;
661
662 static char *kwlist[] = {"x", "xp", "fp", "left", "right", NULL};
663
664 NPY_BEGIN_THREADS_DEF;
665
666 if (!PyArg_ParseTupleAndKeywords(args, kwdict, "OOO|OO:interp_complex",
667 kwlist, &x, &xp, &fp, &left, &right)) {
668 return NULL;
669 }
670
671 afp = (PyArrayObject *)PyArray_ContiguousFromAny(fp, NPY_CDOUBLE, 1, 1);
672
673 if (afp == NULL) {
674 return NULL;
675 }
676
677 axp = (PyArrayObject *)PyArray_ContiguousFromAny(xp, NPY_DOUBLE, 1, 1);
678 if (axp == NULL) {
679 goto fail;
680 }
681 ax = (PyArrayObject *)PyArray_ContiguousFromAny(x, NPY_DOUBLE, 0, 0);
682 if (ax == NULL) {
683 goto fail;
684 }
685 lenxp = PyArray_SIZE(axp);
686 if (lenxp == 0) {
687 PyErr_SetString(PyExc_ValueError,
688 "array of sample points is empty");
689 goto fail;
690 }
691 if (PyArray_SIZE(afp) != lenxp) {
692 PyErr_SetString(PyExc_ValueError,
693 "fp and xp are not of the same length.");
694 goto fail;
695 }
696
697 lenx = PyArray_SIZE(ax);
698 dx = (const npy_double *)PyArray_DATA(axp);
699 dz = (const npy_double *)PyArray_DATA(ax);
700
701 af = (PyArrayObject *)PyArray_SimpleNew(PyArray_NDIM(ax),
702 PyArray_DIMS(ax), NPY_CDOUBLE);
703 if (af == NULL) {
704 goto fail;
705 }
706
707 dy = (const npy_cdouble *)PyArray_DATA(afp);
708 dres = (npy_cdouble *)PyArray_DATA(af);
709 /* Get left and right fill values. */
710 if ((left == NULL) || (left == Py_None)) {
711 lval = dy[0];
712 }
713 else {
714 lval.real = PyComplex_RealAsDouble(left);
715 if (error_converting(lval.real)) {
716 goto fail;
717 }
718 lval.imag = PyComplex_ImagAsDouble(left);
719 if (error_converting(lval.imag)) {
720 goto fail;
721 }
722 }
723
724 if ((right == NULL) || (right == Py_None)) {
725 rval = dy[lenxp - 1];
726 }
727 else {
728 rval.real = PyComplex_RealAsDouble(right);
729 if (error_converting(rval.real)) {
730 goto fail;
731 }
732 rval.imag = PyComplex_ImagAsDouble(right);
733 if (error_converting(rval.imag)) {
734 goto fail;
735 }
736 }
737
738 /* binary_search_with_guess needs at least a 3 item long array */
739 if (lenxp == 1) {
740 const npy_double xp_val = dx[0];
741 const npy_cdouble fp_val = dy[0];
742
743 NPY_BEGIN_THREADS_THRESHOLDED(lenx);
744 for (i = 0; i < lenx; ++i) {
745 const npy_double x_val = dz[i];
746 dres[i] = (x_val < xp_val) ? lval :
747 ((x_val > xp_val) ? rval : fp_val);
748 }
749 NPY_END_THREADS;
750 }
751 else {
752 npy_intp j = 0;
753
754 /* only pre-calculate slopes if there are relatively few of them. */
755 if (lenxp <= lenx) {
756 slopes = PyArray_malloc((lenxp - 1) * sizeof(npy_cdouble));
757 if (slopes == NULL) {
758 PyErr_NoMemory();
759 goto fail;
760 }
761 }
762
763 NPY_BEGIN_THREADS;
764
765 if (slopes != NULL) {
766 for (i = 0; i < lenxp - 1; ++i) {
767 const double inv_dx = 1.0 / (dx[i+1] - dx[i]);
768 slopes[i].real = (dy[i+1].real - dy[i].real) * inv_dx;
769 slopes[i].imag = (dy[i+1].imag - dy[i].imag) * inv_dx;
770 }
771 }
772
773 for (i = 0; i < lenx; ++i) {
774 const npy_double x_val = dz[i];
775
776 if (npy_isnan(x_val)) {
777 dres[i].real = x_val;
778 dres[i].imag = 0.0;
779 continue;
780 }
781
782 j = binary_search_with_guess(x_val, dx, lenxp, j);
783 if (j == -1) {
784 dres[i] = lval;
785 }
786 else if (j == lenxp) {
787 dres[i] = rval;
788 }
789 else if (j == lenxp - 1) {
790 dres[i] = dy[j];
791 }
792 else if (dx[j] == x_val) {
793 /* Avoid potential non-finite interpolation */
794 dres[i] = dy[j];
795 }
796 else {
797 npy_cdouble slope;
798 if (slopes != NULL) {
799 slope = slopes[j];
800 }
801 else {
802 const npy_double inv_dx = 1.0 / (dx[j+1] - dx[j]);
803 slope.real = (dy[j+1].real - dy[j].real) * inv_dx;
804 slope.imag = (dy[j+1].imag - dy[j].imag) * inv_dx;
805 }
806
807 /* If we get nan in one direction, try the other */
808 dres[i].real = slope.real*(x_val - dx[j]) + dy[j].real;
809 if (NPY_UNLIKELY(npy_isnan(dres[i].real))) {
810 dres[i].real = slope.real*(x_val - dx[j+1]) + dy[j+1].real;
811 if (NPY_UNLIKELY(npy_isnan(dres[i].real)) &&
812 dy[j].real == dy[j+1].real) {
813 dres[i].real = dy[j].real;
814 }
815 }
816 dres[i].imag = slope.imag*(x_val - dx[j]) + dy[j].imag;
817 if (NPY_UNLIKELY(npy_isnan(dres[i].imag))) {
818 dres[i].imag = slope.imag*(x_val - dx[j+1]) + dy[j+1].imag;
819 if (NPY_UNLIKELY(npy_isnan(dres[i].imag)) &&
820 dy[j].imag == dy[j+1].imag) {
821 dres[i].imag = dy[j].imag;
822 }
823 }
824 }
825 }
826
827 NPY_END_THREADS;
828 }
829 PyArray_free(slopes);
830
831 Py_DECREF(afp);
832 Py_DECREF(axp);
833 Py_DECREF(ax);
834 return PyArray_Return(af);
835
836 fail:
837 Py_XDECREF(afp);
838 Py_XDECREF(axp);
839 Py_XDECREF(ax);
840 Py_XDECREF(af);
841 return NULL;
842 }
843
844 static const char *EMPTY_SEQUENCE_ERR_MSG = "indices must be integral: the provided " \
845 "empty sequence was inferred as float. Wrap it with " \
846 "'np.array(indices, dtype=np.intp)'";
847
848 static const char *NON_INTEGRAL_ERROR_MSG = "only int indices permitted";
849
850 /* Convert obj to an ndarray with integer dtype or fail */
851 static PyArrayObject *
astype_anyint(PyObject * obj)852 astype_anyint(PyObject *obj) {
853 PyArrayObject *ret;
854
855 if (!PyArray_Check(obj)) {
856 /* prefer int dtype */
857 PyArray_Descr *dtype_guess = NULL;
858 if (PyArray_DTypeFromObject(obj, NPY_MAXDIMS, &dtype_guess) < 0) {
859 return NULL;
860 }
861 if (dtype_guess == NULL) {
862 if (PySequence_Check(obj) && PySequence_Size(obj) == 0) {
863 PyErr_SetString(PyExc_TypeError, EMPTY_SEQUENCE_ERR_MSG);
864 }
865 return NULL;
866 }
867 ret = (PyArrayObject*)PyArray_FromAny(obj, dtype_guess, 0, 0, 0, NULL);
868 if (ret == NULL) {
869 return NULL;
870 }
871 }
872 else {
873 ret = (PyArrayObject *)obj;
874 Py_INCREF(ret);
875 }
876
877 if (!(PyArray_ISINTEGER(ret) || PyArray_ISBOOL(ret))) {
878 /* ensure dtype is int-based */
879 PyErr_SetString(PyExc_TypeError, NON_INTEGRAL_ERROR_MSG);
880 Py_DECREF(ret);
881 return NULL;
882 }
883
884 return ret;
885 }
886
887 /*
888 * Converts a Python sequence into 'count' PyArrayObjects
889 *
890 * seq - Input Python object, usually a tuple but any sequence works.
891 * Must have integral content.
892 * paramname - The name of the parameter that produced 'seq'.
893 * count - How many arrays there should be (errors if it doesn't match).
894 * op - Where the arrays are placed.
895 */
int_sequence_to_arrays(PyObject * seq,char * paramname,int count,PyArrayObject ** op)896 static int int_sequence_to_arrays(PyObject *seq,
897 char *paramname,
898 int count,
899 PyArrayObject **op
900 )
901 {
902 int i;
903
904 if (!PySequence_Check(seq) || PySequence_Size(seq) != count) {
905 PyErr_Format(PyExc_ValueError,
906 "parameter %s must be a sequence of length %d",
907 paramname, count);
908 return -1;
909 }
910
911 for (i = 0; i < count; ++i) {
912 PyObject *item = PySequence_GetItem(seq, i);
913 if (item == NULL) {
914 goto fail;
915 }
916 op[i] = astype_anyint(item);
917 Py_DECREF(item);
918 if (op[i] == NULL) {
919 goto fail;
920 }
921 }
922
923 return 0;
924
925 fail:
926 while (--i >= 0) {
927 Py_XDECREF(op[i]);
928 op[i] = NULL;
929 }
930 return -1;
931 }
932
933 /* Inner loop for ravel_multi_index */
934 static int
ravel_multi_index_loop(int ravel_ndim,npy_intp * ravel_dims,npy_intp * ravel_strides,npy_intp count,NPY_CLIPMODE * modes,char ** coords,npy_intp * coords_strides)935 ravel_multi_index_loop(int ravel_ndim, npy_intp *ravel_dims,
936 npy_intp *ravel_strides,
937 npy_intp count,
938 NPY_CLIPMODE *modes,
939 char **coords, npy_intp *coords_strides)
940 {
941 int i;
942 char invalid;
943 npy_intp j, m;
944
945 /*
946 * Check for 0-dimensional axes unless there is nothing to do.
947 * An empty array/shape cannot be indexed at all.
948 */
949 if (count != 0) {
950 for (i = 0; i < ravel_ndim; ++i) {
951 if (ravel_dims[i] == 0) {
952 PyErr_SetString(PyExc_ValueError,
953 "cannot unravel if shape has zero entries (is empty).");
954 return NPY_FAIL;
955 }
956 }
957 }
958
959 NPY_BEGIN_ALLOW_THREADS;
960 invalid = 0;
961 while (count--) {
962 npy_intp raveled = 0;
963 for (i = 0; i < ravel_ndim; ++i) {
964 m = ravel_dims[i];
965 j = *(npy_intp *)coords[i];
966 switch (modes[i]) {
967 case NPY_RAISE:
968 if (j < 0 || j >= m) {
969 invalid = 1;
970 goto end_while;
971 }
972 break;
973 case NPY_WRAP:
974 if (j < 0) {
975 j += m;
976 if (j < 0) {
977 j = j % m;
978 if (j != 0) {
979 j += m;
980 }
981 }
982 }
983 else if (j >= m) {
984 j -= m;
985 if (j >= m) {
986 j = j % m;
987 }
988 }
989 break;
990 case NPY_CLIP:
991 if (j < 0) {
992 j = 0;
993 }
994 else if (j >= m) {
995 j = m - 1;
996 }
997 break;
998
999 }
1000 raveled += j * ravel_strides[i];
1001
1002 coords[i] += coords_strides[i];
1003 }
1004 *(npy_intp *)coords[ravel_ndim] = raveled;
1005 coords[ravel_ndim] += coords_strides[ravel_ndim];
1006 }
1007 end_while:
1008 NPY_END_ALLOW_THREADS;
1009 if (invalid) {
1010 PyErr_SetString(PyExc_ValueError,
1011 "invalid entry in coordinates array");
1012 return NPY_FAIL;
1013 }
1014 return NPY_SUCCEED;
1015 }
1016
1017 /* ravel_multi_index implementation - see add_newdocs.py */
1018 NPY_NO_EXPORT PyObject *
arr_ravel_multi_index(PyObject * self,PyObject * args,PyObject * kwds)1019 arr_ravel_multi_index(PyObject *self, PyObject *args, PyObject *kwds)
1020 {
1021 int i;
1022 PyObject *mode0=NULL, *coords0=NULL;
1023 PyArrayObject *ret = NULL;
1024 PyArray_Dims dimensions={0,0};
1025 npy_intp s, ravel_strides[NPY_MAXDIMS];
1026 NPY_ORDER order = NPY_CORDER;
1027 NPY_CLIPMODE modes[NPY_MAXDIMS];
1028
1029 PyArrayObject *op[NPY_MAXARGS];
1030 PyArray_Descr *dtype[NPY_MAXARGS];
1031 npy_uint32 op_flags[NPY_MAXARGS];
1032
1033 NpyIter *iter = NULL;
1034
1035 char *kwlist[] = {"multi_index", "dims", "mode", "order", NULL};
1036
1037 memset(op, 0, sizeof(op));
1038 dtype[0] = NULL;
1039
1040 if (!PyArg_ParseTupleAndKeywords(args, kwds,
1041 "OO&|OO&:ravel_multi_index", kwlist,
1042 &coords0,
1043 PyArray_IntpConverter, &dimensions,
1044 &mode0,
1045 PyArray_OrderConverter, &order)) {
1046 goto fail;
1047 }
1048
1049 if (dimensions.len+1 > NPY_MAXARGS) {
1050 PyErr_SetString(PyExc_ValueError,
1051 "too many dimensions passed to ravel_multi_index");
1052 goto fail;
1053 }
1054
1055 if (!PyArray_ConvertClipmodeSequence(mode0, modes, dimensions.len)) {
1056 goto fail;
1057 }
1058
1059 switch (order) {
1060 case NPY_CORDER:
1061 s = 1;
1062 for (i = dimensions.len-1; i >= 0; --i) {
1063 ravel_strides[i] = s;
1064 if (npy_mul_with_overflow_intp(&s, s, dimensions.ptr[i])) {
1065 PyErr_SetString(PyExc_ValueError,
1066 "invalid dims: array size defined by dims is larger "
1067 "than the maximum possible size.");
1068 goto fail;
1069 }
1070 }
1071 break;
1072 case NPY_FORTRANORDER:
1073 s = 1;
1074 for (i = 0; i < dimensions.len; ++i) {
1075 ravel_strides[i] = s;
1076 if (npy_mul_with_overflow_intp(&s, s, dimensions.ptr[i])) {
1077 PyErr_SetString(PyExc_ValueError,
1078 "invalid dims: array size defined by dims is larger "
1079 "than the maximum possible size.");
1080 goto fail;
1081 }
1082 }
1083 break;
1084 default:
1085 PyErr_SetString(PyExc_ValueError,
1086 "only 'C' or 'F' order is permitted");
1087 goto fail;
1088 }
1089
1090 /* Get the multi_index into op */
1091 if (int_sequence_to_arrays(coords0, "multi_index", dimensions.len, op) < 0) {
1092 goto fail;
1093 }
1094
1095 for (i = 0; i < dimensions.len; ++i) {
1096 op_flags[i] = NPY_ITER_READONLY|
1097 NPY_ITER_ALIGNED;
1098 }
1099 op_flags[dimensions.len] = NPY_ITER_WRITEONLY|
1100 NPY_ITER_ALIGNED|
1101 NPY_ITER_ALLOCATE;
1102 dtype[0] = PyArray_DescrFromType(NPY_INTP);
1103 for (i = 1; i <= dimensions.len; ++i) {
1104 dtype[i] = dtype[0];
1105 }
1106
1107 iter = NpyIter_MultiNew(dimensions.len+1, op, NPY_ITER_BUFFERED|
1108 NPY_ITER_EXTERNAL_LOOP|
1109 NPY_ITER_ZEROSIZE_OK,
1110 NPY_KEEPORDER,
1111 NPY_SAME_KIND_CASTING,
1112 op_flags, dtype);
1113 if (iter == NULL) {
1114 goto fail;
1115 }
1116
1117 if (NpyIter_GetIterSize(iter) != 0) {
1118 NpyIter_IterNextFunc *iternext;
1119 char **dataptr;
1120 npy_intp *strides;
1121 npy_intp *countptr;
1122
1123 iternext = NpyIter_GetIterNext(iter, NULL);
1124 if (iternext == NULL) {
1125 goto fail;
1126 }
1127 dataptr = NpyIter_GetDataPtrArray(iter);
1128 strides = NpyIter_GetInnerStrideArray(iter);
1129 countptr = NpyIter_GetInnerLoopSizePtr(iter);
1130
1131 do {
1132 if (ravel_multi_index_loop(dimensions.len, dimensions.ptr,
1133 ravel_strides, *countptr, modes,
1134 dataptr, strides) != NPY_SUCCEED) {
1135 goto fail;
1136 }
1137 } while(iternext(iter));
1138 }
1139
1140 ret = NpyIter_GetOperandArray(iter)[dimensions.len];
1141 Py_INCREF(ret);
1142
1143 Py_DECREF(dtype[0]);
1144 for (i = 0; i < dimensions.len; ++i) {
1145 Py_XDECREF(op[i]);
1146 }
1147 npy_free_cache_dim_obj(dimensions);
1148 NpyIter_Deallocate(iter);
1149 return PyArray_Return(ret);
1150
1151 fail:
1152 Py_XDECREF(dtype[0]);
1153 for (i = 0; i < dimensions.len; ++i) {
1154 Py_XDECREF(op[i]);
1155 }
1156 npy_free_cache_dim_obj(dimensions);
1157 NpyIter_Deallocate(iter);
1158 return NULL;
1159 }
1160
1161
1162 /*
1163 * Inner loop for unravel_index
1164 * order must be NPY_CORDER or NPY_FORTRANORDER
1165 */
1166 static int
unravel_index_loop(int unravel_ndim,npy_intp const * unravel_dims,npy_intp unravel_size,npy_intp count,char * indices,npy_intp indices_stride,npy_intp * coords,NPY_ORDER order)1167 unravel_index_loop(int unravel_ndim, npy_intp const *unravel_dims,
1168 npy_intp unravel_size, npy_intp count,
1169 char *indices, npy_intp indices_stride,
1170 npy_intp *coords, NPY_ORDER order)
1171 {
1172 int i, idx;
1173 int idx_start = (order == NPY_CORDER) ? unravel_ndim - 1: 0;
1174 int idx_step = (order == NPY_CORDER) ? -1 : 1;
1175 char invalid = 0;
1176 npy_intp val = 0;
1177
1178 NPY_BEGIN_ALLOW_THREADS;
1179 /* NPY_KEEPORDER or NPY_ANYORDER have no meaning in this setting */
1180 assert(order == NPY_CORDER || order == NPY_FORTRANORDER);
1181 while (count--) {
1182 val = *(npy_intp *)indices;
1183 if (val < 0 || val >= unravel_size) {
1184 invalid = 1;
1185 break;
1186 }
1187 idx = idx_start;
1188 for (i = 0; i < unravel_ndim; ++i) {
1189 /*
1190 * Using a local seems to enable single-divide optimization
1191 * but only if the / precedes the %
1192 */
1193 npy_intp tmp = val / unravel_dims[idx];
1194 coords[idx] = val % unravel_dims[idx];
1195 val = tmp;
1196 idx += idx_step;
1197 }
1198 coords += unravel_ndim;
1199 indices += indices_stride;
1200 }
1201 NPY_END_ALLOW_THREADS;
1202 if (invalid) {
1203 PyErr_Format(PyExc_ValueError,
1204 "index %" NPY_INTP_FMT " is out of bounds for array with size "
1205 "%" NPY_INTP_FMT,
1206 val, unravel_size
1207 );
1208 return NPY_FAIL;
1209 }
1210 return NPY_SUCCEED;
1211 }
1212
1213 /* unravel_index implementation - see add_newdocs.py */
1214 NPY_NO_EXPORT PyObject *
arr_unravel_index(PyObject * self,PyObject * args,PyObject * kwds)1215 arr_unravel_index(PyObject *self, PyObject *args, PyObject *kwds)
1216 {
1217 PyObject *indices0 = NULL;
1218 PyObject *ret_tuple = NULL;
1219 PyArrayObject *ret_arr = NULL;
1220 PyArrayObject *indices = NULL;
1221 PyArray_Descr *dtype = NULL;
1222 PyArray_Dims dimensions = {0, 0};
1223 NPY_ORDER order = NPY_CORDER;
1224 npy_intp unravel_size;
1225
1226 NpyIter *iter = NULL;
1227 int i, ret_ndim;
1228 npy_intp ret_dims[NPY_MAXDIMS], ret_strides[NPY_MAXDIMS];
1229
1230 char *kwlist[] = {"indices", "shape", "order", NULL};
1231
1232 /*
1233 * TODO: remove this in favor of warning raised in the dispatcher when
1234 * __array_function__ is enabled by default.
1235 */
1236
1237 /*
1238 * Continue to support the older "dims" argument in place
1239 * of the "shape" argument. Issue an appropriate warning
1240 * if "dims" is detected in keywords, then replace it with
1241 * the new "shape" argument and continue processing as usual.
1242 */
1243 if (kwds) {
1244 PyObject *dims_item, *shape_item;
1245 dims_item = _PyDict_GetItemStringWithError(kwds, "dims");
1246 if (dims_item == NULL && PyErr_Occurred()){
1247 return NULL;
1248 }
1249 shape_item = _PyDict_GetItemStringWithError(kwds, "shape");
1250 if (shape_item == NULL && PyErr_Occurred()){
1251 return NULL;
1252 }
1253 if (dims_item != NULL && shape_item == NULL) {
1254 if (DEPRECATE("'shape' argument should be"
1255 " used instead of 'dims'") < 0) {
1256 return NULL;
1257 }
1258 if (PyDict_SetItemString(kwds, "shape", dims_item) < 0) {
1259 return NULL;
1260 }
1261 if (PyDict_DelItemString(kwds, "dims") < 0) {
1262 return NULL;
1263 }
1264 }
1265 }
1266
1267 if (!PyArg_ParseTupleAndKeywords(args, kwds, "OO&|O&:unravel_index",
1268 kwlist,
1269 &indices0,
1270 PyArray_IntpConverter, &dimensions,
1271 PyArray_OrderConverter, &order)) {
1272 goto fail;
1273 }
1274
1275 unravel_size = PyArray_OverflowMultiplyList(dimensions.ptr, dimensions.len);
1276 if (unravel_size == -1) {
1277 PyErr_SetString(PyExc_ValueError,
1278 "dimensions are too large; arrays and shapes with "
1279 "a total size greater than 'intp' are not supported.");
1280 goto fail;
1281 }
1282
1283 indices = astype_anyint(indices0);
1284 if (indices == NULL) {
1285 goto fail;
1286 }
1287
1288 dtype = PyArray_DescrFromType(NPY_INTP);
1289 if (dtype == NULL) {
1290 goto fail;
1291 }
1292
1293 iter = NpyIter_New(indices, NPY_ITER_READONLY|
1294 NPY_ITER_ALIGNED|
1295 NPY_ITER_BUFFERED|
1296 NPY_ITER_ZEROSIZE_OK|
1297 NPY_ITER_DONT_NEGATE_STRIDES|
1298 NPY_ITER_MULTI_INDEX,
1299 NPY_KEEPORDER, NPY_SAME_KIND_CASTING,
1300 dtype);
1301 if (iter == NULL) {
1302 goto fail;
1303 }
1304
1305 /*
1306 * Create the return array with a layout compatible with the indices
1307 * and with a dimension added to the end for the multi-index
1308 */
1309 ret_ndim = PyArray_NDIM(indices) + 1;
1310 if (NpyIter_GetShape(iter, ret_dims) != NPY_SUCCEED) {
1311 goto fail;
1312 }
1313 ret_dims[ret_ndim-1] = dimensions.len;
1314 if (NpyIter_CreateCompatibleStrides(iter,
1315 dimensions.len*sizeof(npy_intp), ret_strides) != NPY_SUCCEED) {
1316 goto fail;
1317 }
1318 ret_strides[ret_ndim-1] = sizeof(npy_intp);
1319
1320 /* Remove the multi-index and inner loop */
1321 if (NpyIter_RemoveMultiIndex(iter) != NPY_SUCCEED) {
1322 goto fail;
1323 }
1324 if (NpyIter_EnableExternalLoop(iter) != NPY_SUCCEED) {
1325 goto fail;
1326 }
1327
1328 ret_arr = (PyArrayObject *)PyArray_NewFromDescr(&PyArray_Type, dtype,
1329 ret_ndim, ret_dims, ret_strides, NULL, 0, NULL);
1330 dtype = NULL;
1331 if (ret_arr == NULL) {
1332 goto fail;
1333 }
1334
1335 if (order != NPY_CORDER && order != NPY_FORTRANORDER) {
1336 PyErr_SetString(PyExc_ValueError,
1337 "only 'C' or 'F' order is permitted");
1338 goto fail;
1339 }
1340 if (NpyIter_GetIterSize(iter) != 0) {
1341 NpyIter_IterNextFunc *iternext;
1342 char **dataptr;
1343 npy_intp *strides;
1344 npy_intp *countptr, count;
1345 npy_intp *coordsptr = (npy_intp *)PyArray_DATA(ret_arr);
1346
1347 iternext = NpyIter_GetIterNext(iter, NULL);
1348 if (iternext == NULL) {
1349 goto fail;
1350 }
1351 dataptr = NpyIter_GetDataPtrArray(iter);
1352 strides = NpyIter_GetInnerStrideArray(iter);
1353 countptr = NpyIter_GetInnerLoopSizePtr(iter);
1354
1355 do {
1356 count = *countptr;
1357 if (unravel_index_loop(dimensions.len, dimensions.ptr,
1358 unravel_size, count, *dataptr, *strides,
1359 coordsptr, order) != NPY_SUCCEED) {
1360 goto fail;
1361 }
1362 coordsptr += count * dimensions.len;
1363 } while (iternext(iter));
1364 }
1365
1366
1367 if (dimensions.len == 0 && PyArray_NDIM(indices) != 0) {
1368 /*
1369 * There's no index meaning "take the only element 10 times"
1370 * on a zero-d array, so we have no choice but to error. (See gh-580)
1371 *
1372 * Do this check after iterating, so we give a better error message
1373 * for invalid indices.
1374 */
1375 PyErr_SetString(PyExc_ValueError,
1376 "multiple indices are not supported for 0d arrays");
1377 goto fail;
1378 }
1379
1380 /* Now make a tuple of views, one per index */
1381 ret_tuple = PyTuple_New(dimensions.len);
1382 if (ret_tuple == NULL) {
1383 goto fail;
1384 }
1385 for (i = 0; i < dimensions.len; ++i) {
1386 PyArrayObject *view;
1387
1388 view = (PyArrayObject *)PyArray_NewFromDescrAndBase(
1389 &PyArray_Type, PyArray_DescrFromType(NPY_INTP),
1390 ret_ndim - 1, ret_dims, ret_strides,
1391 PyArray_BYTES(ret_arr) + i*sizeof(npy_intp),
1392 NPY_ARRAY_WRITEABLE, NULL, (PyObject *)ret_arr);
1393 if (view == NULL) {
1394 goto fail;
1395 }
1396 PyTuple_SET_ITEM(ret_tuple, i, PyArray_Return(view));
1397 }
1398
1399 Py_DECREF(ret_arr);
1400 Py_XDECREF(indices);
1401 npy_free_cache_dim_obj(dimensions);
1402 NpyIter_Deallocate(iter);
1403
1404 return ret_tuple;
1405
1406 fail:
1407 Py_XDECREF(ret_tuple);
1408 Py_XDECREF(ret_arr);
1409 Py_XDECREF(dtype);
1410 Py_XDECREF(indices);
1411 npy_free_cache_dim_obj(dimensions);
1412 NpyIter_Deallocate(iter);
1413 return NULL;
1414 }
1415
1416
1417 /* Can only be called if doc is currently NULL */
1418 NPY_NO_EXPORT PyObject *
arr_add_docstring(PyObject * NPY_UNUSED (dummy),PyObject * args)1419 arr_add_docstring(PyObject *NPY_UNUSED(dummy), PyObject *args)
1420 {
1421 PyObject *obj;
1422 PyObject *str;
1423 #if PY_VERSION_HEX >= 0x030700A2 && (!defined(PYPY_VERSION_NUM) || PYPY_VERSION_NUM > 0x07030300)
1424 const char *docstr;
1425 #else
1426 char *docstr;
1427 #endif
1428 static char *msg = "already has a different docstring";
1429
1430 /* Don't add docstrings */
1431 if (Py_OptimizeFlag > 1) {
1432 Py_RETURN_NONE;
1433 }
1434
1435 if (!PyArg_ParseTuple(args, "OO!:add_docstring", &obj, &PyUnicode_Type, &str)) {
1436 return NULL;
1437 }
1438
1439 docstr = PyUnicode_AsUTF8(str);
1440 if (docstr == NULL) {
1441 return NULL;
1442 }
1443
1444 #define _ADDDOC(doc, name) \
1445 if (!(doc)) { \
1446 doc = docstr; \
1447 Py_INCREF(str); /* hold on to string (leaks reference) */ \
1448 } \
1449 else if (strcmp(doc, docstr) != 0) { \
1450 PyErr_Format(PyExc_RuntimeError, "%s method %s", name, msg); \
1451 return NULL; \
1452 }
1453
1454 if (Py_TYPE(obj) == &PyCFunction_Type) {
1455 PyCFunctionObject *new = (PyCFunctionObject *)obj;
1456 _ADDDOC(new->m_ml->ml_doc, new->m_ml->ml_name);
1457 }
1458 else if (Py_TYPE(obj) == &PyType_Type) {
1459 PyTypeObject *new = (PyTypeObject *)obj;
1460 _ADDDOC(new->tp_doc, new->tp_name);
1461 }
1462 else if (Py_TYPE(obj) == &PyMemberDescr_Type) {
1463 PyMemberDescrObject *new = (PyMemberDescrObject *)obj;
1464 _ADDDOC(new->d_member->doc, new->d_member->name);
1465 }
1466 else if (Py_TYPE(obj) == &PyGetSetDescr_Type) {
1467 PyGetSetDescrObject *new = (PyGetSetDescrObject *)obj;
1468 _ADDDOC(new->d_getset->doc, new->d_getset->name);
1469 }
1470 else if (Py_TYPE(obj) == &PyMethodDescr_Type) {
1471 PyMethodDescrObject *new = (PyMethodDescrObject *)obj;
1472 _ADDDOC(new->d_method->ml_doc, new->d_method->ml_name);
1473 }
1474 else {
1475 PyObject *doc_attr;
1476
1477 doc_attr = PyObject_GetAttrString(obj, "__doc__");
1478 if (doc_attr != NULL && doc_attr != Py_None &&
1479 (PyUnicode_Compare(doc_attr, str) != 0)) {
1480 Py_DECREF(doc_attr);
1481 if (PyErr_Occurred()) {
1482 /* error during PyUnicode_Compare */
1483 return NULL;
1484 }
1485 PyErr_Format(PyExc_RuntimeError, "object %s", msg);
1486 return NULL;
1487 }
1488 Py_XDECREF(doc_attr);
1489
1490 if (PyObject_SetAttrString(obj, "__doc__", str) < 0) {
1491 PyErr_SetString(PyExc_TypeError,
1492 "Cannot set a docstring for that object");
1493 return NULL;
1494 }
1495 Py_RETURN_NONE;
1496 }
1497
1498 #undef _ADDDOC
1499
1500 Py_RETURN_NONE;
1501 }
1502
1503 #if defined NPY_HAVE_SSE2_INTRINSICS
1504 #include <emmintrin.h>
1505 #endif
1506
1507 #ifdef NPY_HAVE_NEON
1508 typedef npy_uint64 uint64_unaligned __attribute__((aligned(16)));
1509 static NPY_INLINE int32_t
sign_mask(uint8x16_t input)1510 sign_mask(uint8x16_t input)
1511 {
1512 int8x8_t m0 = vcreate_s8(0x0706050403020100ULL);
1513 uint8x16_t v0 = vshlq_u8(vshrq_n_u8(input, 7), vcombine_s8(m0, m0));
1514 uint64x2_t v1 = vpaddlq_u32(vpaddlq_u16(vpaddlq_u8(v0)));
1515 return (int)vgetq_lane_u64(v1, 0) + ((int)vgetq_lane_u64(v1, 1) << 8);
1516 }
1517 #endif
1518 /*
1519 * This function packs boolean values in the input array into the bits of a
1520 * byte array. Truth values are determined as usual: 0 is false, everything
1521 * else is true.
1522 */
1523 static NPY_INLINE void
pack_inner(const char * inptr,npy_intp element_size,npy_intp n_in,npy_intp in_stride,char * outptr,npy_intp n_out,npy_intp out_stride,char order)1524 pack_inner(const char *inptr,
1525 npy_intp element_size, /* in bytes */
1526 npy_intp n_in,
1527 npy_intp in_stride,
1528 char *outptr,
1529 npy_intp n_out,
1530 npy_intp out_stride,
1531 char order)
1532 {
1533 /*
1534 * Loop through the elements of inptr.
1535 * Determine whether or not it is nonzero.
1536 * Yes: set corresponding bit (and adjust build value)
1537 * No: move on
1538 * Every 8th value, set the value of build and increment the outptr
1539 */
1540 npy_intp index = 0;
1541 int remain = n_in % 8; /* uneven bits */
1542
1543 #if defined NPY_HAVE_SSE2_INTRINSICS && defined HAVE__M_FROM_INT64
1544 if (in_stride == 1 && element_size == 1 && n_out > 2) {
1545 __m128i zero = _mm_setzero_si128();
1546 /* don't handle non-full 8-byte remainder */
1547 npy_intp vn_out = n_out - (remain ? 1 : 0);
1548 vn_out -= (vn_out & 1);
1549 for (index = 0; index < vn_out; index += 2) {
1550 unsigned int r;
1551 npy_uint64 a = *(npy_uint64*)inptr;
1552 npy_uint64 b = *(npy_uint64*)(inptr + 8);
1553 if (order == 'b') {
1554 a = npy_bswap8(a);
1555 b = npy_bswap8(b);
1556 }
1557
1558 /* note x86 can load unaligned */
1559 __m128i v = _mm_set_epi64(_m_from_int64(b), _m_from_int64(a));
1560 /* false -> 0x00 and true -> 0xFF (there is no cmpneq) */
1561 v = _mm_cmpeq_epi8(v, zero);
1562 v = _mm_cmpeq_epi8(v, zero);
1563 /* extract msb of 16 bytes and pack it into 16 bit */
1564 r = _mm_movemask_epi8(v);
1565 /* store result */
1566 memcpy(outptr, &r, 1);
1567 outptr += out_stride;
1568 memcpy(outptr, (char*)&r + 1, 1);
1569 outptr += out_stride;
1570 inptr += 16;
1571 }
1572 }
1573 #elif defined NPY_HAVE_NEON
1574 if (in_stride == 1 && element_size == 1 && n_out > 2) {
1575 /* don't handle non-full 8-byte remainder */
1576 npy_intp vn_out = n_out - (remain ? 1 : 0);
1577 vn_out -= (vn_out & 1);
1578 for (index = 0; index < vn_out; index += 2) {
1579 unsigned int r;
1580 npy_uint64 a = *((uint64_unaligned*)inptr);
1581 npy_uint64 b = *((uint64_unaligned*)(inptr + 8));
1582 if (order == 'b') {
1583 a = npy_bswap8(a);
1584 b = npy_bswap8(b);
1585 }
1586 uint64x2_t v = vcombine_u64(vcreate_u64(a), vcreate_u64(b));
1587 uint64x2_t zero = vdupq_n_u64(0);
1588 /* false -> 0x00 and true -> 0xFF */
1589 v = vreinterpretq_u64_u8(vmvnq_u8(vceqq_u8(vreinterpretq_u8_u64(v), vreinterpretq_u8_u64(zero))));
1590 /* extract msb of 16 bytes and pack it into 16 bit */
1591 uint8x16_t input = vreinterpretq_u8_u64(v);
1592 r = sign_mask(input);
1593 /* store result */
1594 memcpy(outptr, &r, 1);
1595 outptr += out_stride;
1596 memcpy(outptr, (char*)&r + 1, 1);
1597 outptr += out_stride;
1598 inptr += 16;
1599 }
1600 }
1601 #endif
1602
1603 if (remain == 0) { /* assumes n_in > 0 */
1604 remain = 8;
1605 }
1606 /* Don't reset index. Just handle remainder of above block */
1607 for (; index < n_out; index++) {
1608 unsigned char build = 0;
1609 int i, maxi;
1610 npy_intp j;
1611
1612 maxi = (index == n_out - 1) ? remain : 8;
1613 if (order == 'b') {
1614 for (i = 0; i < maxi; i++) {
1615 build <<= 1;
1616 for (j = 0; j < element_size; j++) {
1617 build |= (inptr[j] != 0);
1618 }
1619 inptr += in_stride;
1620 }
1621 if (index == n_out - 1) {
1622 build <<= 8 - remain;
1623 }
1624 }
1625 else
1626 {
1627 for (i = 0; i < maxi; i++) {
1628 build >>= 1;
1629 for (j = 0; j < element_size; j++) {
1630 build |= (inptr[j] != 0) ? 128 : 0;
1631 }
1632 inptr += in_stride;
1633 }
1634 if (index == n_out - 1) {
1635 build >>= 8 - remain;
1636 }
1637 }
1638 *outptr = (char)build;
1639 outptr += out_stride;
1640 }
1641 }
1642
1643 static PyObject *
pack_bits(PyObject * input,int axis,char order)1644 pack_bits(PyObject *input, int axis, char order)
1645 {
1646 PyArrayObject *inp;
1647 PyArrayObject *new = NULL;
1648 PyArrayObject *out = NULL;
1649 npy_intp outdims[NPY_MAXDIMS];
1650 int i;
1651 PyArrayIterObject *it, *ot;
1652 NPY_BEGIN_THREADS_DEF;
1653
1654 inp = (PyArrayObject *)PyArray_FROM_O(input);
1655
1656 if (inp == NULL) {
1657 return NULL;
1658 }
1659 if (!PyArray_ISBOOL(inp) && !PyArray_ISINTEGER(inp)) {
1660 PyErr_SetString(PyExc_TypeError,
1661 "Expected an input array of integer or boolean data type");
1662 Py_DECREF(inp);
1663 goto fail;
1664 }
1665
1666 new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0);
1667 Py_DECREF(inp);
1668 if (new == NULL) {
1669 return NULL;
1670 }
1671
1672 if (PyArray_NDIM(new) == 0) {
1673 char *optr, *iptr;
1674
1675 out = (PyArrayObject *)PyArray_NewFromDescr(
1676 Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1677 0, NULL, NULL, NULL,
1678 0, NULL);
1679 if (out == NULL) {
1680 goto fail;
1681 }
1682 optr = PyArray_DATA(out);
1683 iptr = PyArray_DATA(new);
1684 *optr = 0;
1685 for (i = 0; i < PyArray_ITEMSIZE(new); i++) {
1686 if (*iptr != 0) {
1687 *optr = 1;
1688 break;
1689 }
1690 iptr++;
1691 }
1692 goto finish;
1693 }
1694
1695
1696 /* Setup output shape */
1697 for (i = 0; i < PyArray_NDIM(new); i++) {
1698 outdims[i] = PyArray_DIM(new, i);
1699 }
1700
1701 /*
1702 * Divide axis dimension by 8
1703 * 8 -> 1, 9 -> 2, 16 -> 2, 17 -> 3 etc..
1704 */
1705 outdims[axis] = ((outdims[axis] - 1) >> 3) + 1;
1706
1707 /* Create output array */
1708 out = (PyArrayObject *)PyArray_NewFromDescr(
1709 Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1710 PyArray_NDIM(new), outdims, NULL, NULL,
1711 PyArray_ISFORTRAN(new), NULL);
1712 if (out == NULL) {
1713 goto fail;
1714 }
1715 /* Setup iterators to iterate over all but given axis */
1716 it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis);
1717 ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis);
1718 if (it == NULL || ot == NULL) {
1719 Py_XDECREF(it);
1720 Py_XDECREF(ot);
1721 goto fail;
1722 }
1723
1724 NPY_BEGIN_THREADS_THRESHOLDED(PyArray_DIM(out, axis));
1725 while (PyArray_ITER_NOTDONE(it)) {
1726 pack_inner(PyArray_ITER_DATA(it), PyArray_ITEMSIZE(new),
1727 PyArray_DIM(new, axis), PyArray_STRIDE(new, axis),
1728 PyArray_ITER_DATA(ot), PyArray_DIM(out, axis),
1729 PyArray_STRIDE(out, axis), order);
1730 PyArray_ITER_NEXT(it);
1731 PyArray_ITER_NEXT(ot);
1732 }
1733 NPY_END_THREADS;
1734
1735 Py_DECREF(it);
1736 Py_DECREF(ot);
1737
1738 finish:
1739 Py_DECREF(new);
1740 return (PyObject *)out;
1741
1742 fail:
1743 Py_XDECREF(new);
1744 Py_XDECREF(out);
1745 return NULL;
1746 }
1747
1748 static PyObject *
unpack_bits(PyObject * input,int axis,PyObject * count_obj,char order)1749 unpack_bits(PyObject *input, int axis, PyObject *count_obj, char order)
1750 {
1751 static int unpack_init = 0;
1752 /*
1753 * lookuptable for bitorder big as it has been around longer
1754 * bitorder little is handled via byteswapping in the loop
1755 */
1756 static union {
1757 npy_uint8 bytes[8];
1758 npy_uint64 uint64;
1759 } unpack_lookup_big[256];
1760 PyArrayObject *inp;
1761 PyArrayObject *new = NULL;
1762 PyArrayObject *out = NULL;
1763 npy_intp outdims[NPY_MAXDIMS];
1764 int i;
1765 PyArrayIterObject *it, *ot;
1766 npy_intp count, in_n, in_tail, out_pad, in_stride, out_stride;
1767 NPY_BEGIN_THREADS_DEF;
1768
1769 inp = (PyArrayObject *)PyArray_FROM_O(input);
1770
1771 if (inp == NULL) {
1772 return NULL;
1773 }
1774 if (PyArray_TYPE(inp) != NPY_UBYTE) {
1775 PyErr_SetString(PyExc_TypeError,
1776 "Expected an input array of unsigned byte data type");
1777 Py_DECREF(inp);
1778 goto fail;
1779 }
1780
1781 new = (PyArrayObject *)PyArray_CheckAxis(inp, &axis, 0);
1782 Py_DECREF(inp);
1783 if (new == NULL) {
1784 return NULL;
1785 }
1786
1787 if (PyArray_NDIM(new) == 0) {
1788 /* Handle 0-d array by converting it to a 1-d array */
1789 PyArrayObject *temp;
1790 PyArray_Dims newdim = {NULL, 1};
1791 npy_intp shape = 1;
1792
1793 newdim.ptr = &shape;
1794 temp = (PyArrayObject *)PyArray_Newshape(new, &newdim, NPY_CORDER);
1795 Py_DECREF(new);
1796 if (temp == NULL) {
1797 return NULL;
1798 }
1799 new = temp;
1800 }
1801
1802 /* Setup output shape */
1803 for (i = 0; i < PyArray_NDIM(new); i++) {
1804 outdims[i] = PyArray_DIM(new, i);
1805 }
1806
1807 /* Multiply axis dimension by 8 */
1808 outdims[axis] *= 8;
1809 if (count_obj != Py_None) {
1810 count = PyArray_PyIntAsIntp(count_obj);
1811 if (error_converting(count)) {
1812 goto fail;
1813 }
1814 if (count < 0) {
1815 outdims[axis] += count;
1816 if (outdims[axis] < 0) {
1817 PyErr_Format(PyExc_ValueError,
1818 "-count larger than number of elements");
1819 goto fail;
1820 }
1821 }
1822 else {
1823 outdims[axis] = count;
1824 }
1825 }
1826
1827 /* Create output array */
1828 out = (PyArrayObject *)PyArray_NewFromDescr(
1829 Py_TYPE(new), PyArray_DescrFromType(NPY_UBYTE),
1830 PyArray_NDIM(new), outdims, NULL, NULL,
1831 PyArray_ISFORTRAN(new), NULL);
1832 if (out == NULL) {
1833 goto fail;
1834 }
1835
1836 /* Setup iterators to iterate over all but given axis */
1837 it = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)new, &axis);
1838 ot = (PyArrayIterObject *)PyArray_IterAllButAxis((PyObject *)out, &axis);
1839 if (it == NULL || ot == NULL) {
1840 Py_XDECREF(it);
1841 Py_XDECREF(ot);
1842 goto fail;
1843 }
1844
1845 /*
1846 * setup lookup table under GIL, 256 8 byte blocks representing 8 bits
1847 * expanded to 1/0 bytes
1848 */
1849 if (unpack_init == 0) {
1850 npy_intp j;
1851 for (j=0; j < 256; j++) {
1852 npy_intp k;
1853 for (k=0; k < 8; k++) {
1854 npy_uint8 v = (j & (1 << k)) == (1 << k);
1855 unpack_lookup_big[j].bytes[7 - k] = v;
1856 }
1857 }
1858 unpack_init = 1;
1859 }
1860
1861 count = PyArray_DIM(new, axis) * 8;
1862 if (outdims[axis] > count) {
1863 in_n = count / 8;
1864 in_tail = 0;
1865 out_pad = outdims[axis] - count;
1866 }
1867 else {
1868 in_n = outdims[axis] / 8;
1869 in_tail = outdims[axis] % 8;
1870 out_pad = 0;
1871 }
1872
1873 in_stride = PyArray_STRIDE(new, axis);
1874 out_stride = PyArray_STRIDE(out, axis);
1875
1876 NPY_BEGIN_THREADS_THRESHOLDED(PyArray_Size((PyObject *)out) / 8);
1877
1878 while (PyArray_ITER_NOTDONE(it)) {
1879 npy_intp index;
1880 unsigned const char *inptr = PyArray_ITER_DATA(it);
1881 char *outptr = PyArray_ITER_DATA(ot);
1882
1883 if (out_stride == 1) {
1884 /* for unity stride we can just copy out of the lookup table */
1885 if (order == 'b') {
1886 for (index = 0; index < in_n; index++) {
1887 npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1888 memcpy(outptr, &v, 8);
1889 outptr += 8;
1890 inptr += in_stride;
1891 }
1892 }
1893 else {
1894 for (index = 0; index < in_n; index++) {
1895 npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1896 if (order != 'b') {
1897 v = npy_bswap8(v);
1898 }
1899 memcpy(outptr, &v, 8);
1900 outptr += 8;
1901 inptr += in_stride;
1902 }
1903 }
1904 /* Clean up the tail portion */
1905 if (in_tail) {
1906 npy_uint64 v = unpack_lookup_big[*inptr].uint64;
1907 if (order != 'b') {
1908 v = npy_bswap8(v);
1909 }
1910 memcpy(outptr, &v, in_tail);
1911 }
1912 /* Add padding */
1913 else if (out_pad) {
1914 memset(outptr, 0, out_pad);
1915 }
1916 }
1917 else {
1918 if (order == 'b') {
1919 for (index = 0; index < in_n; index++) {
1920 for (i = 0; i < 8; i++) {
1921 *outptr = ((*inptr & (128 >> i)) != 0);
1922 outptr += out_stride;
1923 }
1924 inptr += in_stride;
1925 }
1926 /* Clean up the tail portion */
1927 for (i = 0; i < in_tail; i++) {
1928 *outptr = ((*inptr & (128 >> i)) != 0);
1929 outptr += out_stride;
1930 }
1931 }
1932 else {
1933 for (index = 0; index < in_n; index++) {
1934 for (i = 0; i < 8; i++) {
1935 *outptr = ((*inptr & (1 << i)) != 0);
1936 outptr += out_stride;
1937 }
1938 inptr += in_stride;
1939 }
1940 /* Clean up the tail portion */
1941 for (i = 0; i < in_tail; i++) {
1942 *outptr = ((*inptr & (1 << i)) != 0);
1943 outptr += out_stride;
1944 }
1945 }
1946 /* Add padding */
1947 for (index = 0; index < out_pad; index++) {
1948 *outptr = 0;
1949 outptr += out_stride;
1950 }
1951 }
1952
1953 PyArray_ITER_NEXT(it);
1954 PyArray_ITER_NEXT(ot);
1955 }
1956 NPY_END_THREADS;
1957
1958 Py_DECREF(it);
1959 Py_DECREF(ot);
1960
1961 Py_DECREF(new);
1962 return (PyObject *)out;
1963
1964 fail:
1965 Py_XDECREF(new);
1966 Py_XDECREF(out);
1967 return NULL;
1968 }
1969
1970
1971 NPY_NO_EXPORT PyObject *
io_pack(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwds)1972 io_pack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
1973 {
1974 PyObject *obj;
1975 int axis = NPY_MAXDIMS;
1976 static char *kwlist[] = {"in", "axis", "bitorder", NULL};
1977 char c = 'b';
1978 const char * order_str = NULL;
1979
1980 if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&s:pack" , kwlist,
1981 &obj, PyArray_AxisConverter, &axis, &order_str)) {
1982 return NULL;
1983 }
1984 if (order_str != NULL) {
1985 if (strncmp(order_str, "little", 6) == 0)
1986 c = 'l';
1987 else if (strncmp(order_str, "big", 3) == 0)
1988 c = 'b';
1989 else {
1990 PyErr_SetString(PyExc_ValueError,
1991 "'order' must be either 'little' or 'big'");
1992 return NULL;
1993 }
1994 }
1995 return pack_bits(obj, axis, c);
1996 }
1997
1998
1999 NPY_NO_EXPORT PyObject *
io_unpack(PyObject * NPY_UNUSED (self),PyObject * args,PyObject * kwds)2000 io_unpack(PyObject *NPY_UNUSED(self), PyObject *args, PyObject *kwds)
2001 {
2002 PyObject *obj;
2003 int axis = NPY_MAXDIMS;
2004 PyObject *count = Py_None;
2005 static char *kwlist[] = {"in", "axis", "count", "bitorder", NULL};
2006 const char * c = NULL;
2007
2008 if (!PyArg_ParseTupleAndKeywords( args, kwds, "O|O&Os:unpack" , kwlist,
2009 &obj, PyArray_AxisConverter, &axis, &count, &c)) {
2010 return NULL;
2011 }
2012 if (c == NULL) {
2013 c = "b";
2014 }
2015 if (c[0] != 'l' && c[0] != 'b') {
2016 PyErr_SetString(PyExc_ValueError,
2017 "'order' must begin with 'l' or 'b'");
2018 return NULL;
2019 }
2020 return unpack_bits(obj, axis, count, c[0]);
2021 }
2022