Repository: scipy/scipy
Pull request: BUG: sparse: Propagate dtype through CSR/CSC constructors
#### Reference issue
Fixes gh-13329

#### What does this implement/fix?
Propagates the `dtype` parameter through the intermediate COO format constructor.

Problem statement:
Large sparse matrices of big integers lose information
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

When I create a pretty large sparse matrix of big integers, somehow I get slightly different values stored in the matrix, even for values that should fit in an `uint64`. I suspect this is because there's a convertion from int to float and back to int somewhere.

#### Reproducing code example:
<!-- 
If you place your code between the triple backticks below, 
it will be rendered as a code block. 
-->

```python
import random

import numpy as np
import scipy.sparse

size = 100_000
m = n = 10_000
data = [random.randrange(2**64 - 1) for _ in range(size)]
row_ind = [i * m // size for i in range(size)]
col_ind = [i * size // m % n for i in range(size)]
assert len(set(zip(row_ind, col_ind))) == size
assert max(row_ind) < m
assert max(col_ind) < n

matrix = scipy.sparse.csr_matrix((data, (row_ind, col_ind)), dtype=np.uint64)

for i in range(size):
    v = matrix[row_ind[i], col_ind[i]]
    assert v == data[i], (i, data[i], v)
```

#### Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
AssertionError: (2, 18072150216881262918, 18072150216881262592)
```

#### Scipy/Numpy/Python version information:
```
1.4.1 1.18.3 sys.version_info(major=3, minor=8, micro=6, releaselevel='final', serial=0)
```

Hints from issue discussion:
When you create your data like that, you will almost certainly have some repeated occurrences of (row, col) pairs, and when you create the sparse matrix, the values associated with the repeated coordinates will be added. This could result in an overflow of the 64 bit unsigned integers.

Actually, the problem isn't necessarily overflow (although that might happen).  It is just that the repeated occurrences are summed, so you can't directly compare the contents of `matrix` with the original values in `data`.
Thanks Warren, I made a mistake in writing the code I needed to reproduce the bug. Now I patched the example to make sure that the `(row, col)` pairs are all distinct (see the first assertion): I still get an error, where the number I get changes in the 3 least significant digits.

Apparently I _don't_ get an error when `data` is a numpy array rather than a list.
@matteodellamico, thanks for the update.  This looks like a bug.  Even though you specified `dtype=np.uint64`, an intermediate step converts the values in `data` to floating point.  That intermediate result is then cast back to `np.uint64`.  Sufficiently large values in `data` will be changed in the round trip through floating point.

For other devs, here is what is happening:

`csr_matrix` is a subclass of `_cs_matrix`, and in `_cs_matrix.__init__`, when the input is a tuple of the form `(data, ij)`, the argument is passed to `coo_matrix` without specifying the dtype.  In the `coo_matrix` code, `data` eventually gets passed to `np.array` without a dtype being specified.  In the code that shows the problem, `data` is a Python list of Python integers, with some values greater than or equal to `2**63` and some less than that. This is one of those cases where `np.array` surprisingly returns a floating point array, even though all the values can be represented as `uint64`.  A simple example is
```
In [30]: np.array([2**63, 1])                                                                                                        
Out[30]: array([9.22337204e+18, 1.00000000e+00])
```

> Apparently I don't get an error when data is a numpy array rather than a list.

In that case, `data` is already a NumPy array, so `np.array(data)` returns `data`.  This avoids the round trip through `np.float64`. 
@matteodellamico  Has the original reproducer above been edited to correct for the unrelated problem?
Yes, it has.

Commit message: Fix the regression tests, and the underlying bug

Patch to review:
diff --git a/scipy/sparse/compressed.py b/scipy/sparse/compressed.py
index 17b9a23ce7e5..a9e6d114fdd2 100644
--- a/scipy/sparse/compressed.py
+++ b/scipy/sparse/compressed.py
@@ -51,7 +51,7 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
                 if len(arg1) == 2:
                     # (data, ij) format
                     from .coo import coo_matrix
-                    other = self.__class__(coo_matrix(arg1, shape=shape))
+                    other = self.__class__(coo_matrix(arg1, shape=shape, dtype=dtype))
                     self._set_self(other)
                 elif len(arg1) == 3:
                     # (data, indices, indptr) format
diff --git a/scipy/sparse/coo.py b/scipy/sparse/coo.py
index c53a2edf307a..1b6f078ec8e1 100644
--- a/scipy/sparse/coo.py
+++ b/scipy/sparse/coo.py
@@ -131,9 +131,10 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
                 M, N = arg1
                 self._shape = check_shape((M, N))
                 idx_dtype = get_index_dtype(maxval=max(M, N))
+                data_dtype = getdtype(dtype, default=float)
                 self.row = np.array([], dtype=idx_dtype)
                 self.col = np.array([], dtype=idx_dtype)
-                self.data = np.array([], getdtype(dtype, default=float))
+                self.data = np.array([], dtype=data_dtype)
                 self.has_canonical_format = True
             else:
                 try:
@@ -154,11 +155,11 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
                     self._shape = check_shape((M, N))
 
                 idx_dtype = get_index_dtype(maxval=max(self.shape))
+                data_dtype = getdtype(dtype, obj, default=float)
                 self.row = np.array(row, copy=copy, dtype=idx_dtype)
                 self.col = np.array(col, copy=copy, dtype=idx_dtype)
-                self.data = np.array(obj, copy=copy)
+                self.data = np.array(obj, copy=copy, dtype=data_dtype)
                 self.has_canonical_format = False
-
         else:
             if isspmatrix(arg1):
                 if isspmatrix_coo(arg1) and copy:
diff --git a/scipy/sparse/tests/test_base.py b/scipy/sparse/tests/test_base.py
index c720c176abaa..ddcf9e3c2b58 100644
--- a/scipy/sparse/tests/test_base.py
+++ b/scipy/sparse/tests/test_base.py
@@ -3523,6 +3523,11 @@ def test_constructor4(self):
         ij = vstack((row,col))
         csr = csr_matrix((data,ij),(4,3))
         assert_array_equal(arange(12).reshape(4,3),csr.todense())
+        
+        # using Python lists and a specified dtype
+        csr = csr_matrix(([2**63 + 1, 1], ([0, 1], [0, 1])), dtype=np.uint64)
+        dense = array([[2**63 + 1, 0], [0, 1]], dtype=np.uint64)
+        assert_array_equal(dense, csr.toarray())
 
     def test_constructor5(self):
         # infer dimensions from arrays
@@ -4086,13 +4091,16 @@ def test_constructor1(self):
         # unsorted triplet format
         row = array([2, 3, 1, 3, 0, 1, 3, 0, 2, 1, 2])
         col = array([0, 1, 0, 0, 1, 1, 2, 2, 2, 2, 1])
-        data = array([6., 10., 3., 9., 1., 4.,
-                              11., 2., 8., 5., 7.])
+        data = array([6., 10., 3., 9., 1., 4., 11., 2., 8., 5., 7.])
 
         coo = coo_matrix((data,(row,col)),(4,3))
-
         assert_array_equal(arange(12).reshape(4,3),coo.todense())
 
+        # using Python lists and a specified dtype
+        coo = coo_matrix(([2**63 + 1, 1], ([0, 1], [0, 1])), dtype=np.uint64)
+        dense = array([[2**63 + 1, 0], [0, 1]], dtype=np.uint64)
+        assert_array_equal(dense, coo.toarray())
+
     def test_constructor2(self):
         # unsorted triplet format with duplicates (which are summed)
         row = array([0,1,2,2,2,2,0,0,2,2])


Code review comments:
diff_hunk: @@ -51,7 +51,7 @@ def __init__(self, arg1, shape=None, dtype=None, copy=False):
                 if len(arg1) == 2:
                     # (data, ij) format
                     from .coo import coo_matrix
-                    other = self.__class__(coo_matrix(arg1, shape=shape))
+                    other = self.__class__(coo_matrix(arg1, shape=shape, dtype=dtype))
path: scipy/sparse/compressed.py
text: @user1:
PEP-8 line length.
```suggestion
                    other = self.__class__(coo_matrix(arg1, shape=shape,
                                                      dtype=dtype))
```