Merge branch 'develop' into develop

Author: JakeKitchen

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