Unverified Commit 59d9f508 authored by dwuggh's avatar dwuggh
Browse files

update

add measure.py; update
parent 26cf40eb
......@@ -16,10 +16,7 @@ class DensityOperator(QOperator):
if mat is None:
mat = np.zeros((dim, dim))
mat[0, 0] = 1
# the following should be equivalent:
# ground_state = np.zeros(dim)
# ground_state[0] = 1
# mat = np.outer(ground_state, ground_state)
super().__init__(qubits, mat)
# def print(self):
......@@ -27,6 +24,11 @@ class DensityOperator(QOperator):
# print("shape: ", self.operator.shape)
# print(self.operator)
def deepcopy(self):
operator = self.operator
qubits = self.qubits
return DensityOperator(qubits, operator)
def merge(self, other):
op = multiply(self, other)
self.operator = op.operator
......@@ -39,6 +41,9 @@ class DensityOperator(QOperator):
op = self.broadcast_with(qubits)
self.operator = op.operator
self.qubits = op.qubits
def alter_qubits(self, qubits):
self.qubits = np.array(qubits)
def evolution(self, unitary: QOperator, p_g = 0, on_qubits = None):
......@@ -48,9 +53,6 @@ class DensityOperator(QOperator):
result = multiply(unitary, multiply(self, unitary.dagger()))
# big_unitary = unitary.broadcast(self.qnum())
# self.operator = np.matmul(np.matmul(big_unitary, self.operator), big_unitary.transpose())
self.operator = result.operator
self.qubits = result.qubits
......@@ -61,8 +63,6 @@ class DensityOperator(QOperator):
# TODO rename this function
def channel(self, channel: QChannel):
# qmax = np.max([self.qmax(), channel.qmax])
# result = sum(map(lambda E: np.matmul(np.matmul(E.operator, self.operator), E.dagger()), channel.kraus_operators))
result = reduce(add, (multiply(multiply(E, self), E.dagger()) for E in channel.kraus_operators))
self.operator = result.operator
self.qubits = result.qubits
......@@ -74,7 +74,7 @@ class DensityOperator(QOperator):
inversion of the state with probability $p_m$.
NOTE this function may cause porbability loss.
'''
def bell_measure(self, q1, q2, pauli1, pauli2, p_m = 0):
def bell_measure(self, q1, q2, pauli1, pauli2, p_m = 0, normal_probs = True):
channel1 = measure_x(q1) if pauli1 == 'x' else measure_z(q1)
channel2 = measure_x(q2) if pauli2 == 'x' else measure_z(q2)
......@@ -85,7 +85,7 @@ class DensityOperator(QOperator):
P_11 = multiply(channel1.kraus_operators[1], channel2.kraus_operators[1])
# P_00.print()
# the 4 resulting density matrix, with weight of p_ij
# the 4 result density matrix, with weight of p_ij
ρ_00 = multiply(self, P_00).partial_trace([q1, q2])
ρ_01 = multiply(self, P_01).partial_trace([q1, q2])
ρ_10 = multiply(self, P_10).partial_trace([q1, q2])
......@@ -103,14 +103,22 @@ class DensityOperator(QOperator):
p10 = p_m * (1 - p_m) * p_sum
p11 = (1 - p_m) ** 2 * p_11 + p_m ** 2 * p_00
# substitute ρ_00s probability
f = lambda x, y: 0 if y == 0 else x / y
p00 = f(p00, p_00)
p01 = f(p01, p_01)
p10 = f(p10, p_10)
p11 = f(p11, p_11)
# no rescaling: this is intentional
# p00 = p00 / p_sum
# p01 = p01 / p_sum
# p10 = p10 / p_sum
# p11 = p11 / p_sum
if normal_probs:
p00 = p00 / p_sum
p01 = p01 / p_sum
p10 = p10 / p_sum
p11 = p11 / p_sum
print(p_00, p_01, p_10, p_11)
print(p00 , p01 , p10 , p11)
# print(p_00, p_01, p_10, p_11)
# print(p00 , p01 , p10 , p11)
operator = ρ_00.operator * p00 + ρ_01.operator * p10 + ρ_10.operator * p10 +ρ_11.operator * p11
......@@ -133,11 +141,17 @@ def measure_z(q = 0):
return QChannel([QOperator([q], E_0), QOperator([q], E_1)])
# generate bell pair in werner form
def bell_pair(p_n, qubits = [0, 1]):
ρ_0 = np.array([
def bell_pair(p_n, qubits = [0, 1], i = 0):
ρ_00 = np.array([
[1, 0, 0, 1],
[0, 0, 0, 0],
[0, 0, 0, 0],
[1, 0, 0, 1]], dtype=np.float64) / 2
ρ_11 = np.array([
[1, 0, 0, -1],
[0, 0, 0, 0],
[0, 0, 0, 0],
[-1, 0, 0, 1]], dtype=np.float64) / 2
ρ_0 = ρ_00 if i == 0 else ρ_11
mat = (1 - 4 / 3 * p_n) * ρ_0 + (p_n / 3) * np.identity(4)
return DensityOperator(qubits, mat)
......@@ -27,7 +27,7 @@ class QChannel(object):
depolarizing channel for $n$ qubits
$$
ε(ρ) = (1 - p)ρ + \\frac{ρ}{4^n - 1} \sum (⊗ A_i) ρ (⊗ A_i^\dagger)
ε(ρ) = (1 - p)ρ + \\frac{ρ}{4^n - 1} \\sum (⊗ A_i) ρ (⊗ A_i^\\dagger)
$$
'''
def depolarizing_channel(p, qubits = [0]) -> QChannel:
......
......@@ -6,6 +6,11 @@ class QOperator(object):
self.qubits = np.array(qubits)
self.operator = np.array(operator, dtype=np.float64)
def deepcopy(self):
operator = self.operator
qubits = self.qubits
return QOperator(qubits, operator)
def qnum(self):
return self.qubits.size
......@@ -20,6 +25,7 @@ class QOperator(object):
def print(self):
print("qubits: ", self.qubits)
print("shape: ", self.operator.shape)
print("trace: ", self.operator.trace())
print(self.operator)
def alter_qubits(self, qubits):
......@@ -31,8 +37,10 @@ class QOperator(object):
operator = self.operator * scalar
return QOperator(self.qubits, operator)
# construct gate in the big hilbert space
# G = I ⊗ G
'''
construct gate in the big hilbert space
G = I ⊗ G
'''
def broadcast(self, qnum: int):
dim = 2 ** qnum
big_operator = np.zeros((dim, dim), dtype=np.float64)
......@@ -41,7 +49,7 @@ class QOperator(object):
di = np.flip(get_bin_digits(i, qnum))
dj = np.flip(get_bin_digits(j, qnum))
# coordinates in the small HIlbert space
# coordinates in the small Hilbert space
gi = from_bin_digits(np.flip(di[self.qubits]))
gj = from_bin_digits(np.flip(dj[self.qubits]))
......@@ -55,22 +63,17 @@ class QOperator(object):
big_operator[i, j] = val
return big_operator
'''
construct gate in the big hilbert space, with given qubits indices
G = I ⊗ G
'''
def broadcast_with(self, qubits):
# qubits_new = np.fromiter((x for x in qubits if not contains(self.qubits, x)), dtype=np.int32)
qnum = len(qubits)
dim = 2 ** qnum
big_operator = np.zeros((dim, dim), dtype=np.float64)
# occupied qubits of the original operator
# we need to respect the order of self.qubits in `qubits`
# occupy1: original qubits in parameter's order
# occupy1 = []
# for i, q in enumerate(qubits):
# if contains(self.qubits, q):
# occupy1.append(q)
# occupy: original qubits index in `self.qubits` of occupy1
# respect the order of self.qubits in `qubits`
occupy = []
for q1 in self.qubits:
for i, q2 in enumerate(qubits):
......@@ -78,7 +81,6 @@ class QOperator(object):
occupy.append(i)
break
# self.qubits = np.array(qubits)
for i in range(dim):
for j in range(dim):
di = np.flip(get_bin_digits(i, qnum))
......@@ -97,8 +99,6 @@ class QOperator(object):
val = np.float64(id) * self.operator[gi, gj]
big_operator[i, j] = val
return QOperator(qubits, big_operator)
# def rearrange(self, qubits_new):
# pass
# get adjoint operator
def dagger(self):
......@@ -142,77 +142,15 @@ class QOperator(object):
mat_new[i_new, j_new] += self.operator[i, j]
return QOperator(qubits_reserved, mat_new)
'''
extract qubits at given basis('x' or 'z')
'''
def extract(self, qubit, basis, target):
# get index of qubit
index = -1
for (i, q) in enumerate(self.qubits):
if q == qubit:
index = i
break
qnum = self.qnum()
qubits_new = np.delete(self.qubits, index)
dim, _ = self.operator.shape
# get sub matrices of |0><0| and |1><1| (Z basis)
dim_new = 2 ** (qnum - 1)
mat_00 = np.zeros((dim_new, dim_new), dtype=np.float64)
mat_01 = np.zeros((dim_new, dim_new), dtype=np.float64)
mat_10 = np.zeros((dim_new, dim_new), dtype=np.float64)
mat_11 = np.zeros((dim_new, dim_new), dtype=np.float64)
for i in range(dim):
for j in range(dim):
di = np.flip(get_bin_digits(i, qnum))
dj = np.flip(get_bin_digits(j, qnum))
ql = di[index]
qr = dj[index]
di_new = np.delete(di, index)
dj_new = np.delete(dj, index)
i_new = from_bin_digits(np.flip(di_new))
j_new = from_bin_digits(np.flip(dj_new))
if ql == 0 and qr == 0:
mat_00[i_new, j_new] = self.operator[i, j]
if ql == 0 and qr == 1:
mat_01[i_new, j_new] = self.operator[i, j]
if ql == 1 and qr == 0:
mat_10[i_new, j_new] = self.operator[i, j]
if ql == 1 and qr == 1:
mat_11[i_new, j_new] = self.operator[i, j]
operator = None
if basis == 'z':
if target == 0:
operator = mat_00
else:
operator = mat_11
else:
if target == 0:
mat_pp = (mat_00 + mat_01 + mat_10 + mat_11) / 2
operator = mat_pp
else:
mat_nn = (mat_00 - mat_01 - mat_10 + mat_11) / 2
operator = mat_nn
return QOperator(qubits_new, operator)
# multiply 2 QOperators, often needs to broadcast them
def multiply(lhs: QOperator, rhs: QOperator):
q1 = lhs.qubits
q2 = rhs.qubits
qubits = merge(q1, q2)
# print(q1, q2, qubits)
lhs_big = lhs.broadcast_with(qubits)
# print("lhs", lhs.operator.trace())
rhs_big = rhs.broadcast_with(qubits)
# print("rhs", rhs.operator.trace())
result_mat = np.matmul(lhs_big.operator, rhs_big.operator)
# print("result", result_mat.trace())
return QOperator(qubits, result_mat)
def add(lhs: QOperator, rhs: QOperator):
......
import numpy as np
from copy import deepcopy
from functools import reduce
from .utils import *
from .QOperator import *
......@@ -9,6 +10,16 @@ from .fidelity import *
'''
notation in circuit:
q2 ---∩---
q1 ---‖---
data1 ---σ---
|
q4 ---∩---
q3 ---‖---
data2 ---σ---
q1, q2, q3, q4 are the 4 ancilla qubits
q2 ---------∎--⊤--M(x)
q1 ---∎--⊤--|--Z--M(x)
......@@ -41,12 +52,23 @@ def bell_purify_1(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1,
ρ.evolution(c4, err_model.p_g)
# measurement
ρ.bell_measure(q1, q3, 'x', err_model.p_m)
ρ.bell_measure(q2, q4, 'x', err_model.p_m)
ρ.bell_measure(q1, q3, 'x', 'x', err_model.p_m, True)
ρ.bell_measure(q2, q4, 'x', 'x', err_model.p_m, True)
return ρ
'''
notation in circuit:
q2 ---∩---
q1 ---∎---
data1 ---σ---
|
q4 ---∩---
q3 ---∎---
data2 ---σ---
q1, q2, q3, q4 are the 4 ancilla qubits
`stringent` indicates whether the blank-separated region is performed
q2 ------∎--⊤--M(x)--∎--⊤--M(x)---- --∎--⊤--M(x)-- ------
......@@ -57,11 +79,13 @@ q4 ---|--∎--⊤--M(x)--∎--⊤--M(x)---- --∎--⊤--M(x)-- ------
q3 ---∎-----X-----------Z-------⊤-- -----Z-------- --M(x)
data2 -----------------------------σ-- -------------- ------
'''
def bell_purify_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1, q2, q3 ,q4, pauli, stringent = True):
def bell_purify_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1, q2, q3 ,q4, pauli, stringent = True, debug = False):
# for efficiency improvement
ρ1 = bell_pair(err_model.p_n, [q1, q3])
ρ2 = bell_pair(err_model.p_n, [q2, q4])
# ρ2 = bell_pair(err_model.p_n, [q2, q4])
ρ2 = ρ1.deepcopy()
ρ2.alter_qubits([q2, q4])
ρ1.merge(ρ2)
......@@ -70,7 +94,7 @@ def bell_purify_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1,
ρ1.evolution(c1, err_model.p_g)
ρ1.evolution(c2, err_model.p_g)
ρ1.bell_measure(q2, q4, 'x', err_model.p_m)
ρ1.bell_measure(q2, q4, 'x', 'x', err_model.p_m)
ρ3 = bell_pair(err_model.p_n, [q2, q4])
ρ1.merge(ρ3)
......@@ -80,7 +104,10 @@ def bell_purify_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1,
ρ1.evolution(c3, err_model.p_g)
ρ1.evolution(c4, err_model.p_g)
ρ1.bell_measure(q2, q4, 'x', err_model.p_m)
ρ1.bell_measure(q2, q4, 'x', 'x', err_model.p_m)
if debug:
ρ1.print()
# now, take ρ into consideration
ρ.merge(ρ1)
......@@ -99,11 +126,9 @@ def bell_purify_2(ρ: DensityOperator, err_model: ErrorModel, data1, data2, q1,
ρ.evolution(c7, err_model.p_g)
ρ.evolution(c8, err_model.p_g)
ρ.bell_measure(q2, q4, 'x', err_model.p_m)
pass
ρ.bell_measure(q2, q4, 'x', 'x', err_model.p_m)
ρ.bell_measure(q1, q3, 'x', err_model.p_m)
ρ.bell_measure(q1, q3, 'x', 'x', err_model.p_m)
'''
qubit indexing:
......@@ -111,27 +136,88 @@ qubit indexing:
1 4
0 3
'''
def make_bell(err_model: ErrorModel, stringent = True, stringent_plus = True):
def make_bell(err_model: ErrorModel, stringent = True, stringent_plus = True, step = False):
ρ = bell_pair(err_model.p_n, [0, 3])
bell_purify_1(ρ, err_model, 0, 3, 1, 2, 4, 5, 'z')
print(bell_fidelity(ρ))
if step:
print(bell_fidelity(ρ))
bell_purify_1(ρ, err_model, 0, 3, 1, 2, 4, 5, 'x')
print(bell_fidelity(ρ))
if step:
print(bell_fidelity(ρ))
if stringent:
bell_purify_2(ρ, err_model, 0, 3, 1, 2, 4, 5, 'z', stringent)
print(bell_fidelity(ρ))
if step:
print(bell_fidelity(ρ))
bell_purify_2(ρ, err_model, 0, 3, 1, 2, 4, 5, 'x', stringent)
print(bell_fidelity(ρ))
if step:
print(bell_fidelity(ρ))
print("final bell fidelity in make_bell", bell_fidelity(ρ))
return ρ
def make_GHZ(ρ1: DensityOperator, ρ2: DensityOperator, err_model: ErrorModel, stringent = True, stringent_plus = True):
pass
'''
qubit indexing:
2 5
1 4
0 3
---- data qubits
12 13
def stablizer_measurement(err_model: ErrorModel, stringent = True, stringent_plus = True):
14 15
---- data qubits
6 9
7 10
8 11
circuit:
q2 ---∩--∩--
q1 ---∎--∎--
data1 ---z--z--
| |
q4 ---∩--∩--
q3 ---∎--∎--
data2 ---z--z--
'''
def make_GHZ(err_model: ErrorModel, stringent = True, stringent_plus = True):
ρ1 = make_bell(err_model, stringent, stringent_plus)
ρ2 = make_bell(err_model, stringent, stringent_plus)
ρ2 = ρ2.alter_qubits([6, 9])
return make_GHZ_with(ρ1, err_model, stringent, stringent_plus)
def make_GHZ_with(ρ1: DensityOperator, err_model: ErrorModel, stringent = True, stringent_plus = True):
# ρ1 = make_bell(err_model, stringent, stringent_plus)
# ρ2 = make_bell(err_model, stringent, stringent_plus)
ρ2= ρ1.deepcopy()
ρ1.alter_qubits([0, 3])
ρ2.alter_qubits([6, 9])
ρ1.merge(ρ2)
bell_purify_2(ρ1, err_model, 0, 6, 1, 2, 7, 8, 'z', stringent)
# FIXME qubits mismatch for the next line, if we choose ancillas' indices to be 4 5 10 11,
# then the result density matrix would consist of 5 qubits of indices 0, 1, 3, 6, 9.
# However, if the ancilla indices are 1 2 7 8, everything works fine, at least seemingly.
# Don't know why, not even a clue.
# bell_purify_2(ρ1, err_model, 3, 9, 4, 5, 10, 11, 'z', stringent)
bell_purify_2(ρ1, err_model, 3, 9, 1, 2, 7, 8, 'z', stringent)
return ρ1
def stabilizer_measurement(err_model: ErrorModel, stabilizer, stringent = True, stringent_plus = True):
ρ = make_GHZ(err_model, stringent, stringent_plus)
return stabilizer_measurement_with(ρ, err_model, stabilizer, stringent, stringent_plus)
'''
σ is Z or X, indicating stabilizer.
circuit:
GHZ ---⊤--M(x)
data ---σ------
'''
def stabilizer_measurement_with(ρ: DensityOperator, err_model: ErrorModel, stabilizer, stringent = True, stringent_plus = True):
pass
......@@ -158,8 +244,6 @@ q4 ---|--∎--⊤--M(x)
q3 ---∎-----X------
'''
def BBPSSW(err_model: ErrorModel, detail = False):
# ρ1 = bell_pair(err_model.p_n, [1, 3])
# ρ2 = bell_pair(err_model.p_n, [2, 4])
ρ1 = BBPSSW_bell_pair([1, 3], err_model.p_n)
ρ2 = BBPSSW_bell_pair([2, 4], err_model.p_n)
if detail:
......@@ -171,13 +255,10 @@ def BBPSSW(err_model: ErrorModel, detail = False):
ρ1.print()
c1 = cnot([1, 2])
# c1 = cnot([2, 1])
c2 = cnot([4, 3])
c2 = multiply(c1, c2)
# c2.broadcast_with([1, 3, 2, 4]).print()
# ρ1.evolution(c1, err_model.p_g)
ρ1.evolution(c2, err_model.p_g)
ρ1.broadcast_with_self([1, 3, 2, 4])
......@@ -185,12 +266,7 @@ def BBPSSW(err_model: ErrorModel, detail = False):
if detail:
ρ1.print()
# ρ1.bell_measure(2, 4, 'z')
probs = ρ1.operator.trace()
if detail:
print("probs", probs)
ρ1.bell_measure(2, 4, 'z', 'x')
# ρ1 = ρ1.partial_trace([2, 4])
ρ1.bell_measure(2, 4, 'z', 'x', err_model.p_m, False)
probs = ρ1.operator.trace()
ρ1.scale(1 / probs)
......@@ -209,8 +285,8 @@ q4 ---|--∎--⊤--M(z)
q3 ---∎-----X------
'''
def BBPSSW_2(err_model: ErrorModel, detail = False):
ρ1 = bell_pair(err_model.p_n, [1, 3])
ρ2 = bell_pair(err_model.p_n, [2, 4])
ρ1 = bell_pair(err_model.p_n, [1, 3], 3)
ρ2 = bell_pair(err_model.p_n, [2, 4], 3)
if detail:
ρ1.print()
......@@ -234,8 +310,7 @@ def BBPSSW_2(err_model: ErrorModel, detail = False):
if detail:
ρ1.print()
# ρ1.bell_measure(2, 4, 'z')
ρ1.bell_measure(2, 4, 'z', 'z')
ρ1.bell_measure(2, 4, 'z', 'z', err_model.p_m, False)
# ρ1 = ρ1.partial_trace([2, 4])
probs = ρ1.operator.trace()
ρ1.scale(1 / probs)
......
......@@ -7,10 +7,12 @@ if __name__ == "__main__":
np.set_printoptions(edgeitems=16, linewidth=200,
formatter=dict(float=lambda x: "%6.4g" % x)
)
perfect_bell = cc.bell_pair(0, [0, 1])
err_model = cc.ErrorModel(0.05, 0.005, 0.005)
noise_bell = cc.bell_pair(err_model.p_n, [0, 1])
print("no purification: ", cc.entanglement_fidelity(perfect_bell, noise_bell))
ρ = cc.make_bell(err_model, False, False)
print("with purification: ", cc.entanglement_fidelity(perfect_bell, ρ))
err_model = cc.ErrorModel(0.1, 0.006, 0.006)
stringent = False
bell = cc.make_bell(err_model, stringent, False, True)
bell_2 = bell.deepcopy()
bell_2.alter_qubits([6, 9])
bell.merge(bell_2)
cc.bell_purify_2(bell, err_model, 0, 6, 1, 2, 7, 8, 'z', stringent)
cc.bell_purify_2(bell, err_model, 3, 9, 1, 2, 7, 8, 'z', stringent)
bell.print()
import numpy as np
import calc_channel as cc
GHZ = [
[0.4812, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4596],
[0, 0.006503, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0004135, 0],
[0, 0, 0.0008917, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0004135, 0, 0],
[0, 0, 0, 0.0001813, 0, 0, 0, 0, 0, 0, 0, 0, 3.433e-05, 0, 0, 0],
[0, 0, 0, 0, 0.001899, 0, 0, 0, 0, 0, 0, 0.0005545, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0.00805, 0, 0, 0, 0, 0.006679, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 5.584e-05, 0, 0, 4.989e-07, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0.001171, 0.0005545, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0.0005545, 0.001171, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 4.989e-07, 0, 0, 5.584e-05, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0.006679, 0, 0, 0, 0, 0.00805, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0.0005545, 0, 0, 0, 0, 0, 0, 0.001899, 0, 0, 0, 0],
[0, 0, 0, 3.433e-05, 0, 0, 0, 0, 0, 0, 0, 0, 0.0001813, 0, 0, 0],
[0, 0, 0.0004135, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0008917, 0, 0],
[0, 0.0004135, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.006503, 0],
[0.4596, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4812],
]
# 0.940876
GHZ = [
[0.4864, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.4406],
[0, 0.004329, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0003899, 0],
[0, 0, 0.0008959, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0003899, 0, 0],
[0, 0, 0, 0.00034, 0, 0, 0, 0, 0, 0, 0, 0, 0.0001047, 0, 0, 0],
[0, 0, 0, 0, 0.001844, 0, 0, 0, 0, 0, 0, 0.0007361, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0.004601, 0, 0, 0, 0, 0.00274, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 5.125e-05, 0, 0, 6.513e-07, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0.001503, 0.0007361, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0.0007361, 0.001503, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 6.513e-07, 0, 0, 5.125e-05, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0.00274, 0, 0, 0, 0, 0.004601, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0.0007361, 0, 0, 0, 0, 0, 0, 0.001844, 0, 0, 0, 0],
[0, 0, 0, 0.0001047, 0, 0, 0, 0, 0, 0, 0, 0, 0.00034, 0, 0, 0],
[0, 0, 0.0003899,