Circuit.unitary breaks in the presente of two-qubit gate initialized by Gate.controlled_by

[renatomello]

Code to reproduce bug:

from qibo import Circuit, gates

circuit = Circuit(2)
circuit.add(gates.H(0).controlled_by(1))

circuit.unitary()

[rnhmjoj]

gates.Unitary has the same problem, It seems the matrix representation is not updated after calling controlled_by:

>>> q.gates.X(1).controlled_by(0).matrix()
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j]])
>>> q.gates.Unitary(np.array([[0,1],[1,0]])).controlled_by(0).matrix()
array([[0.+0.j, 1.+0.j],
       [1.+0.j, 0.+0.j]])

[rnhmjoj]

To investigate a bit I made a custom gate by subclassing gates.Gate and implemented the .matrix() method myself.
I thus discovered what seems to be an inconsistency with controlled gates in how matrix() is interpreted.
When using the numpy backed it seems that .matrix_fused() calls .matrix() expecting the full 2ⁿ×2ⁿ controlled matrix representation, while .apply_gate() wants the 2×2 matrix.

With this hack I can work around the issue consistently:

class FGate(q.gates.Gate):
    '''
    The F(α, β) gate is an multi-controlled gate defined by:

         F(α,β) |1..10> = α|1..10>
         F(α,β) |1..11> = β|1..10>

    Its matrix representation is diag(1,..,1,α,β).
    '''
    def __init__(self, diag, target, name='F'):
        super().__init__()
        self._values = diag
        self.unitary = True
        self.target_qubits = (target,)
        self.name = 'fgate'
        self.draw_label = name
        controls = range(target)
        if controls:
            self.control_qubits = tuple(controls)
            self.is_controlled_by = True

    def matrix(self, backend=None):
        import inspect
        caller = inspect.stack()[1]

        if caller.function == 'matrix_fused':
            npad = len(self._values) * (2 ** len(self.control_qubits) - 1)
            return np.diag(np.concatenate([np.ones(npad), self._values]))
        elif caller.function == 'apply_gate':
            return np.diag(self._values)

[renatomello]

gates.Unitary has the same problem, It seems the matrix representation is not updated after calling controlled_by:

>>> q.gates.X(1).controlled_by(0).matrix()
array([[1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j],
       [0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j]])
>>> q.gates.Unitary(np.array([[0,1],[1,0]])).controlled_by(0).matrix()
array([[0.+0.j, 1.+0.j],
       [1.+0.j, 0.+0.j]])

This inconsistency was due to these classes defaulting to other gate classes depending on the number of controls added. In your example, X controlled by one qubit was defaulting to the CNOT gate class. This is fixed now in #1367.

[rnhmjoj]

I can confirm the fix in #1367 works for me.
The workaround I posted before can now be replaced with

def Fgate(diag, n, name='F'):
    '''
    The F(α, β) gate is an multi-controlled gate defined by:

         F(α,β) |1..10> = α|1..10>
         F(α,β) |1..11> = β|1..10>

    Its matrix representation is diag(1,..,1,α,β).
    '''
    gate = q.gates.Unitary(np.diag(diag), n-1, name=name)
    gate = gate.controlled_by(*range(n-1))
    gate.draw_label = name
    return gate