Changed __mul__ behviour to be more flexible and consistent.

[AGaliciaMartinez]

Description
I changed mul behaviour following the discussion of issue #1607. However, the implementation details differ slightly from what was discussed there. The current behaviout of __mul__(self, other) is:

• If other is a Qobj: dispatch to __matmul__
• If not try dipatching to mul and return NotImplemented if TypeError is returned.
• The dispatcher now gets other instead of complex(other). This is more flexible and allows specialisations to handle arbitrary scalar like objects (something extremely useful for qutip-tensorflow).
• Infer hermiticity and unitarity when possible trying complex(other)

__rmul__ now directly dispatches to __mul__. Any necessary check (casting to complex included) is done in __mul__.

No changes were required to mul_dense and mul_csr as add_dense(data, value=np.array(1)) works. It internally tries complex(np.array(1)) which is guaranteed to work. This is something I was very happy to find as I do not think that specialisations should return NotImplemented (although they still can if required).

Related issues or PRs
Fixes issue #1611

Changelog

• Qobj __mul__ now handles consistently right and left multiplications of an arbitrary python object.
• __mul__ is now more flexible passing other to the dispatcher instead of complex(other).
• qobj*np.array([1,2]) and qobj*np.array([1,2]) (or any other numpy array that does not represent an scalar) now raise TypeError. This change is not backwards compatible (!!).

Edit: changed changelog and description

[hodgestar]

@AGaliciaMartinez Perhaps it would be better to check for whether other is complex directly using isinstance(other, numbers.Number). This is the canonical approach for detecting scalar numbers and avoids raising and catching an exception and allows objects that support just __complex__ to still be dispatched for special handling.

[AGaliciaMartinez]

There are two reasons I did not go with isinstance(other, numbers.Number):

• This returns False when other is a NumPy array or TensorFlow tensor which makes __mul__ significantly less flexible.
• If we accept objects that are not instances of numbers.Number, we will still have to use complex(other). This is necessary to infer if the output is hermitian or unitary. Hence, instead of doing a set of isinstance checks and shape checks for NumPy TensorFlow or any other future library we may want to be somewhat compatible with, we just check whether other it can be casted to complex. If it can, mul makes sense and we let the specialisation do whatever they "want" (they could always just do complex(value) in any case).

I must admit that the try except with complex also seems not ideal to me. But I could not find a better implementation that fitted the requirements. I mean, we do use this same code pattern at some other points in the code (in __matmul__ we raise not implemented if other can not be casted to Qobj and we catch this with a try) but what worries me is: what happens if complex(other) raises TypeError not because we do not know how to do the cast but because the user made a mistake in the code? An example of this could be:

arrray = np.array([1,2])  # length 2 array
qobj*array

This will return TypeError but with a different message if we catch it in the complex(other) try execpt.

• With a try execpt that returns NotImplemented (current implementation) the error message is: TypeError: operand 'Qobj' does not support ufuncs (__array_ufunc__=None)
• Wehreas without the try expect (not returnt NotImplemented if complex(other) fails), the error message is: TypeError: only length-1 arrays can be converted to Python scalars, which could be a useful error message.

[coveralls]

Coverage decreased (-0.7%) to 64.838% when pulling eccae25 on AGaliciaMartinez:mul into f3dd504 on qutip:dev.major.

[Ericgig]

We don't want the Qobj method to set what is a scalar but let the data-layer do it. So I don't think
complex(other) or isinstance(other, numbers.Number) should be used as the condition to return NotImplimented.

Setting isherm and isunitary is optional, we can set them to None is we can't read other.

[AGaliciaMartinez]

In my previous GSoC meeting we discussed what to do with the changes proposed in this PR. Currently (in this PR) we return NotImplemented if complex(other) raises a TypeError. The major concerns here are that:

• This is not the recommended way to return NotImplemented. It is preferred to do a set of isinstace checks instead.
• This defines what an scalar is at he Qobj level.

It seems that the discussion steers to what a scalar is (which classes are scalars) and at what level do we define it (who is responsible for returning NotImplemented).

Proposed solutions:

1. Let the specialisation return NotImplemented and set isherm only when complex(other) is possible. This defines what an scalar is at the specialisation level. However:

   • It may cause inconsistent behaviour where some operations are allowed for a data layer (TfTensor or any other) but not for another one. I would argue that all specialisations should accept the same input, although operations may be handled differently internally (for example, in tensorflow we want to keep the graph state).
   • NotImplemented only makes sense for some special binary methods (__add__, __lt__, etc) as Python employs NotImplemented to try other's binary method. My concern is that if the dispatcher is used by itself at some points of the code (such as for the solver) returning NotImplemented may cause wrong error messages or even unexpected behaviour with future specialisations ( for instance, although it now raises a deprecation warning, NotImplemented is considered as True when evaluated as a condition). It seems to me that we should avoid having NotImplemented in python space (except for arithmetic operations in a class where they are useful and always make sense).

2. Simon suggested to define what a scalar is at the dispatcher level (or something similar to this but correct me if I am wrong). This would require to have a function that casts any defined scalar to a complex scalar. The dispatcher would be responsible to ensure that the correct scalar gets passed appropriately to the specialisation and for returning not implemented if necessary.

   • This solves the inconsistency problem as you could register a "new" scalar and define how it gets cast to a python complex scalar with a function (in most of the times, if not always, this will be complex(other)).
   • This does not solve having NotImplemented at the python space.
   • It seems to be challenging to implement and will probably make things more difficult to maintain.

3. We can also define more rigorously what a scalar is for Qobj with a list, _ALLOWED_SCALARS = [numbers.Number, ...] that defines what we accept as scalars. We then do: isinstace(other, _ALLOWED_SCALARS) and return NotImplemented if not there. Two options here:

   1. QuTiP allows only numbers.Number:
      • In this case qutip-tensorflow can extend _ALLOWED_SCALARS with tf.Variable and tf.Tensor. However, if the user uses two plug-ins, say qutip-tensorflow and qutip-something-that-does-not-work-well-with-tensorflow (second plug-in), the second plug-in may not work properly.
      • The above point is not a big deal if we make clear that _ALLOWED_SCALARS must accept complex(other) for compatibility reasons. In this way the second plug-in can always default to complex(other) (this is already the default behaviour of mul_dense).
   2. QuTiP allows numbers.Number, ndarray of shape (,), tf.Variable of shape (,) and more in the future.
      • The upside of this method with respect to 3.1 is that users know which types they must accept.
      • This is probably not a good idea as it means qutip has to handle new dependencies (tensorflow and probalby more in the future).
   • This still defines what an scalar is at he Qobj level. I think this is not bad as it seems to be the only way of avoiding NotImplemented at Python space exept for in special methods (__mul__, ...).

I think that 3.i is the closest to an ideal solution.

[jakelishman]

This all seems incredibly complicated. How about all specialisations stored in the dispatcher take object as the scalar type, and throw a TypeError if they can't handle it. Qobj.__mul__ just does a try/except, and returns NotImplemented if it catches a TypeError. If it doesn't catch a type error, it calls a separate function, something like

@functools.singledispatch
def _scalar_properties(x):
    return {
        # if the matrix was Hermitian, is the output Hermitian (True), or not known (False)
        'preserves_hermicity': bool,
        # regardless of the input matrix, is the output Hermitian (True), or not known (False)
        # (i.e. was the scalar zero-like).
        'forces_hermicity': bool,
        # if the matrix was unitary, is it still unitary (True), or not known (False)
        'preserves_unitarity': Optional[bool],
    }

Then just register types with the functools dispatcher?

So then Qobj.__mul__ looks like:

def __mul__(self, other):
    if isinstance(other, (Qobj, QobjEvo)):
        return self.__matmul__(other)
    try:
        out = _data.mul(self.data, other)
    except TypeError:
        return NotImplemented
    properties = _scalar_properties(other)
    isherm = properties['forces_hermicity'] or (self.isherm and properties['preserves_hermicity'] or None)
    isunitary = self.isunitary and properties['preserves_unitarity'] or None
    return Qobj(
        out, copy=False, dims=self.dims, isherm=isherm, isunitary=isunitary,
    )

(I've probably not got the property truthiness checking exactly right, but that's the principle.)

[jakelishman]

If you really wanted, you could even have it so that the _scalar_properties function was an attribute of the specialisation you gave to the dispatcher, then you wouldn't even need the functools single dispatcher, and you could be sure that the specialising was interpreting the scalar in the exact same way when it calculated the properties:

callable = _data.mul.get(self.data, other)
out = callable(self.data, other)
properties = callable.base_operation.scalar_properties(other)

(the _constructed_specialisation.base_operation property doesn't currently exist, but it's pretty trivial to add).

The properties in my example really correspond to the booleans is_real, is_zero and is_unit_norm, but the names aren't super important. The point is that you don't actually care about the number having a solid, known numerical representation. You only care about certain properties of the number, and they may be known without knowing a fixed value. So it makes more sense to query the properties, than to try and query a number, when the latter question doesn't make sense for symbolic quantities.

[AGaliciaMartinez]

I am actually fine with using a try/except with data.mul (definitely better than try/except with complex(other)). Is just that from previous feedback I thought this was not the canonical approach?

Regarding to _scalar_properties, I think it is a very interesting idea that would definitely be helpful for symbolic representations of Qobj. For current extensions (qutip-tensorflow and qutip-cupy) it is not really needed as complex(other) works just fine. I would prefer to not implement it and leave this feature for another PR.

[jakelishman]

The upshot is that magic operator methods should never raise for invalid types, only return NotImplemented (ideally as fast as possible). However, the data-layer dispatched functions aren't only used in Qobj operator methods; they're also used in internal QuTiP functions, where the correct semantics for passing an invalid type is to raise TypeError. So from that perspective, I'd say that a specialisation should raise TypeError, and Qobj.__mul__ and co should convert that into return NotImplemented.

I don't remember saying that the dispatched specialisations shouldn't raise at all - if I did, there may be more considerations that I can't think of right now, but I don't think it's necessarily the case; they already raise on badly shaped inputs, or things like that.

[AGaliciaMartinez]

I After the discussion in the GSoC meeting I changed the implementation to:

• try mul and return NotImplemented if TypeError is raised. (as jake suggested)
• try complex(other) to inger ishemr and isunitary. (as Eric suggested)

I think this is the simplest solution. @jakelishman are you ok with this implementation? I think the ideas of numpy broadcasting Qobj as scalar and _scalar_properties are really good but probably out of the scope of this PR.

[jakelishman]

Fine by me. You lose some additional cases where you would know that the multiplication maintains hermicity or whatever, but that can easily be handled when it's actually needed.

[hodgestar]

Left two small suggestions.

Do we have any tests for whether the isherm and isunitary preserving logic is correct? If not, could we add them in this PR? I know it's not related directly to the fix, but we are changing the code and we might as well add tests while we're looking at it and know what it is supposed to do.

[hodgestar]

I don't know if the test failure is related or not -- it looks unrelated but it is a data layer test, so we should probably check a bit carefully.

[hodgestar]

@AGaliciaMartinez Could you also update the changelog entry in the PR description to say something about the change?

[AGaliciaMartinez]

Do we have any tests for whether the isherm and isunitary preserving logic is correct?

We do have a few tests, see for instance test_QobjUnitaryOper which I modified in this PR to track a few more cases. However, the current tests are not very exhaustive. For example, they do not really ensure that hermiticity is inferred in the operation. some rethinking of these tests may be beneficial but I would prefer leaving that for another PR. Should I open an issue to remember this or do you prefer if I extend these tests in this PR?

[hodgestar]

We do have a few tests, see for instance test_QobjUnitaryOper which I modified in this PR to track a few more cases.

Perfect. Thank you!

However, the current tests are not very exhaustive. For example, they do not really ensure that hermiticity is inferred in the operation. some rethinking of these tests may be beneficial but I would prefer leaving that for another PR. Should I open an issue to remember this or do you prefer if I extend these tests in this PR?

Happy to leave this for another PR.