Merge pull request #148 from minatoyuichiro/master

Photonqat Integrated

README.md

@@ -242,6 +242,59 @@ else:
     print("timeout")
 ```
 
+# Photonic Continuous Variable Circuit
+
+## Fock basis
+
+### Circuit
+```python
+import photonqat as pq
+import numpy as np
+import matplotlib.pyplot as plt
+
+# mode number = 2, cutoff dimension = 15
+F = pq.Fock(2, cutoff = 15)
+```
+
+### Applying gate
+```python
+alpha = (1 + 1j)
+r = -0.5
+
+F.D(0, alpha) # Displacement to mode 0
+F.S(1, r) # Squeezeng to mode 1
+```
+
+method chain is also available
+```python
+F.D(0, alpha).S(1, r)
+```
+
+### run
+```python
+F.run()
+```
+
+### Plot Wigner function
+```python
+# Plot Wigner fucntion for mode 0 using matplotlib
+(x, p, W) = F.Wigner(0, plot = 'y', xrange = 5.0, prange = 5.0)
+```
+
+## Gaussian formula
+
+### Circuit
+```python
+import photonqat as pq
+import numpy as np
+import matplotlib.pyplot as plt
+
+# mode number = 2
+G = pq.Gaussian(2)
+```
+Applying gate, run the circuit, and plotting Wigner function are also available in same fasion as Fock basis.
+But there are differences in availavle getes and measurement methods.
+
 ### Example1: GHZ
 ```python
 from blueqat import Circuit

blueqat/_version.py

@@ -14,4 +14,5 @@
 
 """The version of blueqat."""
 
-__version__ = "0.5.0-dev"
+__version__ = "0.6.0"
+#__version__ = "0.5.0-dev"

blueqat/photonqat/Fock.py

@@ -0,0 +1,133 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from .Fockbase.gates import *
+from .Fockbase.gateOps import homodyneFock
+from .Fockbase.states import *
+from .Fockbase.WignerFunc import *
+
+STATE_SET = {
+    "vacuum": vacuumState,
+    "coherent": coherentState,
+    "cat": catState,
+    "n_photon": photonNumberState
+}
+
+GATE_SET = {
+    "D": Dgate,
+    "BS": BSgate,
+    "S": Sgate,    
+    "Kerr": Kgate,
+    "MeasF": MeasF,
+    "MeasX": MeasX,
+    "MeasP": MeasP,
+    "polyH": polyH,
+}
+
+class Fock():
+    """
+    Class for continuous variable quantum compting in Fock basis.
+    """
+    def __init__(self, N, cutoff = 10):
+        self.N = N
+        self.cutoff = cutoff
+        self.initState = np.zeros([N, self.cutoff + 1]) + 0j
+        self.initState[:, 0] = 1
+        self.state = None
+        self.initStateFlag = False
+        self.ops = []
+        self.creg = [[None, None, None] for i in range(self.N)] # [x, p, n]
+
+    def  __getattr__(self, name):
+        if name in STATE_SET:
+            if self.initStateFlag:
+                raise ValueError("State must be set before gate operation.")
+            self.ops.append(STATE_SET[name])
+            return self._setGateParam
+
+        elif name in GATE_SET:
+            self.ops.append(GATE_SET[name])
+            self.initStateFlag = True
+            return self._setGateParam
+
+        else:
+            raise AttributeError('The method does not exist')
+
+    def _setGateParam(self, *args, **kwargs):
+        self.ops[-1] = self.ops[-1](self, *args, **kwargs)
+        return self
+
+    def Creg(self, idx, var, scale = 1):
+        """
+        Access to classical register.
+        """
+        return CregReader(self.creg, idx, var, scale)
+
+    def run(self):
+        """
+        Run the circuit.
+        """
+        for op in self.ops:
+            if isinstance(op, STATE):
+                self.initState = op.run(self.initState)
+            elif isinstance(op, GATE):
+                if self.state is None:
+                    self.state = self._multiTensordot()
+                self.state = op.run(self.state)
+                sum_of_prob = np.sum(np.abs(self.state)**2)
+                if np.abs(1 - sum_of_prob) < 0.3:
+                    self.state /= np.sqrt(sum_of_prob)
+                # print(op, np.sum(np.abs(self.state)**2))
+        return self
+
+    def _multiTensordot(self):
+        self.state = self.initState[0, :]
+        for i in range(self.N - 1):
+            self.state = np.tensordot(self.state, self.initState[i+1, :], axes = 0)
+        return self.state
+
+    def Wigner(self, mode, method = 'clenshaw', plot = 'y', xrange = 5.0, prange = 5.0):
+        """
+        Calculate the Wigner function of a selected mode.
+        
+        Args:
+            mode (int): Selecting a optical mode.
+            method: "clenshaw" (default) or "moyal".
+            plot: If 'y', the plot of wigner function is output using matplotlib.\
+                 If 'n', only the meshed values are returned.
+            x(p)range: The range in phase space for calculateing Wigner function.
+        """
+        if self.state is None:
+            self.state = self._multiTensordot()
+            self.initState == None
+        x = np.arange(-xrange, xrange, xrange / 50)
+        p = np.arange(-prange, prange, prange / 50)
+        m = len(x)
+        xx, pp = np.meshgrid(x, -p)
+        W = FockWigner(xx, pp, self.state, mode, method)
+        if plot == 'y':
+            h = plt.contourf(x, p, W)
+            plt.show()
+        return (x, p, W)
+
+    def photonSampling(self, mode, ite = 1):
+        """
+        Simulate the result of photon number resolving measurement for a selected mode.
+        
+        Args:
+            mode (int): Selecting a optical mode.
+            ite: The number of sampling.        
+        """
+        if self.state is None:
+            self.state = self._multiTensordot()
+            self.initState == None
+        reducedDensity = reduceState(self.state, mode)
+        probs = np.real(np.diag(reducedDensity))
+        probs = probs / np.sum(probs)
+        return np.random.choice(probs.shape[0], ite, p = probs)
+
+    def homodyneSampling(self, mode, theta, ite = 1):
+        if self.state is None:
+            self.state = self._multiTensordot()
+            self.initState == None
+        res, psi = homodyneFock(self.state, mode, theta, ite = ite)
+        return res

blueqat/photonqat/Fockbase/WignerFunc.py

@@ -0,0 +1,130 @@
+
+
+"""
+`WignerFunc` module implements calculation of Wigner function.
+This module is internally used.
+"""
+
+import numpy as np
+from scipy.special import factorial as fact
+
+def FockWigner(xmat, pmat, fockState, mode, method = 'clenshaw', tol=1e-10):
+    if fockState.ndim < mode + 1:
+        raise  ValueError("The mode is not exist.")
+    if fockState.ndim > 1:
+        rho = reduceState(fockState, mode)
+    else:
+        rho = np.outer(np.conj(fockState), fockState)
+    if method == 'moyal':
+        W = _Wigner_Moyal(rho, xmat, pmat, tol)
+    elif method == 'clenshaw':
+        W = _Wigner_clenshaw(rho, xmat, pmat, tol)
+    else:
+        raise ValueError("method is invalid.")
+    return W
+
+def reduceState(fockState, mode):
+    modeNum = fockState.ndim
+    cutoff = fockState.shape[-1] - 1
+    fockState = np.swapaxes(fockState, mode, -1)
+    fockState = fockState.flatten()
+    rho = np.outer(np.conj(fockState), fockState)
+    for i in range(modeNum - 1):
+        rho = partialTrace(rho, cutoff)
+    return rho
+
+def partialTrace(rho, cutoff):
+    dim1 = np.int(cutoff + 1)
+    dim2 = np.int(rho.shape[0] / dim1)
+    rho_ = np.zeros([dim2, dim2]) + 0j
+    for j in range(dim1):
+        rho_ += rho[(j * dim2):(j * dim2 + dim2), (j * dim2):(j * dim2 + dim2)]
+    return rho_
+
+def _Wigner_Moyal(rho, xmat, pmat, tol):
+    dim = rho.shape[0]
+    [l, m] = np.indices([dim, dim])
+    A = np.max(np.dstack([l, m]), axis=2)
+    B = np.abs(l - m)
+    C = np.min(np.dstack([l, m]), axis=2)
+    R0 = xmat**2 + pmat**2
+    xmat = xmat[:, :, np.newaxis, np.newaxis]
+    pmat = pmat[:, :, np.newaxis, np.newaxis]
+    R = xmat**2 + pmat**2
+    X = xmat - np.sign(l-m) * 1j * pmat
+    W = 2 * (-1)**C * np.sqrt(2**(B) * fact(C) / fact(A))
+    W = W * np.exp(-R) * X**(B)
+    S = _Sonin(C, B, 2 * R0)
+    W = W * S
+    W = rho * W
+    W = np.sum(np.sum(W, axis = -1), axis = -1)
+    if np.max(np.imag(W)) < tol:
+        W = np.real(W)
+    else:
+        raise ValueError("Wigner plot has imaginary value.")
+    return W
+
+# Based on QuTiP
+def _Wigner_clenshaw(rho, xmat, pmat, tol, hbar = 1):
+    g = np.sqrt(2 / hbar)
+    M = rho.shape[0]
+    A2 = g * (xmat + 1.0j * pmat)    
+    B = np.abs(A2)
+    B *= B
+    w0 = (2*rho[0, -1])*np.ones_like(A2)
+    L = M-1
+    rho = rho * (2*np.ones((M,M)) - np.diag(np.ones(M)))
+    while L > 0:
+        L -= 1
+        w0 = _Wigner_laguerre(L, B, np.diag(rho, L)) + w0 * A2 * (L+1)**-0.5
+    W = w0 * np.exp(-B * 0.5) * (g ** 2 * 0.5 / np.pi)
+    # if np.max(np.imag(W)) < tol:
+    #     W = np.real(W)
+    # else:
+    #     raise ValueError("Wigner plot has imaginary value.")
+    W = np.real(W)
+    return W
+
+def _Wigner_laguerre(L, x, c):
+
+    if len(c) == 1:
+        y0 = c[0]
+        y1 = 0
+    elif len(c) == 2:
+        y0 = c[0]
+        y1 = c[1]
+    else:
+        k = len(c)
+        y0 = c[-2]
+        y1 = c[-1]
+        for i in range(3, len(c) + 1):
+            k -= 1
+            y0,    y1 = c[-i] - y1 * (float((k - 1)*(L + k - 1))/((L+k)*k))**0.5, \
+            y0 - y1 * ((L + 2*k -1) - x) * ((L+k)*k)**-0.5
+
+    return y0 - y1 * ((L + 1) - x) * (L + 1)**-0.5
+
+def _to_2d_ndarray(a):
+    if isinstance(a,(np.ndarray)):
+        return a
+    else:
+        return np.array([[a]])
+
+# slow!
+def _Sonin(n, alpha, x):
+    n = _to_2d_ndarray(n)
+    alpha = _to_2d_ndarray(alpha)
+    x = _to_2d_ndarray(x)
+    a = fact(n + alpha)
+    k0 = np.arange(np.max(n) + 1)
+    k0 = k0[:, np.newaxis, np.newaxis]
+    k = k0 * np.ones([np.max(n) + 1, n.shape[0], n.shape[0]], dtype = np.int)
+    mask = np.ones(k.shape, dtype = np.int)
+    for i in range(k.shape[0]):
+        ind = (np.ones(n.shape) * i) > n
+        mask[i, ind] = 0
+    k *= mask
+    S = mask * (-1)**k * a / fact(n - k) / fact(k + alpha) / fact(k)
+    X = x ** k0
+    S = S[:, np.newaxis, np.newaxis, :, :] * X[:, :, :, np.newaxis, np.newaxis]
+    return np.sum(S, axis = 0)

blueqat/photonqat/Fockbase/__init__.py

@@ -0,0 +1,4 @@
+from .bosonicLadder import *
+from .gates import *
+from .states import *
+from .WignerFunc import *