Repository: scipy/scipy
Pull request: BUG: sparse: fix argmin/argmax when all entries are nonzero
#### Reference issue
Fixes gh-16929.

#### What does this implement/fix?
The special case for "all nonzeros" was incorrect, returning the flat index of the `data` array instead of the sparse matrix itself. See the linked issue for examples.


Problem statement:
BUG: `scipy.sparse.csc_matrix.argmin` returns wrong values
### Describe your issue.

`scipy.sparse.csc_matrix.argmin` seems to return wrong value when all elements are positive.

### Reproducing Code Example

```python
>>> import scipy
>>> x = scipy.sparse.coo_matrix(([1, 0.1, 1, 1], ([0, 0, 1, 1], [0, 1, 0, 1])), shape=(2, 2)).tocsc()
>>> x.todense()
matrix([[1. , 0.1],
        [1. , 1. ]])
>>> x.argmin()
2
>>> x.todense().argmin()
1
```


### Error message

```shell
No error messages.
```


### SciPy/NumPy/Python version information

1.9.1 1.23.0 sys.version_info(major=3, minor=9, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
Thanks for the report, @asi1024.  Here are a few more examples:
```
In [75]: a = csc_matrix([[3, 4, 1], [10, 2, 9]])

In [76]: a.A
Out[76]: 
array([[ 3,  4,  1],
       [10,  2,  9]])

In [77]: a.argmin()  # Should be 2
Out[77]: 4

In [78]: a.A.argmin()
Out[78]: 2
```

```
In [79]: b = csr_matrix([[3, 4, 1], [10, 2, 9]])

In [80]: b.A
Out[80]: 
array([[ 3,  4,  1],
       [10,  2,  9]])

In [81]: b.argmin()  # Should be 2
Out[81]: 4

In [82]: b.A.argmin()
Out[82]: 2
```

```
In [83]: c = coo_matrix(([10, 9, 8, 0, 6, 5], ([0, 0, 1, 1, 2, 2], [0, 1, 0, 1, 0, 1])))

In [84]: c.A
Out[84]: 
array([[10,  9],
       [ 8,  0],
       [ 6,  5]])

In [85]: c.argmin()  # Should be 3
Out[85]: 4

In [86]: c.A.argmin()
Out[86]: 3
```
The above was run using the main development branch:
```
In [87]: import scipy

In [88]: scipy.__version__
Out[88]: '1.10.0.dev0+1757.cc63003'
```
The bug is here: https://github.com/scipy/scipy/blob/cc63003326b7ae2e9dd41cd6af0e6e3079c32b95/scipy/sparse/_data.py#L282

If the matrix is entirely nonzero (i.e., `nnz == prod(shape)`) then we exit early with the untransformed argmin result. From a first look, it appears that the correct return value for that case should be `mat.row[am] * mat.shape[1] + mat.col[am]`.
I'll put a fix together soon.

Commit message: Change COO.sum_duplicates to produce a "C-ordered" result

Previously, we were sorting the indices in column-major order, which meant that the linearized form would _not_ be sorted.

Patch to review:
diff --git a/scipy/sparse/_coo.py b/scipy/sparse/_coo.py
index 371329e069c3..652584d2e3f8 100644
--- a/scipy/sparse/_coo.py
+++ b/scipy/sparse/_coo.py
@@ -533,7 +533,7 @@ def _sum_duplicates(self, row, col, data):
         # Assumes (data, row, col) not in canonical format.
         if len(data) == 0:
             return row, col, data
-        order = np.lexsort((row, col))
+        order = np.lexsort((col, row))
         row = row[order]
         col = col[order]
         data = data[order]
diff --git a/scipy/sparse/_data.py b/scipy/sparse/_data.py
index 5f4fedfed49c..6f3b9e0f06cd 100644
--- a/scipy/sparse/_data.py
+++ b/scipy/sparse/_data.py
@@ -261,37 +261,37 @@ def _arg_min_or_max(self, axis, out, op, compare):
 
         validateaxis(axis)
 
-        if axis is None:
-            if 0 in self.shape:
-                raise ValueError("Can't apply the operation to "
-                                 "an empty matrix.")
+        if axis is not None:
+            return self._arg_min_or_max_axis(axis, op, compare)
 
-            if self.nnz == 0:
-                return 0
-            else:
-                zero = self.dtype.type(0)
-                mat = self.tocoo()
-                mat.sum_duplicates()
-                am = op(mat.data)
-                m = mat.data[am]
-
-                if compare(m, zero):
-                    # cast to Python int to avoid overflow
-                    # and RuntimeError
-                    return int(mat.row[am])*mat.shape[1] + int(mat.col[am])
-                else:
-                    size = np.prod(mat.shape)
-                    if size == mat.nnz:
-                        return am
-                    else:
-                        ind = mat.row * mat.shape[1] + mat.col
-                        zero_ind = _find_missing_index(ind, size)
-                        if m == zero:
-                            return min(zero_ind, am)
-                        else:
-                            return zero_ind
-
-        return self._arg_min_or_max_axis(axis, op, compare)
+        if 0 in self.shape:
+            raise ValueError("Can't apply the operation to "
+                             "an empty matrix.")
+
+        if self.nnz == 0:
+            return 0
+
+        zero = self.dtype.type(0)
+        mat = self.tocoo()
+        mat.sum_duplicates()
+        am = op(mat.data)
+        m = mat.data[am]
+        num_row, num_col = mat.shape
+
+        if compare(m, zero):
+            # cast to Python int to avoid overflow and RuntimeError
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        size = num_row * num_col
+        if size == mat.nnz:
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        ind = mat.row * num_col + mat.col
+        ind.sort()
+        zero_ind = _find_missing_index(ind, size)
+        if m == zero:
+            return min(zero_ind, am)
+        return zero_ind
 
     def max(self, axis=None, out=None):
         """
diff --git a/scipy/sparse/tests/test_base.py b/scipy/sparse/tests/test_base.py
index 39014a9e18dc..86f68e2b1ec7 100644
--- a/scipy/sparse/tests/test_base.py
+++ b/scipy/sparse/tests/test_base.py
@@ -3464,6 +3464,7 @@ def test_numpy_minmax(self):
             assert_array_equal(np.max(datsp), np.max(dat))
 
     def test_argmax(self):
+        from scipy.sparse import _data
         D1 = np.array([
             [-1, 5, 2, 3],
             [0, 0, -1, -2],
@@ -3472,9 +3473,14 @@ def test_argmax(self):
             [1, 2, 0, 0],
         ])
         D2 = D1.transpose()
+        D3 = np.array([[4, 3], [7, 5]])
+        D4 = np.array([[4, 3], [7, 0]])
+        D5 = np.array([[5, 5, 3], [4, 9, 10], [3, 4, 9]])
 
-        for D in [D1, D2]:
-            mat = csr_matrix(D)
+        for D in [D1, D2, D3, D4, D5]:
+            mat = self.spmatrix(D)
+            if not isinstance(mat, _data._minmax_mixin):
+                continue
 
             assert_equal(mat.argmax(), np.argmax(D))
             assert_equal(mat.argmin(), np.argmin(D))


Code review comments:
diff_hunk: @@ -261,37 +261,37 @@ def _arg_min_or_max(self, axis, out, op, compare):
 
         validateaxis(axis)
 
-        if axis is None:
-            if 0 in self.shape:
-                raise ValueError("Can't apply the operation to "
-                                 "an empty matrix.")
+        if axis is not None:
+            return self._arg_min_or_max_axis(axis, op, compare)
 
-            if self.nnz == 0:
-                return 0
-            else:
-                zero = self.dtype.type(0)
-                mat = self.tocoo()
-                mat.sum_duplicates()
-                am = op(mat.data)
-                m = mat.data[am]
-
-                if compare(m, zero):
-                    # cast to Python int to avoid overflow
-                    # and RuntimeError
-                    return int(mat.row[am])*mat.shape[1] + int(mat.col[am])
-                else:
-                    size = np.prod(mat.shape)
-                    if size == mat.nnz:
-                        return am
-                    else:
-                        ind = mat.row * mat.shape[1] + mat.col
-                        zero_ind = _find_missing_index(ind, size)
-                        if m == zero:
-                            return min(zero_ind, am)
-                        else:
-                            return zero_ind
-
-        return self._arg_min_or_max_axis(axis, op, compare)
+        if 0 in self.shape:
+            raise ValueError("Can't apply the operation to "
+                             "an empty matrix.")
+
+        if self.nnz == 0:
+            return 0
+
+        zero = self.dtype.type(0)
+        mat = self.tocoo()
+        mat.sum_duplicates()
+        am = op(mat.data)
+        m = mat.data[am]
+        num_row, num_col = mat.shape
+
+        if compare(m, zero):
+            # cast to Python int to avoid overflow and RuntimeError
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        size = num_row * num_col
+        if size == mat.nnz:
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        ind = mat.row * num_col + mat.col
path: scipy/sparse/_data.py
text: @user1:
If this suggestion is accepted, the logic bellow needs to be adapted accordingly.
```suggestion
        # Linearised indices for the C-order.
        indices = mat.row * num_col + mat.col
```

diff_hunk: @@ -261,37 +261,37 @@ def _arg_min_or_max(self, axis, out, op, compare):
 
         validateaxis(axis)
 
-        if axis is None:
-            if 0 in self.shape:
-                raise ValueError("Can't apply the operation to "
-                                 "an empty matrix.")
+        if axis is not None:
+            return self._arg_min_or_max_axis(axis, op, compare)
 
-            if self.nnz == 0:
-                return 0
-            else:
-                zero = self.dtype.type(0)
-                mat = self.tocoo()
-                mat.sum_duplicates()
-                am = op(mat.data)
-                m = mat.data[am]
-
-                if compare(m, zero):
-                    # cast to Python int to avoid overflow
-                    # and RuntimeError
-                    return int(mat.row[am])*mat.shape[1] + int(mat.col[am])
-                else:
-                    size = np.prod(mat.shape)
-                    if size == mat.nnz:
-                        return am
-                    else:
-                        ind = mat.row * mat.shape[1] + mat.col
-                        zero_ind = _find_missing_index(ind, size)
-                        if m == zero:
-                            return min(zero_ind, am)
-                        else:
-                            return zero_ind
-
-        return self._arg_min_or_max_axis(axis, op, compare)
+        if 0 in self.shape:
+            raise ValueError("Can't apply the operation to "
+                             "an empty matrix.")
path: scipy/sparse/_data.py
text: @user1:
```suggestion
            raise ValueError("Can't apply the operation to an empty matrix.")
```

diff_hunk: @@ -261,37 +261,37 @@ def _arg_min_or_max(self, axis, out, op, compare):
 
         validateaxis(axis)
 
-        if axis is None:
-            if 0 in self.shape:
-                raise ValueError("Can't apply the operation to "
-                                 "an empty matrix.")
+        if axis is not None:
+            return self._arg_min_or_max_axis(axis, op, compare)
 
-            if self.nnz == 0:
-                return 0
-            else:
-                zero = self.dtype.type(0)
-                mat = self.tocoo()
-                mat.sum_duplicates()
-                am = op(mat.data)
-                m = mat.data[am]
-
-                if compare(m, zero):
-                    # cast to Python int to avoid overflow
-                    # and RuntimeError
-                    return int(mat.row[am])*mat.shape[1] + int(mat.col[am])
-                else:
-                    size = np.prod(mat.shape)
-                    if size == mat.nnz:
-                        return am
-                    else:
-                        ind = mat.row * mat.shape[1] + mat.col
-                        zero_ind = _find_missing_index(ind, size)
-                        if m == zero:
-                            return min(zero_ind, am)
-                        else:
-                            return zero_ind
-
-        return self._arg_min_or_max_axis(axis, op, compare)
+        if 0 in self.shape:
+            raise ValueError("Can't apply the operation to "
+                             "an empty matrix.")
+
+        if self.nnz == 0:
+            return 0
+
+        zero = self.dtype.type(0)
+        mat = self.tocoo()
+        mat.sum_duplicates()
+        am = op(mat.data)
+        m = mat.data[am]
+        num_row, num_col = mat.shape
+
+        if compare(m, zero):
+            # cast to Python int to avoid overflow and RuntimeError
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        size = num_row * num_col
+        if size == mat.nnz:
+            return int(mat.row[am]) * num_col + int(mat.col[am])
+
+        ind = mat.row * num_col + mat.col
+        ind.sort()
+        zero_ind = _find_missing_index(ind, size)
+        if m == zero:
+            return min(zero_ind, am)
path: scipy/sparse/_data.py
text: @user1:
```suggestion
        # At this stage, the max or min is >= 0.
        first_implicit_zero_index = _find_missing_index(ind, size)
        if m == zero:
            return min(first_implicit_zero_index, am)
```

diff_hunk: @@ -533,7 +533,7 @@ def _sum_duplicates(self, row, col, data):
         # Assumes (data, row, col) not in canonical format.
         if len(data) == 0:
             return row, col, data
-        order = np.lexsort((row, col))
+        order = np.lexsort((col, row))
path: scipy/sparse/_coo.py
text: @user1:
```suggestion
        # First sorting wr.t. to rows, and then on columns.
        # This allows being consistent with NumPy which identifies
        # the argmin or argmax as the first element in the C-order
        # if there are several argmin or argmax candidates.
        order = np.lexsort((col, row))
```