1
/*
2
 * This file implements generic methods for computing reductions on arrays.
3
 *
4
 * Written by Mark Wiebe (mwwiebe@gmail.com)
5
 * Copyright (c) 2011 by Enthought, Inc.
6
 *
7
 * See LICENSE.txt for the license.
8
 */
9
#define _UMATHMODULE
10
#define _MULTIARRAYMODULE
11
#define NPY_NO_DEPRECATED_API NPY_API_VERSION
12

13
#define PY_SSIZE_T_CLEAN
14
#include <Python.h>
15

16
#include "npy_config.h"
17
#include <numpy/arrayobject.h>
18

19
#include "npy_pycompat.h"
20
#include "ctors.h"
21

22
#include "numpy/ufuncobject.h"
23
#include "lowlevel_strided_loops.h"
24
#include "reduction.h"
25
#include "extobj.h"  /* for _check_ufunc_fperr */
26

27

28
/*
29
 * Count the number of dimensions selected in 'axis_flags'
30
 */
31
static int
32
count_axes(int ndim, const npy_bool *axis_flags)
33
{
34
    int idim;
35
    int naxes = 0;
36

37 1
    for (idim = 0; idim < ndim; ++idim) {
38 1
        if (axis_flags[idim]) {
39 1
            naxes++;
40
        }
41
    }
42
    return naxes;
43
}
44

45
/*
46
 * This function initializes a result array for a reduction operation
47
 * which has no identity. This means it needs to copy the first element
48
 * it sees along the reduction axes to result.
49
 *
50
 * If a reduction has an identity, such as 0 or 1, the result should be
51
 * fully initialized to the identity, because this function raises an
52
 * exception when there are no elements to reduce (which is appropriate if,
53
 * and only if, the reduction operation has no identity).
54
 *
55
 * This means it copies the subarray indexed at zero along each reduction axis
56
 * into 'result'.
57
 *
58
 * result  : The array into which the result is computed. This must have
59
 *           the same number of dimensions as 'operand', but for each
60
 *           axis i where 'axis_flags[i]' is True, it has a single element.
61
 * operand : The array being reduced.
62
 * axis_flags : An array of boolean flags, one for each axis of 'operand'.
63
 *              When a flag is True, it indicates to reduce along that axis.
64
 * funcname : The name of the reduction operation, for the purpose of
65
 *            better quality error messages. For example, "numpy.max"
66
 *            would be a good name for NumPy's max function.
67
 *
68
 * Returns -1 if an error occurred, and otherwise the reduce arrays size,
69
 * which is the number of elements already initialized.
70
 */
71
NPY_NO_EXPORT int
72 1
PyArray_CopyInitialReduceValues(
73
                    PyArrayObject *result, PyArrayObject *operand,
74
                    const npy_bool *axis_flags, const char *funcname,
75
                    int keepdims)
76
{
77
    npy_intp shape[NPY_MAXDIMS], strides[NPY_MAXDIMS];
78 1
    npy_intp *shape_orig = PyArray_SHAPE(operand);
79 1
    npy_intp *strides_orig = PyArray_STRIDES(operand);
80 1
    PyArrayObject *op_view = NULL;
81

82 1
    int ndim = PyArray_NDIM(operand);
83

84
    /*
85
     * Copy the subarray of the first element along each reduction axis.
86
     *
87
     * Adjust the shape to only look at the first element along
88
     * any of the reduction axes. If keepdims is False remove the axes
89
     * entirely.
90
     */
91 1
    int idim_out = 0;
92 1
    npy_intp size = 1;
93 1
    for (int idim = 0; idim < ndim; idim++) {
94 1
        if (axis_flags[idim]) {
95 1
            if (NPY_UNLIKELY(shape_orig[idim] == 0)) {
96 1
                PyErr_Format(PyExc_ValueError,
97
                        "zero-size array to reduction operation %s "
98
                        "which has no identity", funcname);
99 1
                return -1;
100
            }
101 1
            if (keepdims) {
102 1
                shape[idim_out] = 1;
103 1
                strides[idim_out] = 0;
104 1
                idim_out++;
105
            }
106
        }
107
        else {
108 1
            size *= shape_orig[idim];
109 1
            shape[idim_out] = shape_orig[idim];
110 1
            strides[idim_out] = strides_orig[idim];
111 1
            idim_out++;
112
        }
113
    }
114

115 1
    PyArray_Descr *descr = PyArray_DESCR(operand);
116 1
    Py_INCREF(descr);
117 1
    op_view = (PyArrayObject *)PyArray_NewFromDescr(
118
            &PyArray_Type, descr, idim_out, shape, strides,
119
            PyArray_DATA(operand), 0, NULL);
120 1
    if (op_view == NULL) {
121
        return -1;
122
    }
123

124
    /*
125
     * Copy the elements into the result to start.
126
     */
127 1
    int res = PyArray_CopyInto(result, op_view);
128 1
    Py_DECREF(op_view);
129 1
    if (res < 0) {
130
        return -1;
131
    }
132

133
    /*
134
     * If there were no reduction axes, we would already be done here.
135
     * Note that if there is only a single reduction axis, in principle the
136
     * iteration could be set up more efficiently here by removing that
137
     * axis before setting up the iterator (simplifying the iteration since
138
     * `skip_first_count` (the returned size) can be set to 0).
139
     */
140 1
    return size;
141
}
142

143
/*
144
 * This function executes all the standard NumPy reduction function
145
 * boilerplate code, just calling the appropriate inner loop function where
146
 * necessary.
147
 *
148
 * operand     : The array to be reduced.
149
 * out         : NULL, or the array into which to place the result.
150
 * wheremask   : NOT YET SUPPORTED, but this parameter is placed here
151
 *               so that support can be added in the future without breaking
152
 *               API compatibility. Pass in NULL.
153
 * operand_dtype : The dtype the inner loop expects for the operand.
154
 * result_dtype : The dtype the inner loop expects for the result.
155
 * casting     : The casting rule to apply to the operands.
156
 * axis_flags  : Flags indicating the reduction axes of 'operand'.
157
 * reorderable : If True, the reduction being done is reorderable, which
158
 *               means specifying multiple axes of reduction at once is ok,
159
 *               and the reduction code may calculate the reduction in an
160
 *               arbitrary order. The calculation may be reordered because
161
 *               of cache behavior or multithreading requirements.
162
 * keepdims    : If true, leaves the reduction dimensions in the result
163
 *               with size one.
164
 * subok       : If true, the result uses the subclass of operand, otherwise
165
 *               it is always a base class ndarray.
166
 * identity    : If Py_None, PyArray_CopyInitialReduceValues is used, otherwise
167
 *               this value is used to initialize the result to
168
 *               the reduction's unit.
169
 * loop        : The loop which does the reduction.
170
 * data        : Data which is passed to the inner loop.
171
 * buffersize  : Buffer size for the iterator. For the default, pass in 0.
172
 * funcname    : The name of the reduction function, for error messages.
173
 * errormask   : forwarded from _get_bufsize_errmask
174
 *
175
 * TODO FIXME: if you squint, this is essentially an second independent
176
 * implementation of generalized ufuncs with signature (i)->(), plus a few
177
 * extra bells and whistles. (Indeed, as far as I can tell, it was originally
178
 * split out to support a fancy version of count_nonzero... which is not
179
 * actually a reduction function at all, it's just a (i)->() function!) So
180
 * probably these two implementation should be merged into one. (In fact it
181
 * would be quite nice to support axis= and keepdims etc. for arbitrary
182
 * generalized ufuncs!)
183
 */
184
NPY_NO_EXPORT PyArrayObject *
185 1
PyUFunc_ReduceWrapper(PyArrayObject *operand, PyArrayObject *out,
186
                      PyArrayObject *wheremask,
187
                      PyArray_Descr *operand_dtype,
188
                      PyArray_Descr *result_dtype,
189
                      NPY_CASTING casting,
190
                      npy_bool *axis_flags, int reorderable,
191
                      int keepdims,
192
                      PyObject *identity,
193
                      PyArray_ReduceLoopFunc *loop,
194
                      void *data, npy_intp buffersize, const char *funcname,
195
                      int errormask)
196
{
197 1
    PyArrayObject *result = NULL;
198 1
    npy_intp skip_first_count = 0;
199

200
    /* Iterator parameters */
201 1
    NpyIter *iter = NULL;
202
    PyArrayObject *op[3];
203
    PyArray_Descr *op_dtypes[3];
204
    npy_uint32 flags, op_flags[3];
205

206
    /* More than one axis means multiple orders are possible */
207 1
    if (!reorderable && count_axes(PyArray_NDIM(operand), axis_flags) > 1) {
208 1
        PyErr_Format(PyExc_ValueError,
209
                     "reduction operation '%s' is not reorderable, "
210
                     "so at most one axis may be specified",
211
                     funcname);
212 1
        return NULL;
213
    }
214
    /* Can only use where with an initial ( from identity or argument) */
215 1
    if (wheremask != NULL && identity == Py_None) {
216 1
        PyErr_Format(PyExc_ValueError,
217
                     "reduction operation '%s' does not have an identity, "
218
                     "so to use a where mask one has to specify 'initial'",
219
                     funcname);
220 1
        return NULL;
221
    }
222

223

224
    /* Set up the iterator */
225 1
    op[0] = out;
226 1
    op[1] = operand;
227 1
    op_dtypes[0] = result_dtype;
228 1
    op_dtypes[1] = operand_dtype;
229

230 1
    flags = NPY_ITER_BUFFERED |
231
            NPY_ITER_EXTERNAL_LOOP |
232
            NPY_ITER_GROWINNER |
233
            NPY_ITER_DONT_NEGATE_STRIDES |
234
            NPY_ITER_ZEROSIZE_OK |
235
            NPY_ITER_REFS_OK |
236
            NPY_ITER_DELAY_BUFALLOC |
237
            NPY_ITER_COPY_IF_OVERLAP;
238 1
    op_flags[0] = NPY_ITER_READWRITE |
239
                  NPY_ITER_ALIGNED |
240
                  NPY_ITER_ALLOCATE |
241
                  NPY_ITER_NO_SUBTYPE;
242 1
    op_flags[1] = NPY_ITER_READONLY |
243
                  NPY_ITER_ALIGNED |
244
                  NPY_ITER_NO_BROADCAST;
245

246 1
    if (wheremask != NULL) {
247 1
        op[2] = wheremask;
248
        /* wheremask is guaranteed to be NPY_BOOL, so borrow its reference */
249 1
        op_dtypes[2] = PyArray_DESCR(wheremask);
250
        assert(op_dtypes[2]->type_num == NPY_BOOL);
251 1
        if (op_dtypes[2] == NULL) {
252
            goto fail;
253
        }
254 1
        op_flags[2] = NPY_ITER_READONLY;
255
    }
256
    /* Set up result array axes mapping, operand and wheremask use default */
257
    int result_axes[NPY_MAXDIMS];
258 1
    int *op_axes[3] = {result_axes, NULL, NULL};
259

260 1
    int curr_axis = 0;
261 1
    for (int i = 0; i < PyArray_NDIM(operand); i++) {
262 1
        if (axis_flags[i]) {
263 1
            if (keepdims) {
264 1
                result_axes[i] = NPY_ITER_REDUCTION_AXIS(curr_axis);
265 1
                curr_axis++;
266
            }
267
            else {
268 1
                result_axes[i] = NPY_ITER_REDUCTION_AXIS(-1);
269
            }
270
        }
271
        else {
272 1
            result_axes[i] = curr_axis;
273 1
            curr_axis++;
274
        }
275
    }
276 1
    if (out != NULL) {
277
        /* NpyIter does not raise a good error message in this common case. */
278 1
        if (NPY_UNLIKELY(curr_axis != PyArray_NDIM(out))) {
279 1
            if (keepdims) {
280 1
                PyErr_Format(PyExc_ValueError,
281
                        "output parameter for reduction operation %s has the "
282
                        "wrong number of dimensions: Found %d but expected %d "
283
                        "(must match the operand's when keepdims=True)",
284
                        funcname, PyArray_NDIM(out), curr_axis);
285
            }
286
            else {
287 1
                PyErr_Format(PyExc_ValueError,
288
                        "output parameter for reduction operation %s has the "
289
                        "wrong number of dimensions: Found %d but expected %d",
290
                        funcname, PyArray_NDIM(out), curr_axis);
291
            }
292
            goto fail;
293
        }
294
    }
295

296 1
    iter = NpyIter_AdvancedNew(wheremask == NULL ? 2 : 3, op, flags,
297
                               NPY_KEEPORDER, casting,
298
                               op_flags,
299
                               op_dtypes,
300
                               PyArray_NDIM(operand), op_axes, NULL, buffersize);
301 1
    if (iter == NULL) {
302
        goto fail;
303
    }
304

305 1
    result = NpyIter_GetOperandArray(iter)[0];
306

307
    /*
308
     * Initialize the result to the reduction unit if possible,
309
     * otherwise copy the initial values and get a view to the rest.
310
     */
311

312 1
    if (identity != Py_None) {
313 1
        if (PyArray_FillWithScalar(result, identity) < 0) {
314
            goto fail;
315
        }
316
    }
317
    else {
318
        /*
319
         * For 1-D skip_first_count could be optimized to 0, but no-identity
320
         * reductions are not super common.
321
         * (see also comment in CopyInitialReduceValues)
322
         */
323 1
        skip_first_count = PyArray_CopyInitialReduceValues(
324
                result, operand, axis_flags, funcname, keepdims);
325 1
        if (skip_first_count < 0) {
326
            goto fail;
327
        }
328
    }
329

330 1
    if (!NpyIter_Reset(iter, NULL)) {
331
        goto fail;
332
    }
333

334
    /* Start with the floating-point exception flags cleared */
335 1
    npy_clear_floatstatus_barrier((char*)&iter);
336

337 1
    if (NpyIter_GetIterSize(iter) != 0) {
338
        NpyIter_IterNextFunc *iternext;
339
        char **dataptr;
340
        npy_intp *strideptr;
341
        npy_intp *countptr;
342
        int needs_api;
343

344 1
        iternext = NpyIter_GetIterNext(iter, NULL);
345 1
        if (iternext == NULL) {
346
            goto fail;
347
        }
348 1
        dataptr = NpyIter_GetDataPtrArray(iter);
349 1
        strideptr = NpyIter_GetInnerStrideArray(iter);
350 1
        countptr = NpyIter_GetInnerLoopSizePtr(iter);
351

352 1
        needs_api = NpyIter_IterationNeedsAPI(iter);
353

354
        /* Straightforward reduction */
355 1
        if (loop == NULL) {
356 0
            PyErr_Format(PyExc_RuntimeError,
357
                    "reduction operation %s did not supply an "
358
                    "inner loop function", funcname);
359 0
            goto fail;
360
        }
361

362 1
        if (loop(iter, dataptr, strideptr, countptr,
363
                        iternext, needs_api, skip_first_count, data) < 0) {
364
            goto fail;
365
        }
366
    }
367

368
    /* Check whether any errors occurred during the loop */
369 1
    if (PyErr_Occurred() ||
370 1
            _check_ufunc_fperr(errormask, NULL, "reduce") < 0) {
371
        goto fail;
372
    }
373

374 1
    if (out != NULL) {
375 1
        result = out;
376
    }
377 1
    Py_INCREF(result);
378

379 1
    if (!NpyIter_Deallocate(iter)) {
380 0
        Py_DECREF(result);
381
        return NULL;
382
    }
383
    return result;
384

385 1
fail:
386 1
    if (iter != NULL) {
387 1
        NpyIter_Deallocate(iter);
388
    }
389

390
    return NULL;
391
}

Read our documentation on viewing source code .

Loading