Merge branch 'master' into mcm-post

doc/introduction/templates.rst

@@ -114,6 +114,10 @@ state preparation is typically used as the first operation.
     :description: :doc:`ArbitraryStatePreparation <../code/api/pennylane.ArbitraryStatePreparation>`
     :figure: _static/templates/subroutines/arbitrarystateprep.png
 
+.. gallery-item::
+    :description: :doc:`CosineWindow <../code/api/pennylane.CosineWindow>`
+    :figure: _static/templates/state_preparations/thumbnail_cosine_window.png
+
 .. raw:: html
 
         <div style='clear:both'></div>

doc/releases/changelog-dev.md

@@ -172,8 +172,59 @@
   [(#4620)](https://github.com/PennyLaneAI/pennylane/pull/4620)
   [(#4632)](https://github.com/PennyLaneAI/pennylane/pull/4632)
 
+* The `CosineWindow` template has been added to prepare an initial state based on a cosine wave function.
+  [(#4683)](https://github.com/PennyLaneAI/pennylane/pull/4683)
+
+  ```python
+  import pennylane as qml
+  import matplotlib.pyplot as plt
+
+  dev = qml.device('default.qubit', wires=4)
+
+  @qml.qnode(dev)
+  def example_circuit():
+        qml.CosineWindow(wires=range(4))
+        return qml.state()
+  output = example_circuit()
+
+  # Graph showing state amplitudes
+  plt.bar(range(len(output)), output)
+  plt.show()
+  ```
+
+  We can show how this operator is built:
+
+  ```python
+  import pennylane as qml
+
+  dev = qml.device("default.qubit", wires=5)
+
+  op = qml.CosineWindow(wires=range(5))
+
+  @qml.qnode(dev)
+  def circuit():
+      op.decomposition()
+      return qml.state()
+
+  print(qml.draw(circuit)())
+
+  ```
+
+  ```pycon
+
+  0: ──────────────╭QFT†──Rϕ(1.57)─┤  State
+  1: ──────────────├QFT†──Rϕ(0.79)─┤  State
+  2: ──────────────├QFT†──Rϕ(0.39)─┤  State
+  3: ──────────────├QFT†──Rϕ(0.20)─┤  State
+  4: ──H──RZ(3.14)─╰QFT†──Rϕ(0.10)─┤  State
+
+  ```  
+
 <h3>Improvements 🛠</h3>
 
+* Multi-controlled operations with a single qubit special unitary target can now automatically decompose.
+  [(#4697)](https://github.com/PennyLaneAI/pennylane/pull/4697)
+
 * `pennylane.devices.preprocess` now offers the transforms `decompose`, `validate_observables`, `validate_measurements`,
   `validate_device_wires`, `validate_multiprocessing_workers`, `warn_about_trainable_observables`,
   and `no_sampling` to assist in the construction of devices under the new `devices.Device` API.
@@ -248,6 +299,11 @@
 * Improve builtin types support with `qml.pauli_decompose`.
   [(#4577)](https://github.com/PennyLaneAI/pennylane/pull/4577)
 
+* The function `integrals.py` is modified to replace indexing with slicing in the computationally
+  expensive functions `electron_repulsion` and `_hermite_coulomb`, for a better compatibility with
+  JAX.
+  [(#4685)](https://github.com/PennyLaneAI/pennylane/pull/4685)
+
 * Various changes to measurements to improve feature parity between the legacy `default.qubit` and
   the new `DefaultQubit2`. This includes not trying to squeeze batched `CountsMP` results and implementing
   `MutualInfoMP.map_wires`.
@@ -515,6 +571,7 @@
 
 This release contains contributions from (in alphabetical order):
 
+Guillermo Alonso,
 Utkarsh Azad,
 Jack Brown,
 Stepan Fomichev,

pennylane/ops/op_math/controlled.py

@@ -30,6 +30,7 @@
 from pennylane.wires import Wires
 
 from .symbolicop import SymbolicOp
+from .controlled_decompositions import ctrl_decomp_bisect, ctrl_decomp_zyz
 
 
 def ctrl(op, control, control_values=None, work_wires=None):
@@ -489,8 +490,8 @@ def has_decomposition(self):
             return True
         if isinstance(self.base, qml.PauliX):
             return True
-        # if len(self.base.wires) == 1 and getattr(self.base, "has_matrix", False):
-        #    return True
+        if _is_single_qubit_special_unitary(self.base):
+            return True
         if self.base.has_decomposition:
             return True
 
@@ -569,6 +570,14 @@ def simplify(self) -> "Controlled":
         )
 
 
+def _is_single_qubit_special_unitary(op):
+    if not op.has_matrix or len(op.wires) != 1:
+        return False
+    mat = op.matrix()
+    det = mat[0, 0] * mat[1, 1] - mat[0, 1] * mat[1, 0]
+    return qmlmath.allclose(det, 1)
+
+
 # pylint: disable=protected-access
 def _decompose_no_control_values(op: "operation.Operator") -> List["operation.Operator"]:
     """Provides a decomposition without considering control values. Returns None if
@@ -584,16 +593,10 @@ def _decompose_no_control_values(op: "operation.Operator") -> List["operation.Op
     if isinstance(op.base, qml.PauliX):
         # has some special case handling of its own for further decomposition
         return [qml.MultiControlledX(wires=op.active_wires, work_wires=op.work_wires)]
-    # if (
-    #    len(op.base.wires) == 1
-    #    and len(op.control_wires) >= 2
-    #    and getattr(op.base, "has_matrix", False)
-    #    and qmlmath.get_interface(*op.data) == "numpy"  # as implemented, not differentiable
-    # ):
-    # Bisect algorithms use CNOTs and single qubit unitary
-    #    return ctrl_decomp_bisect(op.base, op.control_wires)
-    # if len(op.base.wires) == 1 and getattr(op.base, "has_matrix", False):
-    #    return ctrl_decomp_zyz(op.base, op.control_wires)
+    if _is_single_qubit_special_unitary(op.base):
+        if len(op.control_wires) >= 2 and qmlmath.get_interface(*op.data) == "numpy":
+            return ctrl_decomp_bisect(op.base, op.control_wires)
+        return ctrl_decomp_zyz(op.base, op.control_wires)
 
     if not op.base.has_decomposition:
         return None

pennylane/qchem/integrals.py

@@ -774,7 +774,7 @@ def _hermite_coulomb(t, u, v, n, p, dr):
     Returns:
         array[float]: value of the Hermite integral
     """
-    x, y, z = dr[0], dr[1], dr[2]
+    x, y, z = dr[0:3]
     T = p * (dr**2).sum(axis=0)
     r = 0
 
@@ -967,12 +967,17 @@ def electron_repulsion(la, lb, lc, ld, ra, rb, rc, rd, alpha, beta, gamma, delta
         + delta * rd[:, np.newaxis, np.newaxis, np.newaxis, np.newaxis]
     ) / (gamma + delta)
 
-    g_t = [expansion(l1, l2, ra[0], rb[0], alpha, beta, t) for t in range(l1 + l2 + 1)]
-    g_u = [expansion(m1, m2, ra[1], rb[1], alpha, beta, u) for u in range(m1 + m2 + 1)]
-    g_v = [expansion(n1, n2, ra[2], rb[2], alpha, beta, v) for v in range(n1 + n2 + 1)]
-    g_r = [expansion(l3, l4, rc[0], rd[0], gamma, delta, r) for r in range(l3 + l4 + 1)]
-    g_s = [expansion(m3, m4, rc[1], rd[1], gamma, delta, s) for s in range(m3 + m4 + 1)]
-    g_w = [expansion(n3, n4, rc[2], rd[2], gamma, delta, w) for w in range(n3 + n4 + 1)]
+    ra0, ra1, ra2 = ra[0:3]
+    rb0, rb1, rb2 = rb[0:3]
+    rc0, rc1, rc2 = rc[0:3]
+    rd0, rd1, rd2 = rd[0:3]
+
+    g_t = [expansion(l1, l2, ra0, rb0, alpha, beta, t) for t in range(l1 + l2 + 1)]
+    g_u = [expansion(m1, m2, ra1, rb1, alpha, beta, u) for u in range(m1 + m2 + 1)]
+    g_v = [expansion(n1, n2, ra2, rb2, alpha, beta, v) for v in range(n1 + n2 + 1)]
+    g_r = [expansion(l3, l4, rc0, rd0, gamma, delta, r) for r in range(l3 + l4 + 1)]
+    g_s = [expansion(m3, m4, rc1, rd1, gamma, delta, s) for s in range(m3 + m4 + 1)]
+    g_w = [expansion(n3, n4, rc2, rd2, gamma, delta, w) for w in range(n3 + n4 + 1)]
 
     g = 0.0
     lengths = [l1 + l2 + 1, m1 + m2 + 1, n1 + n2 + 1, l3 + l4 + 1, m3 + m4 + 1, n3 + n4 + 1]

pennylane/templates/state_preparations/__init__.py

@@ -20,3 +20,4 @@
 from .basis import BasisStatePreparation
 from .arbitrary_state_preparation import ArbitraryStatePreparation
 from .basis_qutrit import QutritBasisStatePreparation
+from .cosine_window import CosineWindow

pennylane/templates/state_preparations/cosine_window.py

@@ -0,0 +1,129 @@
+# Copyright 2018-2021 Xanadu Quantum Technologies Inc.
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#     http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+r"""
+Contains the CosineWindow template.
+"""
+import numpy as np
+import pennylane as qml
+from pennylane.operation import StatePrepBase
+from pennylane import math
+from pennylane.wires import Wires, WireError
+
+
+class CosineWindow(StatePrepBase):
+    r"""CosineWindow(wires)
+    Prepares an initial state with a cosine wave function.
+
+    The wave function is defined below where :math:`m` is the number of wires.
+
+    .. math::
+
+        |\psi\rangle = \sqrt{2^{1-m}} \sum_{k=0}^{2^m-1} \cos(\frac{\pi k}{2^m} - \frac{\pi}{2}) |k\rangle,
+
+    .. figure:: ../../_static/templates/state_preparations/cosine_window.png
+        :align: center
+        :width: 65%
+        :target: javascript:void(0);
+
+    .. note::
+
+        The wave function is shifted by :math:`\frac{\pi}{2}` units so that the window is centered.
+
+    For more details see `Phys. Rev. D 106 (2022) <https://journals.aps.org/prd/abstract/10.1103/PhysRevD.106.034503>`_.
+
+    .. seealso:: :class:`~.QuantumPhaseEstimation` and :class:`~.QFT`.
+
+    Args:
+        wires (Sequence[int] or int): the wire(s) the operation acts on
+
+    **Example**
+
+    >>> dev = qml.device('default.qubit', wires=2)
+    >>> @qml.qnode(dev)
+    ... def example_circuit():
+    ...     qml.CosineWindow(wires=range(2))
+    ...     return qml.probs()
+    >>> print(example_circuit())
+    [1.87469973e-33 2.50000000e-01 5.00000000e-01 2.50000000e-01]
+    """
+
+    @staticmethod
+    def compute_decomposition(wires):  # pylint: disable=arguments-differ,unused-argument
+        r"""Representation of the operator as a product of other operators (static method).
+        It is efficiently decomposed from one QFT over all qubits and one-qubit rotation gates.
+
+        Args:
+            wires (Iterable, Wires): the wire(s) the operation acts on
+
+        Returns:
+            list[Operator]: decomposition into lower level operations
+        """
+
+        decomp_ops = []
+
+        decomp_ops.append(qml.Hadamard(wires=wires[-1]))
+        decomp_ops.append(qml.RZ(np.pi, wires=wires[-1]))
+        decomp_ops.append(qml.adjoint(qml.QFT)(wires=wires))
+
+        for ind, wire in enumerate(wires):
+            decomp_ops.append(qml.PhaseShift(np.pi * 2 ** (-ind - 1), wires=wire))
+
+        return decomp_ops
+
+    def label(self, decimals=None, base_label=None, cache=None):
+        return "CosineWindow"
+
+    def state_vector(self, wire_order=None):  # pylint: disable=arguments-differ,unused-argument
+        r"""Calculation of the state vector generated by the cosine window.
+
+        Args:
+            wire_order (Iterable, Wires): Custom order of wires for the returned state vector.
+
+        Raises:
+            WireError: Custom wire_order must contain all wires.
+
+        Returns:
+            TensorLike[complex]: output state
+        """
+
+        num_op_wires = len(self.wires)
+        op_vector_shape = (2,) * num_op_wires
+        vector = np.array(
+            [
+                np.sqrt(2 ** (1 - num_op_wires))
+                * np.cos(-np.pi / 2 + np.pi * x / 2**num_op_wires)
+                for x in range(2**num_op_wires)
+            ]
+        )
+        op_vector = math.reshape(vector, op_vector_shape)
+
+        if wire_order is None or Wires(wire_order) == self.wires:
+            return op_vector
+
+        wire_order = Wires(wire_order)
+        if not wire_order.contains_wires(self.wires):
+            raise WireError(f"Custom wire_order must contain all {self.name} wires")
+
+        indices = tuple([Ellipsis] + [slice(None)] * num_op_wires)
+
+        ket_shape = [2] * num_op_wires
+        ket = np.zeros(ket_shape, dtype=np.complex128)
+        ket[indices] = op_vector
+
+        if self.wires != wire_order[:num_op_wires]:
+            current_order = self.wires + list(Wires.unique_wires([wire_order, self.wires]))
+            desired_order = [current_order.index(w) for w in wire_order]
+            ket = ket.transpose(desired_order)
+
+        return math.convert_like(ket, op_vector)

tests/ops/op_math/test_controlled.py

@@ -263,6 +263,12 @@ def test_has_decomposition_true_via_pauli_x(self):
         op = Controlled(qml.PauliX(3), [0, 4])
         assert op.has_decomposition is True
 
+    def test_has_decomposition_multicontrolled_special_unitary(self):
+        """Test that a one qubit special unitary with any number of control
+        wires has a decomposition."""
+        op = Controlled(qml.RX(1.234, wires=0), (1, 2, 3, 4, 5))
+        assert op.has_decomposition
+
     def test_has_decomposition_true_via_base_has_decomp(self):
         """Test that Controlled claims `has_decomposition` to be true if
         the base has a decomposition and indicates this via `has_decomposition`."""
@@ -813,6 +819,38 @@ def test_None_default(self):
         op = Controlled(TempOperator(0), (1, 2))
         assert _decompose_no_control_values(op) is None
 
+    def test_non_differentiable_one_qubit_special_unitary(self):
+        """Assert that a non differentiable on qubit special unitary uses the bisect decomposition."""
+        op = qml.ctrl(qml.RZ(1.2, wires=0), (1, 2, 3, 4))
+        decomp = op.decomposition()
+
+        assert qml.equal(decomp[0], qml.MultiControlledX(wires=(1, 2, 0), work_wires=(3, 4)))
+        assert isinstance(decomp[1], qml.QubitUnitary)
+        assert qml.equal(decomp[2], qml.MultiControlledX(wires=(3, 4, 0), work_wires=(1, 2)))
+        assert isinstance(decomp[3].base, qml.QubitUnitary)
+        assert qml.equal(decomp[4], qml.MultiControlledX(wires=(1, 2, 0), work_wires=(3, 4)))
+        assert isinstance(decomp[5], qml.QubitUnitary)
+        assert qml.equal(decomp[6], qml.MultiControlledX(wires=(3, 4, 0), work_wires=(1, 2)))
+        assert isinstance(decomp[7].base, qml.QubitUnitary)
+
+        decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
+        assert qml.math.allclose(op.matrix(), decomp_mat)
+
+    def test_differentiable_one_qubit_special_unitary(self):
+        """Assert that a differentiable qubit speical unitary uses the zyz decomposition."""
+
+        op = qml.ctrl(qml.RZ(qml.numpy.array(1.2), 0), (1, 2, 3, 4))
+        decomp = op.decomposition()
+
+        assert qml.equal(decomp[0], qml.RZ(qml.numpy.array(1.2), 0))
+        assert qml.equal(decomp[1], qml.MultiControlledX(wires=(1, 2, 3, 4, 0)))
+        assert qml.equal(decomp[2], qml.RZ(qml.numpy.array(-0.6), wires=0))
+        assert qml.equal(decomp[3], qml.MultiControlledX(wires=(1, 2, 3, 4, 0)))
+        assert qml.equal(decomp[4], qml.RZ(qml.numpy.array(-0.6), wires=0))
+
+        decomp_mat = qml.matrix(op.decomposition, wire_order=op.wires)()
+        assert qml.math.allclose(op.matrix(), decomp_mat)
+
 
 @pytest.mark.parametrize("test_expand", (False, True))
 class TestDecomposition: