Skip to content

Commit e563c23

Browse files
yutuer21yutuer21
committed
add thermalrelaxtionchannel and test
1 parent a0c82ec commit e563c23

File tree

2 files changed

+289
-16
lines changed

2 files changed

+289
-16
lines changed

‎tensorcircuit/channels.py‎

Lines changed: 216 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -3,16 +3,18 @@
33
"""
44

55
import sys
6-
from typing import Any, Sequence, Union
6+
from typing import Any, Sequence, Union, Optional, Dict
77
from operator import and_
88
from functools import reduce
99
from functools import partial
1010
import numpy as np
1111

12+
1213
from . import cons
1314
from . import interfaces
1415
from .cons import backend, dtypestr
1516
from . import gates
17+
from .gates import array_to_tensor
1618

1719

1820
thismodule = sys.modules[__name__]
@@ -306,6 +308,128 @@ def phasedampingchannel(gamma: float) -> Sequence[Gate]:
306308
return [m0, m1]
307309

308310

311+
def thermalrelaxationchannel(
312+
t1: float,
313+
t2: float,
314+
time: float,
315+
method: str = "general",
316+
excitedstatepopulation: float = 0.0,
317+
) -> Sequence[Gate]:
318+
r"""
319+
Return a thermal_relaxation_channel
320+
321+
:param t1: the T1 relaxation time.
322+
:type t1: float
323+
:param t2: the T2 dephasing time.
324+
:type t2: float
325+
:param time: gate time
326+
:type time: float
327+
:param method: method to express error (default: "general")
328+
:type time: str
329+
:param excitedstatepopulation: the population of state :math:`|1\rangle` at equilibrium (default: 0)
330+
:type excited_state_population: float, optional
331+
:return: A thermal_relaxation_channel
332+
:rtype: Sequence[Gate]
333+
"""
334+
t1 = backend.cast(array_to_tensor(t1), dtype=cons.dtypestr)
335+
t2 = backend.cast(array_to_tensor(t2), dtype=cons.dtypestr)
336+
time = backend.cast(array_to_tensor(time), dtype=cons.dtypestr)
337+
338+
# T1 relaxation rate: :math:`|1\rangle \rightarrow \0\rangle`
339+
rate1 = 1.0 / t1
340+
p_reset = 1 - backend.exp(-time * rate1)
341+
342+
# T2 dephasing rate: :math:`|+\rangle \rightarrow \-\rangle`
343+
rate2 = 1.0 / t2
344+
exp_t2 = backend.exp(-time * rate2)
345+
346+
# Qubit state equilibrium probabilities
347+
p0 = 1 - excitedstatepopulation
348+
p1 = excitedstatepopulation
349+
p0 = backend.cast(array_to_tensor(p0), dtype=cons.dtypestr)
350+
p1 = backend.cast(array_to_tensor(p1), dtype=cons.dtypestr)
351+
352+
if method == "T1dom":
353+
# git avaliable
354+
m0 = backend.convert_to_tensor(
355+
np.array([[1, 0], [0, 0]], dtype=cons.npdtype)
356+
) # reset channel
357+
m1 = backend.convert_to_tensor(
358+
np.array([[0, 1], [0, 0]], dtype=cons.npdtype)
359+
) # reset channel
360+
m2 = backend.convert_to_tensor(
361+
np.array([[0, 0], [1, 0]], dtype=cons.npdtype)
362+
) # X gate + rest channel
363+
m3 = backend.convert_to_tensor(
364+
np.array([[0, 0], [0, 1]], dtype=cons.npdtype)
365+
) # X gate + rest channel
366+
367+
tup = [
368+
gates.i().tensor, # type: ignore
369+
gates.z().tensor, # type: ignore
370+
m0,
371+
m1,
372+
m2,
373+
m3,
374+
]
375+
376+
p_reset0 = p_reset * p0
377+
p_reset1 = p_reset * p1
378+
p_z = (
379+
(1 - p_reset) * (1 - backend.exp(-time * (rate2 - rate1))) / 2
380+
) # must have rate2 > rate1
381+
p_identity = 1 - p_z - p_reset0 - p_reset1
382+
probs = [p_identity, p_z, p_reset0, p_reset0, p_reset1, p_reset1]
383+
384+
Gkraus = []
385+
for pro, paugate in zip(probs, tup):
386+
Gkraus.append(Gate(_sqrt(pro) * paugate))
387+
return Gkraus
388+
389+
else: # method == "general" or method == "T2dom":
390+
# git avaliable
391+
choi = (1 - p1 * p_reset) * backend.convert_to_tensor(
392+
np.array(
393+
[[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]],
394+
dtype=cons.npdtype,
395+
)
396+
)
397+
choi += (exp_t2) * backend.convert_to_tensor(
398+
np.array(
399+
[[0, 0, 0, 1], [0, 0, 0, 0], [0, 0, 0, 0], [1, 0, 0, 0]],
400+
dtype=cons.npdtype,
401+
)
402+
)
403+
choi += (p1 * p_reset) * backend.convert_to_tensor(
404+
np.array(
405+
[[0, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]],
406+
dtype=cons.npdtype,
407+
)
408+
)
409+
choi += (p0 * p_reset) * backend.convert_to_tensor(
410+
np.array(
411+
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1, 0], [0, 0, 0, 0]],
412+
dtype=cons.npdtype,
413+
)
414+
)
415+
choi += (1 - p0 * p_reset) * backend.convert_to_tensor(
416+
np.array(
417+
[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]],
418+
dtype=cons.npdtype,
419+
)
420+
)
421+
422+
nmax = 4
423+
if (
424+
np.abs(excitedstatepopulation - 0.0) < 1e-3
425+
or np.abs(excitedstatepopulation - 1.0) < 1e-3
426+
):
427+
nmax = 3
428+
429+
listKraus = choi_to_kraus(choi, truncation_rules={"max_singular_values": nmax})
430+
return [Gate(i) for i in listKraus]
431+
432+
309433
def _collect_channels() -> Sequence[str]:
310434
r"""
311435
Return channels names in this module.
@@ -496,7 +620,63 @@ def reshuffle(op: Matrix, order: Sequence[int]) -> Matrix:
496620
argnums=[0],
497621
gate_to_tensor=True,
498622
)
499-
def choi_to_kraus(choi: Matrix, atol: float = 1e-5) -> Matrix:
623+
# def choi_to_kraus(choi: Matrix, atol: float = 1e-5, Maxkraus: Optional[int] = None) -> Matrix:
624+
# r"""
625+
# Convert the Choi matrix representation to Kraus operator representation.
626+
627+
# This can be done by firstly geting eigen-decomposition of Choi-matrix
628+
629+
# .. math::
630+
# \Lambda = \sum_k \gamma_k \vert \phi_k \rangle \langle \phi_k \vert
631+
632+
# Then define Kraus operators
633+
634+
# .. math::
635+
# K_k = \sqrt{\gamma_k} V_k
636+
637+
# where :math:`\gamma_k\geq0` and :math:`\phi_k` is the col-val vectorization of :math:`V_k` .
638+
639+
640+
# :Examples:
641+
642+
643+
# >>> kraus = tc.channels.phasedampingchannel(0.2)
644+
# >>> superop = tc.channels.kraus_to_choi(kraus)
645+
# >>> kraus_new = tc.channels.choi_to_kraus(superop)
646+
647+
648+
# :param choi: Choi matrix
649+
# :type choi: Matrix
650+
# :return: A list of Kraus operators
651+
# :rtype: Sequence[Matrix]
652+
# """
653+
# dim = backend.shape_tuple(choi)
654+
# input_dim = _safe_sqrt(dim[0])
655+
# output_dim = _safe_sqrt(dim[1])
656+
657+
# #if is_hermitian_matrix(choi, atol=atol):
658+
# # Get eigen-decomposition of Choi-matrix
659+
# e, v = backend.eigh(choi) # value of e is from minimal to maxmal
660+
661+
# # CP-map Kraus representation
662+
# kraus = []
663+
# for val, vec in zip(backend.real(e), backend.transpose(v)):
664+
# if val > atol:
665+
# k = backend.sqrt(backend.cast(val, dtypestr)) * backend.transpose(
666+
# backend.reshape(vec, [output_dim, input_dim]), [1, 0]
667+
# )
668+
# kraus.append(k)
669+
670+
# if not kraus:
671+
# kraus.append(backend.zeros([output_dim, input_dim], dtype=dtypestr))
672+
# return kraus
673+
674+
# # raise ValueError("illegal Choi matrix")
675+
676+
677+
def choi_to_kraus(
678+
choi: Matrix, truncation_rules: Optional[Dict[str, Any]] = None
679+
) -> Matrix:
500680
r"""
501681
Convert the Choi matrix representation to Kraus operator representation.
502682
@@ -530,24 +710,44 @@ def choi_to_kraus(choi: Matrix, atol: float = 1e-5) -> Matrix:
530710
input_dim = _safe_sqrt(dim[0])
531711
output_dim = _safe_sqrt(dim[1])
532712

533-
if is_hermitian_matrix(choi, atol=atol):
534-
# Get eigen-decomposition of Choi-matrix
535-
e, v = backend.eigh(choi)
713+
# Get eigen-decomposition of Choi-matrix
714+
e, v = backend.eigh(choi) # value of e is from minimal to maxmal
715+
e = backend.real(e)
716+
v = backend.transpose(v)
717+
718+
# CP-map Kraus representation
719+
kraus = []
720+
721+
if truncation_rules is None:
722+
truncation_rules = {}
723+
724+
if truncation_rules.get("max_singular_values", None) is not None:
536725

537-
# CP-map Kraus representation
538-
kraus = []
539-
for val, vec in zip(backend.real(e), backend.transpose(v)):
540-
if val > atol:
541-
k = backend.sqrt(backend.cast(val, dtypestr)) * backend.transpose(
542-
backend.reshape(vec, [output_dim, input_dim]), [1, 0]
726+
nkraus = truncation_rules["max_singular_values"]
727+
for i in range(nkraus):
728+
k = backend.sqrt(backend.cast(e[-(i + 1)], dtypestr)) * backend.transpose(
729+
backend.reshape(v[-(i + 1)], [output_dim, input_dim]), [1, 0]
730+
)
731+
kraus.append(k)
732+
733+
else:
734+
if truncation_rules.get("max_singular_values", None) is None:
735+
atol = 1e-5
736+
else:
737+
atol = truncation_rules["max_truncattion_err"]
738+
739+
for i in range(dim[0]):
740+
if e[-(i + 1)] > atol:
741+
k = backend.sqrt(
742+
backend.cast(e[-(i + 1)], dtypestr)
743+
) * backend.transpose(
744+
backend.reshape(v[-(i + 1)], [output_dim, input_dim]), [1, 0]
543745
)
544746
kraus.append(k)
545747

546-
if not kraus:
547-
kraus.append(backend.zeros([output_dim, input_dim], dtype=dtypestr))
548-
return kraus
549-
550-
raise ValueError("illegal Choi matrix")
748+
if not kraus:
749+
kraus.append(backend.zeros([output_dim, input_dim], dtype=dtypestr))
750+
return kraus
551751

552752

553753
@partial(

‎tests/test_channels.py‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -78,3 +78,76 @@ def test_rep_transformation(backend):
7878
choi = np.zeros([4, 4])
7979
kraus = tc.channels.choi_to_kraus(choi)
8080
np.testing.assert_allclose(kraus, [np.zeros([2, 2])], atol=1e-5)
81+
82+
83+
@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")])
84+
def test_thermal(backend):
85+
t2 = 100
86+
time = 100
87+
88+
t1 = 180
89+
print(t1, t2)
90+
kraus = tc.channels.thermalrelaxationchannel(t1, t2, time, "general", 0.1)
91+
supop1 = tc.channels.kraus_to_super(kraus)
92+
# print(kraus)
93+
# print(supop1)
94+
95+
kraus = tc.channels.thermalrelaxationchannel(t1, t2, time, "T1dom", 0.1)
96+
supop2 = tc.channels.kraus_to_super(kraus)
97+
# print(kraus)
98+
# print(supop2)
99+
np.testing.assert_allclose(supop1, supop2, atol=1e-5)
100+
101+
t1 = 80
102+
print(t1, t2)
103+
kraus = tc.channels.thermalrelaxationchannel(t1, t2, time, "general", 0.1)
104+
supop1 = tc.channels.kraus_to_super(kraus)
105+
# print(kraus)
106+
# print(supop1)
107+
108+
kraus = tc.channels.thermalrelaxationchannel(t1, t2, time, "T2dom", 0.1)
109+
supop2 = tc.channels.kraus_to_super(kraus)
110+
# print(kraus)
111+
# print(supop2)
112+
np.testing.assert_allclose(supop1, supop2, atol=1e-5)
113+
114+
115+
@pytest.mark.parametrize("backend", [lf("npb"), lf("tfb"), lf("jaxb")])
116+
def test_noisecircuit(backend):
117+
118+
# Monte carlo simulation
119+
def noisecircuit(X):
120+
noise = tc.channels.thermalrelaxationchannel(300, 400, 1000, "T2dom", 0)
121+
122+
n = 1
123+
c = tc.Circuit(n)
124+
c.x(0)
125+
c.general_kraus(noise, 0, status=X)
126+
val = c.expectation_ps(z=[0])
127+
return val
128+
129+
noisec = noisecircuit
130+
noisec_vmap = tc.backend.vmap(noisec, vectorized_argnums=0)
131+
noisec_jit = tc.backend.jit(noisec_vmap)
132+
133+
nmc = 1000
134+
X = tc.backend.implicit_randu(nmc)
135+
valuemc = sum(tc.backend.numpy(noisec_jit(X))) / nmc
136+
137+
# Density matrix simulation
138+
def noisecircuitdm():
139+
n = 1
140+
dmc = tc.DMCircuit(n)
141+
dmc.x(0)
142+
dmc.thermalrelaxation(
143+
0, t1=300, t2=400, time=1000, method="T2dom", excitedstatepopulation=0
144+
)
145+
val = dmc.expectation_ps(z=[0])
146+
return val
147+
148+
noisec = noisecircuitdm
149+
noisec_jit = tc.backend.jit(noisec)
150+
151+
valuedm = noisec_jit()
152+
153+
np.testing.assert_allclose(valuemc, valuedm, atol=1e-1)

0 commit comments

Comments
 (0)