TestDicke fails with scipy 1.6.1: TypeError: can't convert complex to float

[drew-parsons]

Some TestDicke tests in test_piqs.py (test_lindbladian, test_lindbladian_dims, test_liouvillian) fail with the recent scipy 1.6.1 release. They were previously passing with scipy 1.6.0.

To Reproduce

$  cp -r qutip/tests/  /tmp/qutip
$  cd /tmp/qutip
$  pytest-3 -v -k "TestDicke"

The test error message from TestDicke.test_lindbladian is

___________________________________________________________________________________________ TestDicke.test_lindbladian ____________________________________________________________________________________________

self = <tests.test_piqs.TestDicke object at 0x7f55475a0c10>

    def test_lindbladian(self):
        """
        PIQS: Test the generation of the Lindbladian matrix.
        """
        N = 1
        gCE = 0.5
        gCD = 0.5
        gCP = 0.5
        gE = 0.1
        gD = 0.1
        gP = 0.1
    
        system = Dicke(
            N=N,
            emission=gE,
            pumping=gP,
            dephasing=gD,
            collective_emission=gCE,
            collective_pumping=gCP,
            collective_dephasing=gCD,
        )
    
>       lindbladian = system.lindbladian()

tests/test_piqs.py:450: 
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
/usr/lib/python3/dist-packages/qutip/piqs.py:509: in lindbladian
    return cythonized_dicke.lindbladian()
qutip/cy/piqs.pyx:313: in qutip.cy.piqs.Dicke.lindbladian
    ???
qutip/cy/piqs.pyx:431: in qutip.cy.piqs.Dicke.lindbladian
    ???
/usr/lib/python3/dist-packages/scipy/sparse/compressed.py:54: in __init__
    other = self.__class__(coo_matrix(arg1, shape=shape,
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

self = <[AttributeError('dtype not found') raised in repr()] coo_matrix object at 0x7f55475a0ca0>
arg1 = ([(-0.6000000238418579+0j), (0.6000000238418579+0j), (-0.9000000357627869+0j), (-0.9000000357627869+0j), (-0.6000000238418579+0j), (0.6000000014901161+0j)], ([3, 3, 2, 1, 0, 0], [3, 0, 2, 1, 0, 3]))
shape = (4, 4), dtype = None, copy = False

    def __init__(self, arg1, shape=None, dtype=None, copy=False):
        _data_matrix.__init__(self)
    
        if isinstance(arg1, tuple):
            if isshape(arg1):
                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([], dtype=data_dtype)
                self.has_canonical_format = True
            else:
                try:
                    obj, (row, col) = arg1
                except (TypeError, ValueError) as e:
                    raise TypeError('invalid input format') from e
    
                if shape is None:
                    if len(row) == 0 or len(col) == 0:
                        raise ValueError('cannot infer dimensions from zero '
                                         'sized index arrays')
                    M = operator.index(np.max(row)) + 1
                    N = operator.index(np.max(col)) + 1
                    self._shape = check_shape((M, N))
                else:
                    # Use 2 steps to ensure shape has length 2.
                    M, N = shape
                    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, dtype=data_dtype)
E               TypeError: can't convert complex to float

/usr/lib/python3/dist-packages/scipy/sparse/coo.py:161: TypeError

Likewise for TestDicke.test_lindbladian_dims and TestDicke.test_liouvillian.

Your Environment

qutip 4.5.3 built on Debian unstable.

qutip.about()

QuTiP: Quantum Toolbox in Python
================================
Copyright (c) QuTiP team 2011 and later.
Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li and Jake Lishman.
Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng.
Original developers: R. J. Johansson & P. D. Nation.
Previous lead developers: Chris Granade & A. Grimsmo.
Currently developed through wide collaboration. See https://github.com/qutip for details.

QuTiP Version:      4.5.3
Numpy Version:      1.19.5
Scipy Version:      1.6.1
Cython Version:     0.29.21
Matplotlib Version: 3.3.4
Python Version:     3.9.1
Number of CPUs:     4
BLAS Info:          OPENBLAS
OPENMP Installed:   True
INTEL MKL Ext:      False
Platform Info:      Linux (x86_64)
Installation path:  /usr/lib/python3/dist-packages/qutip
================================================================================
Please cite QuTiP in your publication.
================================================================================
For your convenience a bibtex reference can be easily generated using `qutip.cite()`

[drew-parsons]

It's probably relevant that scipy 1.6.1 fixed some problems with sparse matrices (with COO format constructor), see https://docs.scipy.org/doc/scipy-1.6.1/reference/release.1.6.1.html including PR#13403 scipy/scipy#13403

[jakelishman]

Thanks for the bug report and the detailed look into it! As a temporary work-around, in file qutip/cy/piqs.pyx, change lines 431 to 433

qutip/qutip/cy/piqs.pyx

Lines 431 to 433 in 2aaae75

	cdef lindblad_matrix = csr_matrix((lindblad_data,
	(lindblad_row, lindblad_col)),
	shape=(nds**2, nds**2))

to

        cdef lindblad_matrix = csr_matrix((lindblad_data,
                                          (lindblad_row, lindblad_col)),
                                          shape=(nds**2, nds**2),
                                          dtype=np.complex128)

and recompile. If you want to make a PR of something similar against QuTiP, I'll accept it.

I would actually file this against scipy.sparse - I think our usage is completely in line with the contract of scipy.sparse.csr_matrix and they've got a bug in their dtype handling. You currently can't construct a CSR matrix using the COO triplet format for complex data, unless the dtype is made explicit somewhere, but the constructor is meant to correctly infer a suitable dtype if one is not passed.

Basic Scipy reproducer to illustrate the problem:

>>> import scipy.sparse
>>> scipy.__version__
'1.6.1'
>>> scipy.sparse.csr_matrix(([1+1j], ([0], [0])), shape=(2, 2))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/jake/.anaconda3/envs/qutip-dev/lib/python3.8/site-packages/scipy/sparse/compressed.py", line 54, in __init__
    other = self.__class__(coo_matrix(arg1, shape=shape,
  File "/Users/jake/.anaconda3/envs/qutip-dev/lib/python3.8/site-packages/scipy/sparse/coo.py", line 161, in __init__
    self.data = np.array(obj, copy=copy, dtype=data_dtype)
TypeError: can't convert complex to float

This can be fixed either by passing dtype=np.complex128 to the constructor, or passing the data inside a NumPy array, since that'll also fix the dtype.

I imagine with their implicit conversions, SciPy may also need to test the special cases where all list elements are things like 1+0j, which have type complex but can be safely represented by reals - the Python call float(1 + 0j) is forbidden even though the imaginary part is 0. I actually originally thought this was the problem in this issue, since all the Lindbladian data tested is real numbers with complex type. I suppose it's up to SciPy to decide how they want to handle that case - either always maintaining complex or putting in a special-case cast (np.real(x)) for known-safe complex -> float conversions.

[jakelishman]

Note: they're already talking about a change in scipy/scipy#13585 and scipy/scipy#13586. I'll post this repro on the bottom of that.

[drew-parsons]

Thanks for the speedy response @jakelishman . Looks like scipy 1.6.1 is not quite as stable as we would have liked it to have been. I'll wait for the scipy dust to settle.

[jakelishman]

We can still merge the workaround into qutip if you'd like to make the PR (otherwise I'll take care of it). There's no harm in being explicit since we know we will always want complex dtype.

[drew-parsons]

Agreed, safer to make it explicit. I'll prepare a PR.

[drew-parsons]

Hi @jakelishman , cqobjevo.pyx also uses csr_matrix, already with dtype=complex. Is it correct to use dtype=complex128 in piqs.pyx and dtype=complex in cqobjevo.pyx, or should the two cases use the same type?

[jakelishman]

Numpy interprets the Python base type complex as equal to np.complex128 when passed as a dtype. Personally I think np.complex128 is much clearer (since it specifies the size in the name too), but I wouldn't worry too much about changing everything everywhere. cqobjevo.pyx in particular is going to be nearly entirely rewritten in the next major QuTiP release anyway.

[drew-parsons]

No worries, I'll stick with just piq.pyx then.