An initial implementation of the phase noise channel

mrmustard/lab/gates.py

@@ -933,7 +933,7 @@ class Amplifier(Channel):
     """
 
     is_gaussian = True
-    short_name = "Amp"
+    short_name = "Amp~"
     parallelizable = True
 
     def __init__(

mrmustard/lab_dev/states/base.py

@@ -52,6 +52,7 @@
 from ..circuit_components import CircuitComponent
 from ..circuit_components_utils import BtoPS
 
+
 __all__ = ["State"]
 
 # ~~~~~~~
@@ -333,6 +334,7 @@ def phase_space(self, s: float) -> tuple:
             Returns:
                 The covariance matrix, the mean vector and the coefficient of the state in s-parametrized phase space.
         """
+
         if not isinstance(self.ansatz, PolyExpAnsatz):
             raise ValueError("Can calculate phase space only for Bargmann states.")
 

mrmustard/lab_dev/states/dm.py

@@ -36,6 +36,7 @@
 from .base import State, _validate_operator, OperatorType
 from ..circuit_components import CircuitComponent
 from ..circuit_components_utils import BtoQ, TraceOut
+
 from ..utils import shape_check
 
 __all__ = ["DM"]
@@ -196,6 +197,7 @@ def from_quadrature(
             ValueError: If the given triple has shapes that are inconsistent
                 with the number of modes.
         """
+
         QtoB = BtoQ(modes, phi).inverse()
         Q = DM.from_ansatz(modes, PolyExpAnsatz(*triple))
         return DM.from_ansatz(modes, (Q >> QtoB).ansatz, name)
@@ -302,6 +304,7 @@ def expectation(self, operator: CircuitComponent):
             ValueError: If ``operator`` is defined over a set of modes that is not a subset of the
                 modes of this state.
         """
+
         op_type, msg = _validate_operator(operator)
         if op_type is OperatorType.INVALID_TYPE:
             raise ValueError(msg)
@@ -413,6 +416,7 @@ def __rshift__(self, other: CircuitComponent) -> CircuitComponent:
         Returns a ``DM`` when the wires of the resulting components are compatible with
         those of a ``DM``, a ``CircuitComponent`` otherwise, and a scalar if there are no wires left.
         """
+
         result = super().__rshift__(other)
         if not isinstance(result, CircuitComponent):
             return result  # scalar case handled here

mrmustard/lab_dev/states/ket.py

@@ -152,6 +152,7 @@ def from_quadrature(
         phi: float = 0.0,
         name: str | None = None,
     ) -> State:
+
         QtoB = BtoQ(modes, phi).inverse()
         Q = Ket.from_ansatz(modes, PolyExpAnsatz(*triple))
         return Ket.from_ansatz(modes, (Q >> QtoB).ansatz, name)
@@ -263,6 +264,7 @@ def expectation(self, operator: CircuitComponent):
             ValueError: If ``operator`` is defined over a set of modes that is not a subset of the
                 modes of this state.
         """
+
         op_type, msg = _validate_operator(operator)
         if op_type is OperatorType.INVALID_TYPE:
             raise ValueError(msg)
@@ -371,6 +373,7 @@ def __rshift__(self, other: CircuitComponent | Scalar) -> CircuitComponent | Bat
         with those of a ``DM`` or of a ``Ket``. Returns a ``CircuitComponent`` in general,
         and a (batched) scalar if there are no wires left, for convenience.
         """
+
         result = super().__rshift__(other)
         if not isinstance(result, CircuitComponent):
             return result  # scalar case handled here

mrmustard/lab_dev/transformations/__init__.py

@@ -26,6 +26,7 @@
 from .ggate import *
 from .gaussrandnoise import *
 from .identity import *
+from .phasenoise import *
 from .rgate import *
 from .s2gate import *
 from .sgate import *

mrmustard/lab_dev/transformations/amplifier.py

@@ -84,7 +84,7 @@ def __init__(
         gain_trainable: bool = False,
         gain_bounds: tuple[float | None, float | None] = (1.0, None),
     ):
-        super().__init__(name="Amp")
+        super().__init__(name="Amp~")
         (gs,) = list(reshape_params(len(modes), gain=gain))
         self._add_parameter(
             make_parameter(

mrmustard/lab_dev/transformations/attenuator.py

@@ -75,7 +75,7 @@ class Attenuator(Channel):
             c &= 1\:.
     """
 
-    short_name = "Att"
+    short_name = "Att~"
 
     def __init__(
         self,
@@ -84,7 +84,7 @@ def __init__(
         transmissivity_trainable: bool = False,
         transmissivity_bounds: tuple[float | None, float | None] = (0.0, 1.0),
     ):
-        super().__init__(name="Att")
+        super().__init__(name="Att~")
         (etas,) = list(reshape_params(len(modes), transmissivity=transmissivity))
         self._add_parameter(
             make_parameter(

mrmustard/lab_dev/transformations/gaussrandnoise.py

@@ -74,7 +74,7 @@ def __init__(
         if (math.real(math.eigvals(Y)) >= -settings.ATOL).min() == 0:
             raise ValueError("The input Y matrix has negative eigen-values.")
 
-        super().__init__(name="GRN")
+        super().__init__(name="GRN~")
         self._add_parameter(make_parameter(Y_trainable, value=Y, name="Y", bounds=(None, None)))
         self._representation = self.from_ansatz(
             modes_in=modes,

mrmustard/lab_dev/transformations/phasenoise.py

@@ -0,0 +1,86 @@
+# Copyright 2024 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.
+
+"""
+The class representing a Phase noise channel.
+"""
+
+from __future__ import annotations
+from typing import Sequence
+from mrmustard.lab_dev.circuit_components import CircuitComponent
+from mrmustard.physics.ansatz.array_ansatz import ArrayAnsatz
+from mrmustard.physics.representations import Representation
+from mrmustard import math
+import numpy as np
+from .base import Channel
+from ..utils import make_parameter
+
+__all__ = ["PhaseNoise"]
+
+
+class PhaseNoise(Channel):
+    r"""
+    The Phase noise channel.
+
+    This class represents the application of a random phase. The distributiuon of the phase
+    is assumed to be a Gaussian with mean zero, and standard deviation `phase_stdev`.
+
+    Args:
+        modes: The modes the channel is applied to
+        phase_stdev: The standard deviation of the random phase noise.
+
+    ..details::
+        The Fock representation is connected to the Fourier coefficients of the distribution.
+    """
+
+    short_name = "P~"
+
+    def __init__(
+        self,
+        modes: Sequence[int],
+        phase_stdev: float | Sequence[float],
+        phase_stdev_trainable: bool = False,
+        phase_stdev_bounds: tuple[float | None, float | None] = (0.0, None),
+    ):
+        super().__init__(name="PhaseNoise")
+        self._add_parameter(
+            make_parameter(phase_stdev_trainable, phase_stdev, "phase_stdev", phase_stdev_bounds)
+        )
+        self._representation = self.from_ansatz(
+            modes_in=modes, modes_out=modes, ansatz=None
+        ).representation
+
+    def __custom_rrshift__(self, other: CircuitComponent) -> CircuitComponent:
+        r"""
+        Since PhaseNoise admits a particularly nice form in the Fock basis, we have implemented its right-shift operation separately.
+
+        Args:
+            other: the component other than the PhaseNoise object that is present in the contraction
+
+        Output:
+            the result of the contraction.
+        """
+
+        if not other.wires.bra or not other.wires.ket:
+            other = other @ other.adjoint
+        array = math.asnumpy(other.fock_array())
+        mode_indices = np.indices(array.shape)
+        for mode in self.modes:
+            phase_factors = np.exp(
+                -0.5
+                * (mode_indices[mode] - mode_indices[other.n_modes + mode]) ** 2
+                * self.phase_stdev.value**2
+            )
+            array *= phase_factors
+        return CircuitComponent(Representation(ArrayAnsatz(array, False), other.wires), self.name)

mrmustard/physics/triples.py

@@ -835,7 +835,7 @@ def XY_to_channel_Abc(X: RealMatrix, Y: RealMatrix, d: Vector | None = None) ->
 
     A = math.Xmat(2 * m) @ R @ xi_inv_in_blocks @ math.conj(R).T
     temp = math.block([[(xi_inv @ d).reshape(2 * m, 1)], [(-X.T @ xi_inv @ d).reshape((2 * m, 1))]])
-    b = 1 / math.sqrt(settings.HBAR) * math.conj(R) @ temp
+    b = 1 / math.sqrt(complex(settings.HBAR)) * math.conj(R) @ temp
     c = math.exp(-0.5 / settings.HBAR * d @ xi_inv @ d) / math.sqrt(math.det(xi))
 
     return A, b, c

tests/test_lab_dev/test_transformations/test_amplifier.py

@@ -36,7 +36,7 @@ class TestAmplifier:
     def test_init(self, modes, gain):
         gate = Amplifier(modes, gain)
 
-        assert gate.name == "Amp"
+        assert gate.name == "Amp~"
         assert gate.modes == [modes] if not isinstance(modes, list) else sorted(modes)
 
     def test_init_error(self):

tests/test_lab_dev/test_transformations/test_attenuator.py

@@ -35,7 +35,7 @@ class TestAttenuator:
     def test_init(self, modes, transmissivity):
         gate = Attenuator(modes, transmissivity)
 
-        assert gate.name == "Att"
+        assert gate.name == "Att~"
         assert gate.modes == [modes] if not isinstance(modes, list) else sorted(modes)
 
     def test_init_error(self):

tests/test_lab_dev/test_transformations/test_gaussrandnoise.py

@@ -33,7 +33,7 @@ def test_init(self):
 
         a = np.random.random((2, 2))
         grn = GaussRandNoise([0], a @ a.T)
-        assert grn.name == "GRN"
+        assert grn.name == "GRN~"
         assert grn.modes == [0]
 
     def test_grn(self):

tests/test_lab_dev/test_transformations/test_phasenoise.py

@@ -0,0 +1,75 @@
+# Copyright 2024 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.
+
+"""Tests for the ``PhaseNoise`` class."""
+
+# pylint: disable=missing-function-docstring, expression-not-assigned
+
+import numpy as np
+import pytest
+
+from mrmustard import math
+from mrmustard.lab_dev.circuit_components import CircuitComponent
+from mrmustard.lab_dev.states import DM, Coherent, Ket, Number
+from mrmustard.lab_dev.transformations import Dgate, FockDamping, PhaseNoise
+
+
+class TestPhaseNoise:
+    r"""
+    Tests for the ``PhaseNoise`` class.
+    """
+
+    def test_init(self):
+        "Tests the PhaseNoise initialization."
+        ch = PhaseNoise([0, 1], 0.2)
+        assert ch.name == "PhaseNoise"
+        assert ch.phase_stdev.value == 0.2
+        assert ch.modes == [0, 1]
+        assert ch.ansatz is None
+
+    def test_application(self):
+        "Tests application of PhaseNoise on Ket and DM"
+        psi_1 = Ket.random([0, 1]) >> Dgate([0], 0.5, 0.5) >> PhaseNoise([0], 0.2)
+        assert isinstance(psi_1, DM)
+        assert psi_1.purity < 1
+
+        psi_2 = Coherent([0], 2)
+        phi = psi_2 >> PhaseNoise([0], 10)
+        after_noise_array = phi.fock_array(10)
+        assert math.allclose(
+            np.diag(after_noise_array), np.diag(psi_2.dm().fock_array(10))
+        )  # the diagonal entries must remain unchanged
+        mask = ~math.eye(after_noise_array.shape[0], dtype=bool)
+        assert math.allclose(
+            after_noise_array[mask], math.zeros_like(after_noise_array[mask])
+        )  # the off-diagonal entries must vanish
+
+        rho = DM.random([0, 1]) >> Dgate([0], 0.5, 0.5) >> PhaseNoise([0], 0.2)
+        assert isinstance(rho, DM)
+        assert rho.purity < 1
+
+    @pytest.mark.parametrize("sigma", [0.2, 0.5, 0.7])
+    def test_numeric(self, sigma):
+        r"""
+        A numeric example
+        """
+        psi = Number([0], 0) + Number([0], 1)
+        phi = psi >> PhaseNoise([0], sigma)
+        assert phi.fock_array(2)[0, 1] == math.exp(-(complex(sigma) ** 2) / 2)
+
+    def test_check_adding_adjoint(self):
+        r"""
+        Tests if the PhaseNoise custum rrshift correcly adds the adjoint.
+        """
+        assert isinstance(FockDamping([0], 0.5) >> PhaseNoise([0], 0.2), CircuitComponent)

tests/test_lab_dev/test_transformations/test_transformations_base.py

@@ -171,8 +171,8 @@ def test_repr(self):
         channel1 = Attenuator([0, 1], 0.9)
         ch_component = CircuitComponent(channel1.representation, channel1.name)
 
-        assert repr(channel1) == "Attenuator(modes=[0, 1], name=Att, repr=PolyExpAnsatz)"
-        assert repr(ch_component) == "CircuitComponent(modes=[0, 1], name=Att, repr=PolyExpAnsatz)"
+        assert repr(channel1) == "Attenuator(modes=[0, 1], name=Att~, repr=PolyExpAnsatz)"
+        assert repr(ch_component) == "CircuitComponent(modes=[0, 1], name=Att~, repr=PolyExpAnsatz)"
 
     def test_inverse_channel(self):
         gate = Sgate([0], 0.1, 0.2) >> Dgate([0], 0.1, 0.2) >> Attenuator([0], 0.5)