Operator.compose should allow indices (qiskit-community/qiskit-aqua#1144

)

* 1) modified PrimitiveOp._check_zero_for_composition_and_expand, so that it expands shorter operator with identity
2) implemented OperatorBase.expand
3) defined abstract method OperatorBase.identity, and implemented in child classes of OperatorBase (TODO check if implementation makes sense for every class)
4) created tests for composing Operators of different dimensions

* CircuitStateFn.tensor did not set is_measurement parameter (hence always default to False), even if tensored instances were is_measurement=True

* Added test for expand method on StateFn subclasses

* 1) PauliOp.permute implemented
2) test_permute_on_primitive_op implemented. It checks correctness and consistency of permute implementations.

* 1) permute implemented for MatrixOp
2) added testcase for MatrixOp permute, to check consistency with  PauliOp.permute and CircuitOp.permute

* identity renamed to identity_operator in OperatorBase

* OperatorBase.expand renamed to expand_to_dim

* expand_to_dim overridden for PrimitiveOps

* 1) expand_to_dim implemented for each subclass of OperatorBase
2) removed identity_operator

* OperatorStateFn permute implemented

* 1) permute for TensoredOp
2) extended documentation

* 1) TensoredOp.permute moved to ListOp.permute
2) ListOp.permute uses CircuitOps instead of MatrixOps for permutations

* permute defined abstract in OperatorBase

* composition of PrimitiveOp with ComposedOp prepends the ComposedOp.oplist with PrimitiveOp

* Revert "composition of PrimitiveOp with ComposedOp prepends the ComposedOp.oplist with PrimitiveOp"

* 1) compose in MatrixOp, PauliOp and CircuitOp enhanced, to be consistent with ComposedOp.compose
2) test_compose_consistency added to verify the consistency

* DRY applied

* Test if ListOp.permute is consistent with PrimitiveOp permute methods

* test for expand_to_dim on ListOps

* 1) changed signature of compose to allow permutations on operators
2) permutation on operators included in compose method

* refactoring

* unit test for compose with indices

* refactoring and fixed linting

* 1) StateFn.compose expands the shorter operators with identity
2) refactored _check_zero_for_composition_and_expand
3) unit test for new composition features of StateFns
4) fixing style

* 1) expand_to_dim renamed to expand_with_identities (expand_to_dim was misleading, because its parameter is not the target dimension, but the number or qubits to add to operator)
2) documentation for the new features improved

* permute implemented for DictStateFn

* 1) VectorStateFn.to_dict_fn and DictStateFn.to_vector_state_fn implemented
2) VectorStateFn permute implemented
3) unit tests for the new functionality

* implemented custom function to decompose permutations into transpositions, to avoid import of sympy

* explaining comment for CircuitStateFn.expand_with_identity

* expand_with_identity made private and renamed to _expand_dim (expand_with_identity is misleading for implementations in DictStateFn and VectorStateFn, where zeros are used for expansion)

* 1) Compose method has only one set of permutation indices (for other operator).
2) Added new argument front=False in compose method. If front==True, return other.compose(self), but only after expansion and permutation is performed.

* 1) modified test_op_construction because of changes in compose signature
2) fix in CircuitOp.compose

* _check_zero_for_composition_and_expand renamed to _expand_shorter_operator_and_permute, to be more self-explanatory

* VectorStateFn.permute reimplemented for better performance

* removed duplicate method to_vector_state_fn (to_matrix_op was already defined)

* new_self instead of self in EvolvedOp

* test_op_construction refactored

* raise AquaError in ListOp.permute, if ListOp contains operators with different num_qubit

Co-authored-by: Manoel Marques <manoelmrqs@gmail.com>
Co-authored-by: Julien Gacon <gaconju@gmail.com>
Co-authored-by: Manoel Marques <Manoel.Marques@ibm.com>
Co-authored-by: Steve Wood <40241007+woodsp-ibm@users.noreply.github.com>

qiskit/aqua/operators/evolutions/evolved_op.py

@@ -86,13 +86,23 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def _expand_dim(self, num_qubits: int) -> OperatorBase:
+        from qiskit.aqua.operators import I
+        return self.tensor(I ^ num_qubits)
 
+    def permute(self, permutation: List[int]) -> OperatorBase:
+        return EvolvedOp(self.primitive.permute(permutation), coeff=self.coeff)  # type: ignore
+
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        if front:
+            return other.compose(new_self)
         if isinstance(other, ComposedOp):
-            return ComposedOp([self] + other.oplist)  # type: ignore
+            return ComposedOp([new_self] + other.oplist)
 
-        return ComposedOp([self, other])
+        return ComposedOp([new_self, other])
 
     def __str__(self) -> str:
         prim_str = str(self.primitive)

qiskit/aqua/operators/list_ops/composed_op.py

@@ -12,7 +12,7 @@
 
 """ ComposedOp Class """
 
-from typing import List, Union, cast
+from typing import List, Union, cast, Optional
 from functools import reduce, partial
 import numpy as np
 
@@ -82,24 +82,31 @@ def to_circuit(self) -> QuantumCircuit:
     def adjoint(self) -> OperatorBase:
         return ComposedOp([op.adjoint() for op in reversed(self.oplist)], coeff=self.coeff)
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(ComposedOp, new_self)
+
+        if front:
+            return other.compose(new_self)
         # Try composing with last element in list
         if isinstance(other, ComposedOp):
-            return ComposedOp(self.oplist + other.oplist, coeff=self.coeff * other.coeff)
+            return ComposedOp(new_self.oplist + other.oplist, coeff=new_self.coeff * other.coeff)
 
         # Try composing with last element of oplist. We only try
         # this if that last element isn't itself an
         # ComposedOp, so we can tell whether composing the
         # two elements directly worked. If it doesn't,
         # continue to the final return statement below, appending other to the oplist.
-        if not isinstance(self.oplist[-1], ComposedOp):
-            comp_with_last = self.oplist[-1].compose(other)
+        if not isinstance(new_self.oplist[-1], ComposedOp):
+            comp_with_last = new_self.oplist[-1].compose(other)
             # Attempt successful
             if not isinstance(comp_with_last, ComposedOp):
-                new_oplist = self.oplist[0:-1] + [comp_with_last]
-                return ComposedOp(new_oplist, coeff=self.coeff)
+                new_oplist = new_self.oplist[0:-1] + [comp_with_last]
+                return ComposedOp(new_oplist, coeff=new_self.coeff)
 
-        return ComposedOp(self.oplist + [other], coeff=self.coeff)
+        return ComposedOp(new_self.oplist + [other], coeff=new_self.coeff)
 
     def eval(self,
              front: Union[str, dict, np.ndarray,

qiskit/aqua/operators/list_ops/list_op.py

@@ -13,15 +13,18 @@
 """ ListOp Operator Class """
 
 from functools import reduce
-from typing import List, Union, Optional, Callable, Iterator, Set, Dict
+from typing import List, Union, Optional, Callable, Iterator, Set, Dict, cast
 from numbers import Number
 
 import numpy as np
 from scipy.sparse import spmatrix
 
-from qiskit.circuit import ParameterExpression
+from qiskit.circuit import QuantumCircuit, ParameterExpression
+
 from ..legacy.base_operator import LegacyBaseOperator
 from ..operator_base import OperatorBase
+from ... import AquaError
+from ...utils import arithmetic
 
 
 class ListOp(OperatorBase):
@@ -212,11 +215,61 @@ def tensorpower(self, other: int) -> Union[OperatorBase, int]:
         from .tensored_op import TensoredOp
         return TensoredOp([self] * other)
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
+    def _expand_dim(self, num_qubits: int) -> 'ListOp':
+        return ListOp([op._expand_dim(num_qubits + self.num_qubits - op.num_qubits)
+                       for op in self.oplist], combo_fn=self.combo_fn, coeff=self.coeff)
+
+    def permute(self, permutation: List[int]) -> 'ListOp':
+        """Permute the qubits of the operator.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j should be permuted to position permutation[j].
+
+        Returns:
+            A new ListOp representing the permuted operator.
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        new_self = self
+        circuit_size = max(permutation) + 1
+
+        try:
+            if self.num_qubits != len(permutation):
+                raise AquaError("New index must be defined for each qubit of the operator.")
+        except ValueError:
+            raise AquaError("Permute is only possible if all operators in the ListOp have the "
+                            "same number of qubits.") from ValueError
+        if self.num_qubits < circuit_size:
+            # pad the operator with identities
+            new_self = self._expand_dim(circuit_size - self.num_qubits)
+        qc = QuantumCircuit(circuit_size)
+        # extend the indices to match the size of the circuit
+        permutation \
+            = list(filter(lambda x: x not in permutation, range(circuit_size))) + permutation
+
+        # decompose permutation into sequence of transpositions
+        transpositions = arithmetic.transpositions(permutation)
+        for trans in transpositions:
+            qc.swap(trans[0], trans[1])
+
+        from qiskit.aqua.operators import CircuitOp
+
+        return CircuitOp(qc.reverse_ops()) @ new_self @ CircuitOp(qc)
+
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(ListOp, new_self)
+
+        if front:
+            return other.compose(new_self)
         # Avoid circular dependency
         # pylint: disable=cyclic-import,import-outside-toplevel
         from .composed_op import ComposedOp
-        return ComposedOp([self, other])
+        return ComposedOp([new_self, other])
 
     def power(self, exponent: int) -> OperatorBase:
         if not isinstance(exponent, int) or exponent <= 0:
@@ -258,8 +311,8 @@ def to_spmatrix(self) -> Union[spmatrix, List[spmatrix]]:
             [op.to_spmatrix() for op in self.oplist]) * self.coeff  # type: ignore
 
     def eval(self,
-             front: Optional[Union[str, Dict[str, complex], 'OperatorBase']] = None
-             ) -> Union['OperatorBase', float, complex, list]:
+             front: Optional[Union[str, Dict[str, complex], OperatorBase]] = None
+             ) -> Union[OperatorBase, float, complex, list]:
         """
         Evaluate the Operator's underlying function, either on a binary string or another Operator.
         A square binary Operator can be defined as a function taking a binary function to another

qiskit/aqua/operators/list_ops/tensored_op.py

@@ -54,6 +54,16 @@ def num_qubits(self) -> int:
     def distributive(self) -> bool:
         return False
 
+    def _expand_dim(self, num_qubits: int) -> 'TensoredOp':
+        """Appends I ^ num_qubits to ``oplist``. Choice of PauliOp as
+        identity is arbitrary and can be substituted for other PrimitiveOp identity.
+
+        Returns:
+            TensoredOp expanded with identity operator.
+        """
+        from qiskit.aqua.operators import I
+        return TensoredOp(self.oplist + [I ^ num_qubits], coeff=self.coeff)
+
     def tensor(self, other: OperatorBase) -> OperatorBase:
         if isinstance(other, TensoredOp):
             return TensoredOp(self.oplist + other.oplist, coeff=self.coeff * other.coeff)

qiskit/aqua/operators/operator_base.py

@@ -12,7 +12,7 @@
 
 """ OperatorBase Class """
 
-from typing import Set, Union, Dict, Optional, List, cast
+from typing import Set, Union, Dict, Optional, List, cast, Tuple
 from numbers import Number
 from abc import ABC, abstractmethod
 import numpy as np
@@ -452,6 +452,32 @@ def assign_parameters(self,
         """
         raise NotImplementedError
 
+    @abstractmethod
+    def _expand_dim(self, num_qubits: int) -> 'OperatorBase':
+        """Expands the operator with identity operator of dimension 2**num_qubits.
+
+        Returns:
+            Operator corresponding to self.tensor(identity_operator), where dimension of identity
+            operator is 2 ** num_qubits.
+        """
+        raise NotImplementedError
+
+    @abstractmethod
+    def permute(self, permutation: List[int]) -> 'OperatorBase':
+        """Permutes the qubits of the operator.
+
+        Args:
+            permutation: A list defining where each qubit should be permuted. The qubit at index
+                j should be permuted to position permutation[j].
+
+        Returns:
+            A new OperatorBase containing the permuted operator.
+
+        Raises:
+            AquaError: if indices do not define a new index for each qubit.
+        """
+        raise NotImplementedError
+
     def bind_parameters(self,
                         param_dict: Dict[ParameterExpression,
                                          Union[Number,
@@ -501,6 +527,24 @@ def _get_param_dict_for_index(unrolled_dict: Dict[ParameterExpression, List[Numb
         """ Gets a single non-list-nested param_dict for a given list index from a nested one. """
         return {k: v[i] for (k, v) in unrolled_dict.items()}
 
+    def _expand_shorter_operator_and_permute(self, other: 'OperatorBase',
+                                             permutation: Optional[List[int]] = None) \
+            -> Tuple['OperatorBase', 'OperatorBase']:
+        if permutation is not None:
+            other = other.permute(permutation)
+        new_self = self
+        if not self.num_qubits == other.num_qubits:
+            # pylint: disable=cyclic-import,import-outside-toplevel
+            from .operator_globals import Zero
+            if other == Zero:
+                # Zero is special - we'll expand it to the correct qubit number.
+                other = Zero.__class__('0' * self.num_qubits)
+            elif other.num_qubits < self.num_qubits:
+                other = other._expand_dim(self.num_qubits - other.num_qubits)
+            elif other.num_qubits > self.num_qubits:
+                new_self = self._expand_dim(other.num_qubits - self.num_qubits)
+        return new_self, other
+
     # Composition
 
     def __matmul__(self, other: 'OperatorBase') -> 'OperatorBase':
@@ -515,7 +559,8 @@ def __matmul__(self, other: 'OperatorBase') -> 'OperatorBase':
         return self.compose(other)
 
     @abstractmethod
-    def compose(self, other: 'OperatorBase') -> 'OperatorBase':
+    def compose(self, other: 'OperatorBase',
+                permutation: Optional[List[int]] = None, front: bool = False) -> 'OperatorBase':
         r""" Return Operator Composition between self and other (linear algebra-style:
         A@B(x) = A(B(x))), overloaded by ``@``.
 
@@ -529,6 +574,8 @@ def compose(self, other: 'OperatorBase') -> 'OperatorBase':
 
         Args:
             other: The ``OperatorBase`` with which to compose self.
+            permutation: ``List[int]`` which defines permutation on other operator.
+            front: If front==True, return ``other.compose(self)``.
 
         Returns:
             An ``OperatorBase`` equivalent to the function composition of self and other.

qiskit/aqua/operators/primitive_ops/circuit_op.py

@@ -108,30 +108,37 @@ def tensor(self, other: OperatorBase) -> OperatorBase:
 
         return TensoredOp([self, other])
 
-    def compose(self, other: OperatorBase) -> OperatorBase:
-        other = self._check_zero_for_composition_and_expand(other)
+    def compose(self, other: OperatorBase,
+                permutation: Optional[List[int]] = None, front: bool = False) -> OperatorBase:
+
+        new_self, other = self._expand_shorter_operator_and_permute(other, permutation)
+        new_self = cast(CircuitOp, new_self)
+
+        if front:
+            return other.compose(new_self)
+        # ignore
         # pylint: disable=cyclic-import,import-outside-toplevel
         from ..operator_globals import Zero
         from ..state_fns import CircuitStateFn
         from .pauli_op import PauliOp
         from .matrix_op import MatrixOp
 
-        if other == Zero ^ self.num_qubits:
-            return CircuitStateFn(self.primitive, coeff=self.coeff)
+        if other == Zero ^ new_self.num_qubits:
+            return CircuitStateFn(new_self.primitive, coeff=new_self.coeff)
 
         if isinstance(other, (PauliOp, CircuitOp, MatrixOp)):
             other = other.to_circuit_op()
 
         if isinstance(other, (CircuitOp, CircuitStateFn)):
-            new_qc = other.primitive.compose(self.primitive)  # type: ignore
+            new_qc = other.primitive.compose(new_self.primitive)  # type: ignore
             if isinstance(other, CircuitStateFn):
                 return CircuitStateFn(new_qc,
                                       is_measurement=other.is_measurement,
-                                      coeff=self.coeff * other.coeff)
+                                      coeff=new_self.coeff * other.coeff)
             else:
-                return CircuitOp(new_qc, coeff=self.coeff * other.coeff)
+                return CircuitOp(new_qc, coeff=new_self.coeff * other.coeff)
 
-        return super().compose(other)
+        return super(CircuitOp, new_self).compose(other)
 
     def to_matrix(self, massive: bool = False) -> np.ndarray:
         if self.num_qubits > 16 and not massive:
@@ -219,6 +226,9 @@ def reduce(self) -> OperatorBase:
                     del self.primitive.data[i]  # type: ignore
         return self
 
+    def _expand_dim(self, num_qubits: int) -> 'CircuitOp':
+        return self.permute(list(range(num_qubits, num_qubits + self.num_qubits)))
+
     def permute(self, permutation: List[int]) -> 'CircuitOp':
         r"""
         Permute the qubits of the circuit.
@@ -230,5 +240,5 @@ def permute(self, permutation: List[int]) -> 'CircuitOp':
         Returns:
             A new CircuitOp containing the permuted circuit.
         """
-        new_qc = QuantumCircuit(self.num_qubits).compose(self.primitive, qubits=permutation)
+        new_qc = QuantumCircuit(max(permutation) + 1).compose(self.primitive, qubits=permutation)
         return CircuitOp(new_qc, coeff=self.coeff)