qutip: QuTiP incompatible with numpy 1.20 due to Qobj.__array__

QuTiP is currently incompatible with numpy 1.20. Tests will fail to even collect with an error such as

AttributeError: 'numpy.ndarray' object has no attribute 'dag'

and a lot of functionality will break - anything that requires Qobj.eigenstates() for example.

This is because Qobj defines __array__, one of numpy’s “array interface” functions, intended for classes that can be safely converted implicitly into an ndarray. This isn’t really the case for Qobj - it loses all sorts of information when you do that, which is why we’ve maintained the separate Qobj.full() for explicitly getting the dense matrix representation of a Qobj. This is not to mention that numpy ufuncs probably should not be able to implicitly convert Qobj - I’d strongly argue that np.sin(qutip.basis(2, 1)) should be TypeError, not array([[0. +0.j], [0.84147098+0.j]]) (like it is right now). We actually already removed this “functionality” in dev.major.

The most pressing incompatibility is that a few points in QuTiP put a few Qobj into a np.array(dtype=object). In numpy 1.20, this no longer produces a 1D array of Qobj, but a 3D array of complex. This breaks Qobj.eigenstates, and prevents test collection due to it being present in states.py::qutrit_basis(), which is called during parametrisation.

There are a two possible ways to solve this, and we ought to release a fix with one of them in a patch ASAP:

1. remove Qobj.__array__
2. remove all use of Qobj in np.array

I’m personally in favour of “explicit is better than implicit” in this case, i.e. removing Qobj.__array__ and relying on Qobj.full(). As another example along this vein, note scipy.sparse matrices don’t implement this either, and they’re arguably closer to being safely coerced to ndarray than we are.

Related issues

#938: feature request for implementing __array__. This only asks for np.array(qobj) as a convenience, acknowledging the availability of Qobj.full(). #1017: includes a comment on buggy behaviour caused by __array__. Note that the solution given there (np.asarray(..., dtype=object)) will no longer work with numpy 1.20.