Add hypothesis tests for data operators.

.gitignore

@@ -43,3 +43,5 @@ benchmark/benchmark_data.js
 *-tasks.txt
 *compiled_coeff*
 result_images/
+
+.hypothesis/

doc/changes/1957.misc

@@ -0,0 +1 @@
+Add hypothesis test for Data class operations (#1957).

qutip/core/data/base.pxd

@@ -10,7 +10,7 @@ cdef int idxint_DTYPE
 cdef class Data:
     cdef readonly (idxint, idxint) shape
     cpdef object to_array(self)
-    cpdef double complex trace(self)
+    cpdef double complex trace(self) except *
     cpdef Data adjoint(self)
     cpdef Data conj(self)
     cpdef Data transpose(self)

qutip/core/data/base.pyx

@@ -27,7 +27,7 @@ cdef class Data:
     cpdef object to_array(self):
         raise NotImplementedError
 
-    cpdef double complex trace(self):
+    cpdef double complex trace(self) except *:
         return NotImplementedError
 
     cpdef Data adjoint(self):

qutip/core/data/csr.pxd

@@ -23,7 +23,7 @@ cdef class CSR(base.Data):
     cpdef CSR copy(CSR self)
     cpdef object as_scipy(CSR self, bint full=*)
     cpdef CSR sort_indices(CSR self)
-    cpdef double complex trace(CSR self)
+    cpdef double complex trace(CSR self) except *
     cpdef CSR adjoint(CSR self)
     cpdef CSR conj(CSR self)
     cpdef CSR transpose(CSR self)

qutip/core/data/csr.pyx

@@ -229,7 +229,7 @@ cdef class CSR(base.Data):
             sort.inplace(self, ptr, diff)
         return self
 
-    cpdef double complex trace(self):
+    cpdef double complex trace(self) except *:
         return trace_csr(self)
 
     cpdef CSR adjoint(self):

qutip/core/data/dense.pxd

@@ -15,7 +15,7 @@ cdef class Dense(base.Data):
     cpdef Dense copy(Dense self)
     cpdef object as_ndarray(Dense self)
     cpdef object to_array(Dense self)
-    cpdef double complex trace(Dense self)
+    cpdef double complex trace(Dense self) except *
     cpdef Dense adjoint(Dense self)
     cpdef Dense conj(Dense self)
     cpdef Dense transpose(Dense self)

qutip/core/data/dense.pyx

@@ -192,7 +192,7 @@ cdef class Dense(base.Data):
         self._deallocate = False
         return self._np
 
-    cpdef double complex trace(self):
+    cpdef double complex trace(self) except *:
         return trace_dense(self)
 
     cpdef Dense adjoint(self):

qutip/core/data/mul.pyx

@@ -4,6 +4,8 @@
 from qutip.core.data cimport idxint, csr, CSR, dense, Dense, Data
 from scipy.linalg.cython_blas cimport zscal
 
+from libc.math cimport isfinite
+
 __all__ = [
     'mul', 'mul_csr', 'mul_dense',
     'imul', 'imul_csr', 'imul_dense', 'imul_data',
@@ -15,15 +17,24 @@ cpdef CSR imul_csr(CSR matrix, double complex value):
     """Multiply this CSR `matrix` by a complex scalar `value`."""
     cdef idxint l = csr.nnz(matrix)
     cdef int ONE=1
+    if not (isfinite(value.real) and isfinite(value.imag)):
+        raise ValueError(
+            f"Multiplying CSR matrix by non-finite value {value!r}"
+            " would generate a dense matrix of NaNs."
+        )
     zscal(&l, &value, matrix.data, &ONE)
     return matrix
 
 cpdef CSR mul_csr(CSR matrix, double complex value):
     """Multiply this CSR `matrix` by a complex scalar `value`."""
-    if value == 0:
-        return csr.zeros(matrix.shape[0], matrix.shape[1])
-    cdef CSR out = csr.copy_structure(matrix)
+    cdef CSR out
     cdef idxint ptr
+    if not (isfinite(value.real) and isfinite(value.imag)):
+        raise ValueError(
+            f"Multiplying CSR matrix by non-finite value {value!r}"
+            " would generate a dense matrix of NaNs."
+        )
+    out = csr.copy_structure(matrix)
     with nogil:
         for ptr in range(csr.nnz(matrix)):
             out.data[ptr] = value * matrix.data[ptr]

qutip/core/data/properties.pyx

@@ -10,6 +10,8 @@ from qutip.core.data.base cimport idxint
 from qutip.core.data cimport csr, dense, CSR, Dense
 from qutip.core.data.adjoint cimport transpose_csr
 
+from libc.math cimport isnan
+
 cdef extern from *:
     # Not defined in cpython.mem for some reason, but is in pymem.h.
     void *PyMem_Calloc(size_t nelem, size_t elsize)
@@ -201,23 +203,27 @@ cpdef bint isdiag_dense(Dense matrix) nogil:
 
 cpdef bint iszero_csr(CSR matrix, double tol=-1) nogil:
     cdef size_t ptr
+    cdef double complex v
     if tol < 0:
         with gil:
             tol = settings.core["atol"]
     tolsq = tol*tol
     for ptr in range(csr.nnz(matrix)):
-        if _abssq(matrix.data[ptr]) > tolsq:
+        v = matrix.data[ptr]
+        if _abssq(v) > tolsq or isnan(v.real) or isnan(v.imag):
             return False
     return True
 
 cpdef bint iszero_dense(Dense matrix, double tol=-1) nogil:
     cdef size_t ptr
+    cdef double complex v
     if tol < 0:
         with gil:
             tol = settings.core["atol"]
     tolsq = tol*tol
     for ptr in range(matrix.shape[0]*matrix.shape[1]):
-        if _abssq(matrix.data[ptr]) > tolsq:
+        v = matrix.data[ptr]
+        if _abssq(matrix.data[ptr]) > tolsq or isnan(v.real) or isnan(v.imag):
             return False
     return True
 

qutip/tests/core/data/test_operators_hypothesis.py

@@ -0,0 +1,195 @@
+import numpy
+
+from hypothesis import given, strategies as st
+
+import qutip
+from qutip.core import data as _data
+from qutip.tests import strategies as qst
+
+same_shape = st.shared(qst.qobj_shapes())
+matmul_shape_a, matmul_shape_b = qst.qobj_shared_shapes([
+    ("n", "k"),
+    ("k", "m"),
+])
+
+
+@given(qst.qobj_np(), qst.qobj_dtypes())
+def test_data_create(np_array, dtype):
+    data = _data.to(dtype, _data.create(np_array))
+    qst.note(data=data, np_array=np_array)
+    qst.assert_allclose(data.to_array(), np_array)
+
+
+@given(qst.qobj_datas())
+def test_data_neg_operator(data):
+    neg = -data
+    qst.note(neg=neg, data=data)
+    qst.assert_allclose(neg.to_array(), -data.to_array())
+
+
+@given(qst.qobj_datas())
+def test_data_iszero(data):
+    result = _data.iszero(data, tol=1e-15)
+    qst.note(result=result, data=data)
+    with qst.ignore_arithmetic_warnings():
+        np_array = data.to_array()
+        expected = numpy.allclose(
+            np_array, numpy.zeros_like(np_array), atol=1e-15, rtol=0,
+        )
+    assert result is expected
+
+
+@given(qst.qobj_datas(shape=same_shape), qst.qobj_datas(shape=same_shape))
+def test_data_add_operator(a, b):
+    result = a + b
+    qst.note(result=result, a=a, b=b)
+    with qst.ignore_arithmetic_warnings():
+        expected = a.to_array() + b.to_array()
+    qst.assert_allclose(result.to_array(), expected, treat_inf_as_nan=True)
+
+
+@given(qst.qobj_datas(shape=same_shape), qst.qobj_datas(shape=same_shape))
+def test_data_minus_operator(a, b):
+    result = a - b
+    qst.note(result=result, a=a, b=b)
+    with qst.ignore_arithmetic_warnings():
+        expected = a.to_array() - b.to_array()
+    qst.assert_allclose(result.to_array(), expected, treat_inf_as_nan=True)
+
+
+@given(st.complex_numbers(), qst.qobj_datas(shape=same_shape))
+@qst.raises_when(
+    ValueError,
+    "Multiplying CSR matrix by non-finite value {x!r}"
+    " would generate a dense matrix of NaNs.",
+    when=lambda a, x: not numpy.isfinite(x) and isinstance(a, _data.CSR),
+)
+def test_data_scalar_multiplication_left_operator(x, a):
+    result = x * a
+    qst.note(result=result, x=x, a=a)
+    with qst.ignore_arithmetic_warnings():
+        expected = x * a.to_array()
+    qst.assert_allclose(result.to_array(), expected, treat_inf_as_nan=True)
+
+
+@given(qst.qobj_datas(shape=same_shape), st.complex_numbers())
+@qst.raises_when(
+    ValueError,
+    "Multiplying CSR matrix by non-finite value {x!r}"
+    " would generate a dense matrix of NaNs.",
+    when=lambda a, x: not numpy.isfinite(x) and isinstance(a, _data.CSR),
+)
+def test_data_scalar_multiplication_right_operator(a, x):
+    result = a * x
+    qst.note(result=result, a=a, x=x)
+    with qst.ignore_arithmetic_warnings():
+        expected = a.to_array() * x
+    qst.assert_allclose(result.to_array(), expected, treat_inf_as_nan=True)
+
+
+@given(
+    qst.qobj_datas(shape=same_shape),
+    st.complex_numbers(min_magnitude=1e-12),
+)
+@qst.raises_when(
+    ValueError,
+    "Multiplying CSR matrix by non-finite value {inv_x!r}"
+    " would generate a dense matrix of NaNs.",
+    when=lambda a, x: not numpy.isfinite(1 / x) and isinstance(a, _data.CSR),
+    format_args=lambda a, x: {"inv_x": 1 / x},
+)
+def test_data_scalar_division_operator(a, x):
+    result = a / x
+    qst.note(result=result, a=a, x=x)
+    with qst.ignore_arithmetic_warnings():
+        expected = a.to_array() * (1 / x)
+    qst.assert_allclose(result.to_array(), expected, treat_inf_as_nan=True)
+
+
+@given(qst.qobj_datas(shape=same_shape), qst.qobj_datas(shape=same_shape))
+def test_data_equality_operator_same_shapes(a, b):
+    with qutip.CoreOptions(atol=1e-15):
+        result = (a == b)
+    qst.note(result=result, a=a, b=b)
+    with qst.ignore_arithmetic_warnings():
+        expected = numpy.allclose(
+            numpy.zeros(a.shape),
+            a.to_array() - b.to_array(),
+            atol=1e-15, rtol=0,
+        )
+    assert result == expected
+
+
+@given(qst.qobj_datas(), qst.qobj_datas())
+def test_data_equality_operator_different_shapes(a, b):
+    with qutip.CoreOptions(atol=1e-15):
+        result = (a == b)
+    qst.note(result=result, a=a, b=b)
+    if a.shape != b.shape:
+        assert result is False
+    else:
+        with qst.ignore_arithmetic_warnings():
+            expected = numpy.allclose(
+                numpy.zeros(a.shape),
+                a.to_array() - b.to_array(),
+                atol=1e-15
+            )
+        assert result == expected
+
+
+@given(
+    qst.qobj_datas_matmul(shape=matmul_shape_a),
+    qst.qobj_datas_matmul(shape=matmul_shape_b),
+)
+def test_data_matmul_operator(a, b):
+    result = a @ b
+    qst.note(result=result, a=a, b=b)
+    with qst.ignore_arithmetic_warnings():
+        expected = a.to_array() @ b.to_array()
+    qst.assert_allclose(
+        result.to_array(), expected,
+        atol=1e-10, rtol=1e-10, treat_inf_as_nan=True,
+    )
+
+
+@given(qst.qobj_datas())
+@qst.raises_when(
+    ValueError,
+    "matrix shape {shape} is not square.",
+    when=lambda data: data.shape[0] != data.shape[1],
+    format_args=lambda data: {"shape": data.shape},
+)
+def test_trace(data):
+    result = data.trace()
+    with qst.ignore_arithmetic_warnings():
+        expected = data.to_array().trace()
+    qst.note(result=result, data=data)
+    qst.assert_allclose(result, expected)
+
+
+@given(qst.qobj_datas())
+def test_adjoint(data):
+    result = data.adjoint()
+    qst.note(result=result, data=data)
+    qst.assert_allclose(result.to_array(), data.to_array().T.conj())
+
+
+@given(qst.qobj_datas())
+def test_conj(data):
+    result = data.conj()
+    qst.note(result=result, data=data)
+    qst.assert_allclose(result.to_array(), data.to_array().conj())
+
+
+@given(qst.qobj_datas())
+def test_transpose(data):
+    result = data.transpose()
+    qst.note(result=result, data=data)
+    qst.assert_allclose(result.to_array(), data.to_array().T)
+
+
+@given(qst.qobj_datas())
+def test_copy(data):
+    result = data.copy()
+    qst.note(result=result, data=data)
+    qst.assert_allclose(result.to_array(), data.to_array())