Version 4.1.0 (#85)

* debug 'to' function; resolve conflict

* debug photonic UAnyGate

* remove old sample_sc_mcmc

* debug omega matrix

* update pyproject

* debug quadrature_to_ladder and ladder_to_quadrature

* update pytest

* update

* version 4.1.0

Author: sansiro77

src/deepquantum/__init__.py

@@ -3,7 +3,7 @@
 DeepQuantum can be directly imported.
 """
 
-__version__ = '4.0.0'
+__version__ = '4.1.0'
 
 
 from . import ansatz

src/deepquantum/circuit.py

@@ -119,8 +119,14 @@ def __add__(self, rhs: 'QubitCircuit') -> 'QubitCircuit':
     def to(self, arg: Any) -> 'QubitCircuit':
         """Set dtype or device of the ``QubitCircuit``."""
         self.init_state.to(arg)
-        self.operators.to(arg)
-        self.observables.to(arg)
+        if arg in (torch.float, torch.double):
+            for op in self.operators:
+                op.to(arg)
+            for ob in self.observables:
+                ob.to(arg)
+        else:
+            self.operators.to(arg)
+            self.observables.to(arg)
         return self
 
     # pylint: disable=arguments-renamed

src/deepquantum/operation.py

@@ -411,6 +411,12 @@ def __init__(
         # MBQC
         self.nodes = copy(self.wires)
 
+    def to(self, arg: Any) -> 'Layer':
+        """Set dtype or device of the ``Layer``."""
+        for gate in self.gates:
+            gate.to(arg)
+        return self
+
     def get_unitary(self) -> torch.Tensor:
         """Get the global unitary matrix."""
         u = None

src/deepquantum/photonic/circuit.py

@@ -132,8 +132,7 @@ def set_init_state(self, init_state: Any) -> None:
             if self.mps:
                 assert self.backend == 'fock' and not self.basis, \
                     'Only support MPS for Fock backend with Fock state tensor.'
-                assert self.cutoff is not None, \
-                    'Please set the cutoff.'
+                assert self.cutoff is not None, 'Please set the cutoff.'
                 self.init_state = MatrixProductState(nsite=self.nmode, state=init_state, chi=self.chi,
                                                      qudit=self.cutoff, normalize=False)
             else:
@@ -191,8 +190,14 @@ def __add__(self, rhs: 'QumodeCircuit') -> 'QumodeCircuit':
     def to(self, arg: Any) -> 'QumodeCircuit':
         """Set dtype or device of the ``QumodeCircuit``."""
         self.init_state.to(arg)
-        self.operators.to(arg)
-        self.measurements.to(arg)
+        if arg in (torch.float, torch.double):
+            for op in self.operators:
+                op.to(arg)
+            for op_m in self.measurements:
+                op_m.to(arg)
+        else:
+            self.operators.to(arg)
+            self.measurements.to(arg)
         if self.backend == 'bosonic' and isinstance(self._bosonic_states, list):
             for bs in self._bosonic_states:
                 bs.to(arg)
@@ -1058,12 +1063,8 @@ def _get_prob_gaussian_base(
             prob = p_vac * torontonian(sub_mat, sub_gamma)
         return abs(prob.real).squeeze()
 
-    def _get_prob_mps(
-        self,
-        final_state: Any,
-        wires: Union[int, List[int], None] = None
-    ) -> torch.Tensor:
-        """Get the probability for the given bit string in mps mode.
+    def _get_prob_mps(self, final_state: Any, wires: Union[int, List[int], None] = None) -> torch.Tensor:
+        """Get the probability of the given bit string for MPS.
 
         Args:
             final_state (Any): The final Fock basis state.
@@ -1077,13 +1078,11 @@ def _get_prob_mps(
         else:
             wires = self._convert_indices(wires)
         assert len(final_state) == len(wires)
-        idx = 0
         state = copy(self.state)
         if self.state[0].ndim == 3:
             state = [site.unsqueeze(0) for site in state]
-        for i in wires:
-            state[i] = state[i][:, :, [final_state[idx]], :]
-            idx += 1
+        for i, wire in enumerate(wires):
+            state[wire] = state[wire][..., [final_state[i]], :]
         return inner_product_mps(state, state).real
 
     def measure(
@@ -1145,6 +1144,7 @@ def _measure_fock(
             assert not mcmc, "Final states have been calculated, we don't need mcmc!"
             return self._measure_dict(shots, with_prob, wires)
         elif isinstance(self.state, List):
+            assert not mcmc, "Final states have been calculated, we don't need mcmc!"
             return self._measure_mps(shots, with_prob, wires)
         else:
             assert False, 'Check your forward function or input!'
@@ -1305,7 +1305,7 @@ def _measure_mps(
         with_prob: bool = False,
         wires: Union[int, List[int], None] = None
     ) -> List[Dict]:
-        """Measure the final state according to mps state."""
+        """Measure the final state according to MPS."""
         if wires is None:
             wires = self.wires
         wires = sorted(self._convert_indices(wires))
@@ -1430,29 +1430,30 @@ def _generate_rand_sample(self, detector: str = 'pnrd'):
             sample = torch.randint(0, self.cutoff, [nmode])
         return sample
 
-    def _generate_chain_sample(
-        self,
-        wires: Union[int, List[int], None] = None
-    ) -> torch.Tensor:
-        """Generate random sample via chain rule.
+    def _generate_chain_sample(self, wires: Union[int, List[int], None] = None) -> torch.Tensor:
+        """Generate random samples via chain rule.
+
+        Args:
+            wires (int, List[int] or None, optional): The wires to measure. It can be an integer or a list of
+                integers specifying the indices of the wires. Default: ``None`` (which means all wires are
+                measured)
 
         Returns:
-                torch.Tensor: sample tensor of shape (batch, len(wires)).
+            torch.Tensor: Tensor of shape (batch, nwire).
         """
         if wires is None:
             wires = self.wires
         wires = sorted(self._convert_indices(wires))
         sample = []
-        mps_state = copy(self.state)
-        if mps_state[0].ndim == 3:
-            mps_state = [site.unsqueeze(0) for site in mps_state]
+        mps = copy(self.state)
+        if mps[0].ndim == 3:
+            mps = [site.unsqueeze(0) for site in mps]
         for i in wires:
-            p = vmap(get_prob_mps)(mps_state, wire=i)
+            p = vmap(get_prob_mps)(mps, wire=i)
             sample_single_wire = torch.multinomial(p, num_samples=1)
             sample.append(sample_single_wire)
-            index = sample_single_wire.reshape(-1).view(-1, 1, 1, 1)\
-                    .expand(-1, mps_state[i].shape[-3], -1, mps_state[i].shape[-1])
-            mps_state[i] = torch.gather(mps_state[i], dim=2, index=index)
+            index = sample_single_wire.reshape(-1, 1, 1, 1).expand(-1, mps[i].shape[-3], -1, mps[i].shape[-1])
+            mps[i] = torch.gather(mps[i], dim=2, index=index)
         sample = torch.stack(sample, dim=-1).squeeze(1)
         return sample
 

src/deepquantum/photonic/gate.py

@@ -164,6 +164,7 @@ def inputs_to_tensor(self, inputs: Any = None) -> torch.Tensor:
 
     def get_matrix(self, theta: Any) -> torch.Tensor:
         """Get the local unitary matrix acting on creation operators."""
+        # correspond to: U a^+ U^+ = u^T @ a^+
         theta = self.inputs_to_tensor(theta)
         if self.inv_mode:
             theta = -theta
@@ -309,6 +310,7 @@ def _add_noise(self, theta: torch.Tensor, phi: torch.Tensor) -> Tuple[torch.Tens
 
     def get_matrix(self, theta: Any, phi: Any) -> torch.Tensor:
         """Get the local unitary matrix acting on creation operators."""
+        # correspond to: U a^+ U^+ = u^T @ a^+
         theta, phi = self.inputs_to_tensor([theta, phi])
         cos = torch.cos(theta)
         sin = torch.sin(theta)
@@ -359,7 +361,6 @@ def get_transform_xp(self, theta: Any, phi: Any) -> Tuple[torch.Tensor, torch.Te
         # correspond to: U a U^+ = (u^*)^T @ a and U^+ a^+ U = u^* @ a^+
         matrix_xp = torch.cat([torch.cat([matrix.real, -matrix.imag], dim=-1),
                                torch.cat([matrix.imag,  matrix.real], dim=-1)], dim=-2).reshape(4, 4)
-        matrix_xp = matrix_xp.to(theta.device, theta.dtype)
         vector_xp = torch.zeros(4, 1, dtype=theta.dtype, device=theta.device)
         return matrix_xp, vector_xp
 
@@ -472,7 +473,8 @@ def __init__(
         self.name = 'MZI'
 
     def get_matrix(self, theta: Any, phi: Any) -> torch.Tensor:
-        """Get the local unitary matrix acting on operators."""
+        """Get the local unitary matrix acting on creation operators."""
+        # correspond to: U a^+ U^+ = u^T @ a^+
         theta, phi = self.inputs_to_tensor([theta, phi])
         cos = torch.cos(theta / 2)
         sin = torch.sin(theta / 2)
@@ -765,6 +767,7 @@ def inputs_to_tensor(self, inputs: Any = None) -> torch.Tensor:
 
     def get_matrix(self, theta: Any) -> torch.Tensor:
         """Get the local unitary matrix acting on creation operators."""
+        # correspond to: U a^+ U^+ = u^T @ a^+
         theta = self.inputs_to_tensor(theta)
         cos = torch.cos(theta / 2) + 0j
         sin = torch.sin(theta / 2) + 0j
@@ -848,8 +851,6 @@ def __init__(
             wires = list(range(minmax[0], minmax[1] + 1))
         super().__init__(name=name, nmode=nmode, wires=wires, cutoff=cutoff, den_mat=den_mat, noise=False)
         self.minmax = [min(self.wires), max(self.wires)]
-        # for i in range(len(self.wires) - 1):
-        #     assert self.wires[i] + 1 == self.wires[i + 1], 'The wires should be consecutive integers'
         if not isinstance(unitary, torch.Tensor):
             unitary = torch.tensor(unitary, dtype=torch.cfloat).reshape(-1, len(self.wires))
         assert unitary.dtype in (torch.cfloat, torch.cdouble)
@@ -879,7 +880,7 @@ def get_matrix_state(self, matrix: torch.Tensor) -> torch.Tensor:
         """
         nt = len(self.wires)
         sqrt = torch.sqrt(torch.arange(self.cutoff, dtype=torch.double, device=matrix.device))
-        tran_mat = matrix.new_zeros([self.cutoff] *  2 * nt)
+        tran_mat = matrix.new_zeros([self.cutoff] * 2 * nt)
         tran_mat[tuple([0] * 2 * nt)] = 1.0
         for rank in range(nt + 1, 2 * nt + 1):
             col_num = rank - nt - 1
@@ -916,10 +917,10 @@ def get_transform_xp(self, matrix: torch.Tensor) -> Tuple[torch.Tensor, torch.Te
         """Get the local affine symplectic transformation acting on quadrature operators in ``xxpp`` order."""
         # correspond to: U a^+ U^+ = u^T @ a^+ and U^+ a U = u @ a
         # correspond to: U a U^+ = (u^*)^T @ a and U^+ a^+ U = u^* @ a^+
+        n = len(self.wires)
         matrix_xp = torch.cat([torch.cat([matrix.real, -matrix.imag], dim=-1),
-                               torch.cat([matrix.imag,  matrix.real], dim=-1)], dim=-2)
-        matrix_xp = matrix_xp.reshape(2 * self.nmode, 2 * self.nmode)
-        vector_xp = torch.zeros(2 * self.nmode, 1, dtype=matrix.real.dtype, device=matrix.real.device)
+                               torch.cat([matrix.imag,  matrix.real], dim=-1)], dim=-2).reshape(2 * n, 2 * n)
+        vector_xp = torch.zeros(2 * n, 1, dtype=matrix.real.dtype, device=matrix.real.device)
         return matrix_xp, vector_xp
 
     def update_transform_xp(self) -> Tuple[torch.Tensor, torch.Tensor]:
@@ -1000,6 +1001,7 @@ def inputs_to_tensor(self, inputs: Any = None) -> Tuple[torch.Tensor, torch.Tens
 
     def get_matrix(self, r: Any, theta: Any) -> torch.Tensor:
         """Get the local matrix acting on annihilation and creation operators."""
+        # correspond to: U^+ (a a^+) U = u @ (a a^+)
         r, theta = self.inputs_to_tensor([r, theta])
         ch = torch.cosh(r)
         sh = torch.sinh(r)
@@ -1155,6 +1157,7 @@ def inputs_to_tensor(self, inputs: Any = None) -> Tuple[torch.Tensor, torch.Tens
 
     def get_matrix(self, r: Any, theta: Any) -> torch.Tensor:
         """Get the local matrix acting on annihilation and creation operators."""
+        # correspond to: U^+ (a a^+) U = u @ (a a^+)
         r, theta = self.inputs_to_tensor([r, theta])
         ch = torch.cosh(r)
         sh = torch.sinh(r)
@@ -1321,6 +1324,7 @@ def _add_noise(self, r: torch.Tensor, theta: torch.Tensor) -> Tuple[torch.Tensor
 
     def get_matrix(self, r: Any, theta: Any) -> torch.Tensor:
         """Get the local unitary matrix acting on annihilation and creation operators."""
+        # correspond to: U^+ (a a^+) U = u @ (a a^+)
         r, theta = self.inputs_to_tensor([r, theta])
         return torch.eye(2, dtype=r.dtype, device=r.device) + 0j
 

src/deepquantum/photonic/operation.py

@@ -372,6 +372,12 @@ def __init__(
         self.ntau = ntau
         self.gates = nn.Sequential()
 
+    def to(self, arg: Any) -> 'Delay':
+        """Set dtype or device of the ``Delay``."""
+        for gate in self.gates:
+            gate.to(arg)
+        return self
+
     def init_para(self, inputs: Any = None) -> None:
         """Initialize the parameters."""
         count = 0

src/deepquantum/photonic/qmath.py

@@ -4,15 +4,11 @@
 
 import itertools
 import warnings
-from collections import defaultdict
-from copy import deepcopy
-from typing import Callable, Dict, Generator, List, Optional, Tuple
+from typing import Dict, Generator, List, Optional, Tuple
 
-import numpy as np
 import torch
 from torch import vmap
 from torch.distributions.multivariate_normal import MultivariateNormal
-from tqdm import tqdm
 
 import deepquantum.photonic as dqp
 from ..qmath import is_unitary
@@ -199,7 +195,7 @@ def xxpp_to_xpxp(matrix: torch.Tensor) -> torch.Tensor:
     """Transform the representation in ``xxpp`` ordering to the representation in ``xpxp`` ordering."""
     nmode = matrix.shape[-2] // 2
     # transformation matrix
-    t = torch.zeros([2 * nmode] * 2, dtype=matrix.dtype, device=matrix.device)
+    t = matrix.new_zeros([2 * nmode] * 2)
     for i in range(2 * nmode):
         if i % 2 == 0:
             t[i][i // 2] = 1
@@ -216,7 +212,7 @@ def xpxp_to_xxpp(matrix: torch.Tensor) -> torch.Tensor:
     """Transform the representation in ``xpxp`` ordering to the representation in ``xxpp`` ordering."""
     nmode = matrix.shape[-2] // 2
     # transformation matrix
-    t = torch.zeros([2 * nmode] * 2, dtype=matrix.dtype, device=matrix.device)
+    t = matrix.new_zeros([2 * nmode] * 2)
     for i in range(2 * nmode):
         if i < nmode:
             t[i][2 * i] = 1
@@ -229,30 +225,48 @@ def xpxp_to_xxpp(matrix: torch.Tensor) -> torch.Tensor:
         return t @ matrix
 
 
-def quadrature_to_ladder(matrix: torch.Tensor) -> torch.Tensor:
-    """Transform the representation in ``xxpp`` ordering to the representation in ``aa^+`` ordering."""
-    nmode = matrix.shape[-2] // 2
-    matrix = matrix + 0j
-    identity = torch.eye(nmode, dtype=matrix.dtype, device=matrix.device)
+def quadrature_to_ladder(tensor: torch.Tensor, symplectic: bool = False) -> torch.Tensor:
+    """Transform the representation in ``xxpp`` ordering to the representation in ``aaa^+a^+`` ordering.
+
+    Args:
+        tensor (torch.Tensor): The input tensor in ``xxpp`` ordering.
+        symplectic (bool, optional): Whether the transformation is applied for symplectic matrix or Gaussian state.
+            Default: ``False`` (which means covariance matrix or displacement vector)
+    """
+    nmode = tensor.shape[-2] // 2
+    tensor = tensor + 0j
+    identity = torch.eye(nmode, dtype=tensor.dtype, device=tensor.device)
     omega = torch.cat([torch.cat([identity, identity * 1j], dim=-1),
-                       torch.cat([identity, identity * -1j], dim=-1)]) * dqp.kappa / dqp.hbar ** 0.5
-    if matrix.shape[-1] == 2 * nmode:
-        return omega @ matrix @ omega.mH
-    elif matrix.shape[-1] == 1:
-        return omega @ matrix
+                       torch.cat([identity, identity * -1j], dim=-1)])
+    if tensor.shape[-1] == 2 * nmode:
+        if symplectic:
+            return omega @ tensor @ omega.mH / 2 # inversed omega
+        else:
+            return omega @ tensor @ omega.mH * dqp.kappa**2 / dqp.hbar
+    elif tensor.shape[-1] == 1:
+        return omega @ tensor * dqp.kappa / dqp.hbar**0.5
 
 
-def ladder_to_quadrature(matrix: torch.Tensor) -> torch.Tensor:
-    """Transform the representation in ``aa^+`` ordering to the representation in ``xxpp`` ordering."""
-    nmode = matrix.shape[-2] // 2
-    matrix = matrix + 0j
-    identity = torch.eye(nmode, dtype=matrix.dtype, device=matrix.device)
+def ladder_to_quadrature(tensor: torch.Tensor, symplectic: bool = False) -> torch.Tensor:
+    """Transform the representation in ``aaa^+a^+`` ordering to the representation in ``xxpp`` ordering.
+
+    Args:
+        tensor (torch.Tensor): The input tensor in ``aaa^+a^+`` ordering.
+        symplectic (bool, optional): Whether the transformation is applied for symplectic matrix or Gaussian state.
+            Default: ``False`` (which means covariance matrix or displacement vector)
+    """
+    nmode = tensor.shape[-2] // 2
+    tensor = tensor + 0j
+    identity = torch.eye(nmode, dtype=tensor.dtype, device=tensor.device)
     omega = torch.cat([torch.cat([identity, identity], dim=-1),
-                       torch.cat([identity * -1j, identity * 1j], dim=-1)]) * dqp.hbar ** 0.5 / (2 * dqp.kappa)
-    if matrix.shape[-1] == 2 * nmode:
-        return (omega @ matrix @ omega.mH).real
-    elif matrix.shape[-1] == 1:
-        return (omega @ matrix).real
+                       torch.cat([identity * -1j, identity * 1j], dim=-1)])
+    if tensor.shape[-1] == 2 * nmode:
+        if symplectic:
+            return (omega @ tensor @ omega.mH).real / 2 # inversed omega
+        else:
+            return (omega @ tensor @ omega.mH).real * dqp.hbar / (4 * dqp.kappa**2)
+    elif tensor.shape[-1] == 1:
+        return (omega @ tensor).real * dqp.hbar**0.5 / (2 * dqp.kappa)
 
 
 def _photon_number_mean_var_gaussian(cov: torch.Tensor, mean: torch.Tensor) -> Tuple[torch.Tensor, torch.Tensor]: