Repository: numpy/numpy
Pull request: MAINT: Reduce allocation size of empty (0 size) arrays to 1 byte
Restart #15788 from @eric-wieser.
- Allocate a minimum 1 byte when the data chunk should be 0 bytes, to avoid undefined behaviour
- Make allocation match the free calculation of how much memory to release. This closes #16410 and #20289

Problem statement:
empty, zeros, ones, full allocate product of non-zero size elements, even if final array is zero-sized
<!-- Please describe the issue in detail here, and fill in the field
s below -->

Functions like `np.zeros`, `np.empty`, `np.ones` and `np.full` can be used to allocate large arrays easily. However, they always allocate memory, even if the array contains no data because one of the shape elements is 0. They, in particular, seem to allocate memory ignoring the 0s in the shape (`np.product([x for x in shape if x != 0])`).

If any of the shape elements in zero, there's no possibility of any memory accesses ever, so the allocation seems unnecessary?

If the shape has a zero, one can get an equivalent array using `broadcast_to`, which won't allocate.

### Reproducing code example:

<!-- A short code example that reproduces the problem/missing feature. It should be
self-contained, i.e., possible to run as-is via 'python myproblem.py' -->

```python
import numpy as np
import tracemalloc as tm

def test(f):
    tm.start()
    arr = f((0, 100, 0, 200, 0, 300, 0), dtype=np.uint8)
    mem, _ = tm.get_traced_memory()
    tm.stop()
    print(f"{f.__name__}: allocated memory: {mem} bytes, array size: {arr.size}")

def full(shape, dtype):
    return np.full(shape, 123, dtype=dtype)

test(np.zeros)
test(np.empty)
test(np.ones)
test(full)

# workaround, use broadcasting
def broadcast(shape, dtype):
    return np.broadcast_to(dtype(), shape)
test(broadcast)
```

```
zeros: allocated memory: 6000248 bytes, array size: 0
empty: allocated memory: 6000080 bytes, array size: 0
ones: allocated memory: 6000080 bytes, array size: 0
full: allocated memory: 6000080 bytes, array size: 0
broadcast: allocated memory: 3680 bytes, array size: 0
```

### Numpy/Python version information:

- `numpy.__version__`: 1.18.4 
- `sys.version`: 3.6.9 (default, Jul 10 2019, 12:25:55) [GCC 4.2.1 Compatible Apple LLVM 10.0.1 (clang-1001.0.46.4)]

Hints from issue discussion:
I think @mattip and I both have (had?) open PRs about this.
Right, gh-15780 was rejected since it hit undefined behaviour in the C spec, and gh-15788 is still marked "Draft"
Nice!

I noticed that this applies to (fancy) indexing along a non-zero sized axis too, e.g.:

```python
import numpy as np
import tracemalloc as tm

base = np.broadcast_to(np.uint8(0), (0, 100, 0, 200, 0, 300, 0))
idx = np.arange(100)

tm.start()
arr = base[:, idx, ...]
mem, _ = tm.get_traced_memory()
tm.stop()

print(f"allocated memory: {mem} bytes, array size: {arr.size}")
```
```
allocated memory: 6000480 bytes, array size: 0
```

Can this be considered as part of this issue/those PRs, or should I open a separate one?
Yes, there is no difference with fancy indexing, since in the end it also allocates a normal array. 
So it is a duplicate of gh-15753 and those PRs address it.
I fine with either proposal (just because at some point this was likely intentionally, does not mean it was good), but to be honest I (EDIT *+not*) very enthusiastic about trying to figure out whether or not its a good idea. (I suppose right now I fail to see a bit that this memory loss is actually a big deal. The memory is never touched and in most use cases I would expect there to be much larger arrays flying around somewhere, it seems more like something you run into when running crazy test suits)

Commit message: BUG: Fix FromIter stride resetting logic

We have to reset all strides in the case of a subarray dtype which
may have caused the result array to not be 1-D.  Since the logic
is a bit complicated now, restore the logic to only do this when it
really is necessary.

Patch to review:
diff --git a/numpy/core/src/multiarray/ctors.c b/numpy/core/src/multiarray/ctors.c
index c0e80d1eeee0..1e6f8ee249cb 100644
--- a/numpy/core/src/multiarray/ctors.c
+++ b/numpy/core/src/multiarray/ctors.c
@@ -763,14 +763,16 @@ PyArray_NewFromDescr_int(
          * calls, `PyArray_NBYTES_ALLOCATED` is a private helper matching this
          * behaviour, but without overflow checking.
          */
+        int is_zero = 0;
         for (int i = 0; i < nd; i++) {
             fa->dimensions[i] = dims[i];
 
             if (fa->dimensions[i] == 0) {
                 /*
-                 * Compare to PyArray_OverflowMultiplyList that
-                 * returns 0 in this case. See also `PyArray_NBYTES_ALLOCATED`.
+                 * Continue calculating the max size "as if" this were 1
+                 * to get the proper overflow error
                  */
+                is_zero = 1;
                 continue;
             }
 
@@ -791,6 +793,9 @@ PyArray_NewFromDescr_int(
                 goto fail;
             }
         }
+        if (is_zero) {
+            nbytes = 0;
+        }
 
         /* Fill the strides (or copy them if they were passed in) */
         if (strides == NULL) {
@@ -825,11 +830,13 @@ PyArray_NewFromDescr_int(
          * Allocate something even for zero-space arrays
          * e.g. shape=(0,) -- otherwise buffer exposure
          * (a.data) doesn't work as it should.
-         * Could probably just allocate a few bytes here. -- Chuck
-         * Note: always sync this with calls to PyDataMem_UserFREE
          */
         if (nbytes == 0) {
-            nbytes = descr->elsize ? descr->elsize : 1;
+            nbytes = 1;
+            /* Make sure all the strides are 0 */
+            for (int i = 0; i < nd; i++) {
+                fa->strides[i] = 0;
+            }
         }
         /*
          * It is bad to have uninitialized OBJECT pointers
@@ -3936,20 +3943,13 @@ PyArray_FromIter(PyObject *obj, PyArray_Descr *dtype, npy_intp count)
     if (ret == NULL) {
         goto done;
     }
-#ifdef NPY_RELAXED_STRIDES_DEBUG
-    /* Incompatible with NPY_RELAXED_STRIDES_DEBUG due to growing */
-    if (elcount == 1) {
-        PyArray_STRIDES(ret)[0] = elsize;
-    }
-#endif /* NPY_RELAXED_STRIDES_DEBUG */
-
 
     char *item = PyArray_BYTES(ret);
     for (i = 0; i < count || count == -1; i++, item += elsize) {
         PyObject *value = PyIter_Next(iter);
         if (value == NULL) {
             if (PyErr_Occurred()) {
-                /* Fetching next item failed rather than exhausting iterator */
+                /* Fetching next item failed perhaps due to exhausting iterator */
                 goto done;
             }
             break;
@@ -3995,7 +3995,6 @@ PyArray_FromIter(PyObject *obj, PyArray_Descr *dtype, npy_intp count)
         Py_DECREF(value);
     }
 
-
     if (i < count) {
         PyErr_Format(PyExc_ValueError,
                 "iterator too short: Expected %zd but iterator had only %zd "
@@ -4022,6 +4021,21 @@ PyArray_FromIter(PyObject *obj, PyArray_Descr *dtype, npy_intp count)
             goto done;
         }
         ((PyArrayObject_fields *)ret)->data = new_data;
+
+        if (count < 0 || NPY_RELAXED_STRIDES_DEBUG) {
+            /*
+             * If the count was smaller than zero or NPY_RELAXED_STRIDES_DEBUG
+             * was active, the strides may be all 0 or intentionally mangled
+             * (even in the later dimensions for `count < 0`!
+             * Thus, fix all strides here again for C-contiguity.
+             */
+            int oflags;
+            _array_fill_strides(
+                    PyArray_STRIDES(ret), PyArray_DIMS(ret), PyArray_NDIM(ret),
+                    PyArray_ITEMSIZE(ret), NPY_ARRAY_C_CONTIGUOUS, &oflags);
+            PyArray_STRIDES(ret)[0] = elsize;
+            assert(oflags & NPY_ARRAY_C_CONTIGUOUS);
+        }
     }
     PyArray_DIMS(ret)[0] = i;
 
diff --git a/numpy/core/src/multiarray/getset.c b/numpy/core/src/multiarray/getset.c
index d6406845a186..dded9a031c40 100644
--- a/numpy/core/src/multiarray/getset.c
+++ b/numpy/core/src/multiarray/getset.c
@@ -513,6 +513,7 @@ array_descr_set(PyArrayObject *self, PyObject *arg, void *NPY_UNUSED(ignored))
         /* resize on last axis only */
         int axis = PyArray_NDIM(self) - 1;
         if (PyArray_DIMS(self)[axis] != 1 &&
+                PyArray_STRIDES(self)[axis] != 0 &&
                 PyArray_STRIDES(self)[axis] != PyArray_DESCR(self)->elsize) {
             PyErr_SetString(PyExc_ValueError,
                     "To change to a dtype of a different size, the last axis "
diff --git a/numpy/core/tests/test_numeric.py b/numpy/core/tests/test_numeric.py
index 165fcbce695d..a7f5df472ffc 100644
--- a/numpy/core/tests/test_numeric.py
+++ b/numpy/core/tests/test_numeric.py
@@ -1218,18 +1218,26 @@ def test_empty_not_structured(self, dtype):
         with pytest.raises(ValueError, match="Must specify length"):
             np.fromiter([], dtype=dtype)
 
-    @pytest.mark.parametrize("dtype",
-            # Note that `np.dtype(("O", (10, 5)))` is a subarray dtype
-            ["d", "i,O", np.dtype(("O", (10, 5))), "O"])
-    def test_growth_and_complicated_dtypes(self, dtype):
+    @pytest.mark.parametrize(["dtype", "data"],
+            [("d", [1, 2, 3, 4, 5, 6, 7, 8, 9]),
+             ("O", [1, 2, 3, 4, 5, 6, 7, 8, 9]),
+             ("i,O", [(1, 2), (5, 4), (2, 3), (9, 8), (6, 7)]),
+             # subarray dtypes (important because their dimensions end up
+             # in the result arrays dimension:
+             ("2i", [(1, 2), (5, 4), (2, 3), (9, 8), (6, 7)]),
+             (np.dtype(("O", (2, 3))),
+              [((1, 2, 3), (3, 4, 5)), ((3, 2, 1), (5, 4, 3))])])
+    @pytest.mark.parametrize("length_hint", [0, 1])
+    def test_growth_and_complicated_dtypes(self, dtype, data, length_hint):
         dtype = np.dtype(dtype)
-        data = [1, 2, 3, 4, 5, 6, 7, 8, 9] * 100  # make sure we realloc a bit
+
+        data = data * 100  # make sure we realloc a bit
 
         class MyIter:
             # Class/example from gh-15789
             def __length_hint__(self):
                 # only required to be an estimate, this is legal
-                return 1
+                return length_hint  # 0 or 1
 
             def __iter__(self):
                 return iter(data)
diff --git a/numpy/core/tests/test_regression.py b/numpy/core/tests/test_regression.py
index 0ebdfec612c0..e578491b31a7 100644
--- a/numpy/core/tests/test_regression.py
+++ b/numpy/core/tests/test_regression.py
@@ -430,7 +430,6 @@ def __getitem__(self, key):
     def test_lexsort_zerolen_custom_strides(self):
         # Ticket #14228
         xs = np.array([], dtype='i8')
-        assert xs.strides == (8,)
         assert np.lexsort((xs,)).shape[0] == 0 # Works
 
         xs.strides = (16,)


Code review comments:
diff_hunk: @@ -513,6 +513,7 @@ array_descr_set(PyArrayObject *self, PyObject *arg, void *NPY_UNUSED(ignored))
         /* resize on last axis only */
         int axis = PyArray_NDIM(self) - 1;
         if (PyArray_DIMS(self)[axis] != 1 &&
+                PyArray_STRIDES(self)[axis] != 0 &&
path: numpy/core/src/multiarray/getset.c
text: @user1:
I think it would be worth understanding the failure. The `PyArray_DIMS(self)[axis] != 1 && PyArray_STRIDES(self)[axis] != 0` logic feels wrong to me, as it means zero-size arrays may or may not reach this error depending on their stride (which can be anything as long as their allocation is ok).