Repository: pydata/sparse (python)

Issue:
Constructing GCXS from non-canonical `scipy.sparse.csr_matrix` results in wrong results
**Describe the bug**

GCXS seems to require a canonical CSR data structure, i.e. column indices sorted within a row and without duplicates.  Slicing a GCXS array with non-canonical internal data structure gives wrong results. The requirement for canonical data structures is not documented: https://sparse.pydata.org/en/stable/generated/sparse.GCXS.html

`scipy.sparse.csr_matrix` doesn't require its data structure to be canonical. The GCXS constructor doesn't seem to check if a passed `scipy.sparse.csr_matrix` is canonical. In effect, constructing a GCXS array from a perfectly valid `scipy.sparse.csr_matrix` may create a broken GCXS array which gives wrong numerical results. Using the explicit `GCXS.from_scipy_sparse()` gives the same behavior.

**To Reproduce**

This example uses `cupyx.scipy.sparse.csr_matrix` since that produces non-canonical CSR data structures when constructed from dense arrays. This highlights the severity of this bug where normal use of the APIs of these packages gives numerically wrong results without causing any warning or error. In my understanding the bug is mostly in GCXS since `scipy.sparse.csr_matrix` and  `cupyx.scipy.sparse.csr_matrix` work perfectly fine with non-canonical data structures.

```python
import numpy as np
import cupy
import cupyx.scipy.sparse as csp
import sparse

# Just some data
source = np.random.random((23, 42))
# Reference slice
ref = source[2:4]

# We make a cupyx.scipy.sparse.csr_matrix from the source data
cupyx_sparse = csp.csr_matrix(cupy.array(source))
# It is not canonical, i.e. has duplicates or unsorted indices
print("Canonical?",  cupyx_sparse.has_canonical_format)

# We slice it
r1 = cupy.asnumpy(cupyx_sparse[2:4].todense())
# Seems to work correctly
print("max difference cupyx.scipy.sparse", np.max(np.abs(ref - r1)))

# We make a GCXS matrix by first transferring from the GPU and then
# constructing GCXS
gcxs = sparse.GCXS(cupyx_sparse.get())
# We slice it
r2 = gcxs[2:4].todense()
# It is wrong
print("max difference GCXS", np.max(np.abs(ref - r2)))
```

Edit: Include output

```
Canonical? False
max difference cupyx.scipy.sparse 0.0
max difference GCXS 0.9852408955033491
```

**Expected behavior**

The expected behavior can be split into two aspects, of which one is preventing a bug with serious impact on downstream code, and the other a feature request.

First, constructing a GCXS array from a valid `scipy.sparse.csr_matrix` that GCXS can't handle correctly should throw an error instead of producing a GCXS array that gives wrong results. 

Second, it would be good if GCXS could be constructed successfully from a non-canonical `scipy.sparse.csr_matrix`. 

**System**
 - OS and version: Win 11
 - `sparse` version (`sparse.__version__` '0.14.0')
 - NumPy version (`np.__version__` '1.24.4')
 - Numba version (`numba.__version__` '0.57.1')
 - CuPy version (`cupy.__version__` '12.1.0')

**Additional context**

At https://github.com/LiberTEM/sparseconverter/pull/14 I am currently working on tests and workarounds. Perhaps the test matrix in https://github.com/LiberTEM/sparseconverter/blob/main/tests/test_sparseconverters.py could be useful here?


Pull request:
fix: Fix conversions from non-canonical `scipy.sparse` arrays
<!--
# Thanks for contributing a pull request!
## Please make sure you see our contribution guidelines: https://github.com/pydata/sparse/blob/main/docs/contributing.md

Your PR title should start with any of these abbreviatons: `build`, `chore`, `ci`, `depr`, `docs`, `feat`, `fix`, `perf`, `refactor`, `release`, `test`. Add a `!`at the end, if it is a breaking change.
-->

## What type of PR is this? (check all applicable)

- [ ] 💾 Refactor
- [ ] 🪄 Feature
- [x] 🐞 Bug Fix
- [ ] 🔧 Optimization
- [ ] 📚 Documentation
- [ ] 🧪 Test
- [ ] 🛠️ Other

## Related issues

- Closes #602

## Checklist

- [x] Code follows style guide
- [x] Tests added
- [ ] Documented the changes

***


Patch:
diff --git a/sparse/numba_backend/_compressed/compressed.py b/sparse/numba_backend/_compressed/compressed.py
index 6309df0..0b87b9a 100644
--- a/sparse/numba_backend/_compressed/compressed.py
+++ b/sparse/numba_backend/_compressed/compressed.py
@@ -209,11 +209,14 @@ class GCXS(SparseArray, NDArrayOperatorsMixin):
 
     @classmethod
     def from_scipy_sparse(cls, x, /, *, fill_value=None):
-        if x.format == "csc":
-            return cls((x.data, x.indices, x.indptr), shape=x.shape, compressed_axes=(1,), fill_value=fill_value)
-
-        x = x.asformat("csr")
-        return cls((x.data, x.indices, x.indptr), shape=x.shape, compressed_axes=(0,), fill_value=fill_value)
+        is_csc = x.format == "csc"
+        ca = (1,) if is_csc else (0,)
+        if not is_csc:
+            x = x.asformat("csr")
+        if not x.has_canonical_format:
+            x.eliminate_zeros()
+            x.sum_duplicates()
+        return cls((x.data, x.indices, x.indptr), shape=x.shape, compressed_axes=ca, fill_value=fill_value)
 
     @classmethod
     def from_iter(cls, x, shape=None, compressed_axes=None, fill_value=None, idx_dtype=None):
@@ -903,6 +906,10 @@ class CSR(_Compressed2d):
     @classmethod
     def from_scipy_sparse(cls, x, /, *, fill_value=None):
         x = x.asformat("csr", copy=False)
+        if not x.has_canonical_format:
+            x.eliminate_zeros()
+            x.sum_duplicates()
+
         return cls((x.data, x.indices, x.indptr), shape=x.shape, fill_value=fill_value)
 
     def transpose(self, axes: None = None, copy: bool = False) -> Union["CSC", "CSR"]:
@@ -935,6 +942,10 @@ class CSC(_Compressed2d):
     @classmethod
     def from_scipy_sparse(cls, x, /, *, fill_value=None):
         x = x.asformat("csc", copy=False)
+        if not x.has_canonical_format:
+            x.eliminate_zeros()
+            x.sum_duplicates()
+
         return cls((x.data, x.indices, x.indptr), shape=x.shape, fill_value=fill_value)
 
     def transpose(self, axes: None = None, copy: bool = False) -> Union["CSC", "CSR"]:
diff --git a/sparse/numba_backend/_coo/core.py b/sparse/numba_backend/_coo/core.py
index 2db83fd..fe3b374 100644
--- a/sparse/numba_backend/_coo/core.py
+++ b/sparse/numba_backend/_coo/core.py
@@ -447,6 +447,10 @@ class COO(SparseArray, NDArrayOperatorsMixin):  # lgtm [py/missing-equals]
         True
         """
         x = x.asformat("coo")
+        if not x.has_canonical_format:
+            x.eliminate_zeros()
+            x.sum_duplicates()
+
         coords = np.empty((2, x.nnz), dtype=x.row.dtype)
         coords[0, :] = x.row
         coords[1, :] = x.col


Test patch:
diff --git a/sparse/numba_backend/tests/test_conversion.py b/sparse/numba_backend/tests/test_conversion.py
index 9f45ab8..7c343ce 100644
--- a/sparse/numba_backend/tests/test_conversion.py
+++ b/sparse/numba_backend/tests/test_conversion.py
@@ -3,6 +3,9 @@ from sparse.numba_backend._utils import assert_eq
 
 import pytest
 
+import numpy as np
+import scipy.sparse as sps
+
 FORMATS_ND = [
     sparse.COO,
     sparse.DOK,
@@ -38,3 +41,20 @@ def test_conversion_scalar(format1, format2):
     x = sparse.random((), format=format1, fill_value=0.5)
     y = x.asformat(format2)
     assert_eq(x, y)
+
+
+def test_non_canonical_conversion():
+    """
+    Regression test for gh-602.
+
+    Adapted from https://github.com/LiberTEM/sparseconverter/blob/4cfc0ee2ad4c37b07742db8f3643bcbd858a4e85/src/sparseconverter/__init__.py#L154-L183
+    """
+    data = np.array((2.0, 1.0, 3.0, 3.0, 1.0))
+    indices = np.array((1, 0, 0, 1, 1), dtype=int)
+    indptr = np.array((0, 2, 5), dtype=int)
+
+    x = sps.csr_matrix((data, indices, indptr), shape=(2, 2))
+    ref = np.array(((1.0, 2.0), (3.0, 4.0)))
+
+    gcxs_check = sparse.GCXS(x)
+    assert np.all(gcxs_check[:1].todense() == ref[:1]) and np.all(gcxs_check[1:].todense() == ref[1:])