CatTelecorrection/utils.py at main · polyquantique/CatTelecorrection · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
"""
various tools for the simulation of the telecorrection of lossy cats and derived quantities
"""

import numpy as np


def recursive_poisson_list(lmda, up_to):
    """
    a recurrence relation for the Poisson to aid numerical stability at high k and lmda
    slower than SciPy/math factorial, but more stable for high lmda and k

    f(lmda,k) = lmda**k / k!

    will output the sequence of f(lmda,k) up to k=up_to

    """

    p = [1]
    for k in range(1, up_to + 1):
        p += [(lmda / k) * p[k - 1]]

    return np.array(p)


def recursive_poisson_sqrt_list(lmda, up_to):
    """
    a recurrence relation for the Poisson to aid numerical stability at high k and lmda
    slower than SciPy/math factorial, but more stable for high lmda and k

    f(lmda,k) = lmda**k / sqrt(k!)

    will output the sequence of f(lmda,k) up to k=up_to

    """

    p = [1]
    for k in range(1, up_to + 1):
        p += [(lmda / k**0.5) * p[k - 1]]

    return np.array(p)


def Fock_cutoff_calc(alpha, tolerance, limit=400):
    """
    calculates a Fock cutoff based on alpha and tolerance
    sets cutoff at n when cumulative values of

        alpha**n / sqrt(n!)

    are within 10**(-tolerance) of 1
    or if n equals limit

    Args:
        alpha (float): coherent state amplitude
        tolerance (int): cumulative tolerance threshold
        limit (int): maximum cutoff
    Returns:
        (int): Fock cutoff
    """

    value = [np.exp(-abs(alpha) ** 2 / 2)]
    total = value[0] ** 2
    n = 0
    while ((value[-1] / total > 10 ** (-tolerance)) or (n < abs(alpha) ** 2)) and (n < limit):
        n += 1

        value += [value[-1] * abs(alpha) / n**0.5]

        total += value[-1] ** 2

    return n


def Pauli_eigenstates():
    """
    returns the logical Pauli eigenstates in X+, X-, Y+, Y-, Z+, Z- order
    """

    dmX0 = np.array([[1, 1], [1, 1]]) / 2
    dmX1 = np.array([[1, -1], [-1, 1]]) / 2

    dmY0 = np.array([[1, -1j], [1j, 1]]) / 2
    dmY1 = np.array([[1, 1j], [-1j, 1]]) / 2

    dmZ0 = np.array([[1, 0], [0, 0]])
    dmZ1 = np.array([[0, 0], [0, 1]])

    return np.array([dmX0, dmX1, dmY0, dmY1, dmZ0, dmZ1], dtype=np.complex128)