# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import itertools
from typing import cast, Tuple
import cirq
import cirq.contrib.acquaintance as cca
import numpy as np
import openfermion
import pytest
import scipy.linalg as la
import sympy
import openfermioncirq as ofc
from openfermioncirq.gates.fermionic_simulation import (
sum_of_interaction_operator_gate_generators,
state_swap_eigen_component,
)
def test_state_swap_eigen_component_args():
with pytest.raises(TypeError):
state_swap_eigen_component(0, '12', 1)
with pytest.raises(ValueError):
state_swap_eigen_component('01', '01', 1)
with pytest.raises(ValueError):
state_swap_eigen_component('01', '10', 0)
with pytest.raises(ValueError):
state_swap_eigen_component('01', '100', 1)
with pytest.raises(ValueError):
state_swap_eigen_component('01', 'ab', 1)
@pytest.mark.parametrize('index_pair,n_qubits', [
((0, 1), 2),
((0, 3), 2),
])
def test_state_swap_eigen_component(index_pair, n_qubits):
state_pair = tuple(format(i, '0' + str(n_qubits) + 'b') for i in index_pair)
i, j = index_pair
dim = 2**n_qubits
for sign in (-1, 1):
actual_component = state_swap_eigen_component(state_pair[0],
state_pair[1], sign)
expected_component = np.zeros((dim, dim))
expected_component[i, i] = expected_component[j, j] = 0.5
expected_component[i, j] = expected_component[j, i] = sign * 0.5
assert np.allclose(actual_component, expected_component)
@pytest.mark.parametrize('n_modes, seed',
[(7, np.random.randint(1 << 30)) for _ in range(2)])
def test_interaction_operator_interconversion(n_modes, seed):
operator = openfermion.random_interaction_operator(n_modes,
real=False,
seed=seed)
gates = ofc.fermionic_simulation_gates_from_interaction_operator(operator)
other_operator = sum_of_interaction_operator_gate_generators(n_modes, gates)
operator = openfermion.normal_ordered(operator)
other_operator = openfermion.normal_ordered(other_operator)
assert operator == other_operator
def test_interaction_operator_from_bad_gates():
for gates in [{(): 'bad'}, {(0,): cirq.X}]:
with pytest.raises(TypeError):
sum_of_interaction_operator_gate_generators(5, gates)
def random_real(size=None, mag=20):
return np.random.uniform(-mag, mag, size)
def random_complex(size=None, mag=20):
return random_real(size, mag) + 1j * random_real(size, mag)
def random_fermionic_simulation_gate(order):
exponent = random_real()
if order == 2:
weights = (random_complex(), random_real())
return ofc.QuadraticFermionicSimulationGate(weights, exponent=exponent)
weights = random_complex(3)
if order == 3:
return ofc.CubicFermionicSimulationGate(weights, exponent=exponent)
if order == 4:
return ofc.QuarticFermionicSimulationGate(weights, exponent=exponent)
def assert_symbolic_decomposition_consistent(gate):
expected_unitary = cirq.unitary(gate)
weights = tuple(sympy.Symbol(f'w{i}') for i in range(gate.num_weights()))
exponent = sympy.Symbol('t')
symbolic_gate = type(gate)(weights, exponent=exponent)
qubits = cirq.LineQubit.range(gate.num_qubits())
circuit = cirq.Circuit(symbolic_gate._decompose_(qubits))
resolver = {'t': gate.exponent}
for i, w in enumerate(gate.weights):
resolver[f'w{i}'] = w
resolved_circuit = cirq.resolve_parameters(circuit, resolver)
decomp_unitary = resolved_circuit.unitary(qubit_order=qubits)
assert np.allclose(expected_unitary, decomp_unitary)
def assert_generators_consistent(gate):
qubit_generator = gate.qubit_generator_matrix
qubit_generator_from_fermion_generator = (super(
type(gate), gate).qubit_generator_matrix)
assert np.allclose(qubit_generator, qubit_generator_from_fermion_generator)
def assert_resolution_consistent(gate):
weight_names = [f'w{i}' for i in range(gate.num_weights())]
weight_params = [sympy.Symbol(w) for w in weight_names]
resolver = dict(zip(weight_names, gate.weights))
resolver['s'] = gate._global_shift
resolver['e'] = gate._exponent
symbolic_gate = type(gate)(weight_params,
exponent=sympy.Symbol('e'),
global_shift=sympy.Symbol('s'))
resolved_gate = cirq.resolve_parameters(symbolic_gate, resolver)
assert resolved_gate == gate
def assert_fswap_consistent(gate):
gate = gate.__copy__()
n_qubits = gate.num_qubits()
for i in range(n_qubits - 1):
fswap = cirq.kron(np.eye(1 << i), cirq.unitary(ofc.FSWAP),
np.eye(1 << (n_qubits - i - 2)))
assert fswap.shape == (1 << n_qubits,) * 2
generator = gate.qubit_generator_matrix
fswapped_generator = np.linalg.multi_dot([fswap, generator, fswap])
gate.fswap(i)
assert np.allclose(gate.qubit_generator_matrix, fswapped_generator)
for i in (-1, n_qubits):
with pytest.raises(ValueError):
gate.fswap(i)
def assert_permute_consistent(gate):
gate = gate.__copy__()
n_qubits = gate.num_qubits()
qubits = cirq.LineQubit.range(n_qubits)
for pos in itertools.permutations(range(n_qubits)):
permuted_gate = gate.__copy__()
gate.permute(pos)
assert permuted_gate.permuted(pos) == gate
actual_unitary = cirq.unitary(permuted_gate)
ops = [
cca.LinearPermutationGate(n_qubits, dict(zip(range(n_qubits), pos)),
ofc.FSWAP)(*qubits),
gate(*qubits),
cca.LinearPermutationGate(n_qubits, dict(zip(pos, range(n_qubits))),
ofc.FSWAP)(*qubits)
]
circuit = cirq.Circuit(ops)
expected_unitary = cirq.unitary(circuit)
assert np.allclose(actual_unitary, expected_unitary)
with pytest.raises(ValueError):
gate.permute(range(1, n_qubits))
with pytest.raises(ValueError):
gate.permute([1] * n_qubits)
def assert_interaction_operator_consistent(gate):
interaction_op = gate.interaction_operator_generator()
other_gate = gate.from_interaction_operator(operator=interaction_op)
if other_gate is None:
assert np.allclose(gate.weights, 0)
else:
assert cirq.approx_eq(gate, other_gate)
interaction_op = openfermion.normal_ordered(interaction_op)
other_interaction_op = openfermion.InteractionOperator.zero(
interaction_op.n_qubits)
super(type(gate),
gate).interaction_operator_generator(operator=other_interaction_op)
other_interaction_op = openfermion.normal_ordered(interaction_op)
assert interaction_op == other_interaction_op
other_interaction_op = super(type(gate),
gate).interaction_operator_generator()
other_interaction_op = openfermion.normal_ordered(interaction_op)
assert interaction_op == other_interaction_op
random_quadratic_gates = [random_fermionic_simulation_gate(2) for _ in range(5)]
manual_quadratic_gates = [
ofc.QuadraticFermionicSimulationGate(weights)
for weights in [cast(Tuple[float, float], (1, 1)), (1, 0), (0, 1), (0, 0)]
]
quadratic_gates = random_quadratic_gates + manual_quadratic_gates
cubic_gates = ([ofc.CubicFermionicSimulationGate()] +
[random_fermionic_simulation_gate(3) for _ in range(5)])
quartic_gates = ([ofc.QuarticFermionicSimulationGate()] +
[random_fermionic_simulation_gate(4) for _ in range(5)])
gates = quadratic_gates + cubic_gates + quartic_gates
@pytest.mark.parametrize('gate', gates)
def test_fermionic_simulation_gate(gate):
ofc.testing.assert_implements_consistent_protocols(gate)
generator = gate.qubit_generator_matrix
expected_unitary = la.expm(-1j * gate.exponent * generator)
actual_unitary = cirq.unitary(gate)
assert np.allclose(expected_unitary, actual_unitary)
assert_fswap_consistent(gate)
assert_permute_consistent(gate)
assert_generators_consistent(gate)
assert_resolution_consistent(gate)
assert_interaction_operator_consistent(gate)
assert gate.num_weights() == super(type(gate), gate).num_weights()
@pytest.mark.parametrize('weights', list(np.random.rand(10, 3)) + [(1, 0, 1)])
def test_weights_and_exponent(weights):
exponents = np.linspace(-1, 1, 8)
gates = tuple(
ofc.QuarticFermionicSimulationGate(
weights / exponent, exponent=exponent, absorb_exponent=True)
for exponent in exponents)
for g1, g2 in itertools.combinations(gates, 2):
assert cirq.approx_eq(g1, g2, atol=1e-100)
for i, (gate, exponent) in enumerate(zip(gates, exponents)):
assert gate.exponent == 1
new_exponent = exponents[-i]
new_gate = gate._with_exponent(new_exponent)
assert new_gate.exponent == new_exponent
def test_zero_weights():
for gate_type in [
ofc.QuadraticFermionicSimulationGate,
ofc.CubicFermionicSimulationGate, ofc.QuarticFermionicSimulationGate
]:
weights = (0,) * gate_type.num_weights()
gate = gate_type(weights)
n_qubits = gate.num_qubits()
assert np.allclose(cirq.unitary(gate), np.eye(2**n_qubits))
cirq.testing.assert_decompose_is_consistent_with_unitary(gate)
operator = openfermion.InteractionOperator.zero(n_qubits)
assert gate_type.from_interaction_operator(operator=operator) is None
@pytest.mark.parametrize(
'weights,exponent',
[((np.random.uniform(-5, 5) + 1j * np.random.uniform(-5, 5),
np.random.uniform(-5, 5)), np.random.uniform(-5, 5)) for _ in range(5)])
def test_quadratic_fermionic_simulation_gate_unitary(weights, exponent):
generator = np.zeros((4, 4), dtype=np.complex128)
# w0 |10><01| + h.c.
generator[2, 1] = weights[0]
generator[1, 2] = weights[0].conjugate()
# w1 |11><11|
generator[3, 3] = weights[1]
expected_unitary = la.expm(-1j * exponent * generator)
gate = ofc.QuadraticFermionicSimulationGate(weights, exponent=exponent)
actual_unitary = cirq.unitary(gate)
assert np.allclose(expected_unitary, actual_unitary)
symbolic_gate = (ofc.QuadraticFermionicSimulationGate(
(sympy.Symbol('w0'), sympy.Symbol('w1')), exponent=sympy.Symbol('t')))
qubits = cirq.LineQubit.range(2)
circuit = cirq.Circuit(symbolic_gate._decompose_(qubits))
resolver = {'w0': weights[0], 'w1': weights[1], 't': exponent}
resolved_circuit = cirq.resolve_parameters(circuit, resolver)
decomp_unitary = resolved_circuit.unitary(qubit_order=qubits)
assert np.allclose(expected_unitary, decomp_unitary)
@pytest.mark.parametrize('gate', random_quadratic_gates)
def test_quadratic_fermionic_simulation_gate_symbolic_decompose(gate):
assert_symbolic_decomposition_consistent(gate)
def test_cubic_fermionic_simulation_gate_equality():
eq = cirq.testing.EqualsTester()
eq.add_equality_group(
ofc.CubicFermionicSimulationGate()**0.5,
ofc.CubicFermionicSimulationGate((1,) * 3, exponent=0.5),
ofc.CubicFermionicSimulationGate((0.5,) * 3),
)
eq.add_equality_group(ofc.CubicFermionicSimulationGate((1j, 0, 0)),)
eq.add_equality_group(
ofc.CubicFermionicSimulationGate((sympy.Symbol('s'), 0, 0), exponent=2),
ofc.CubicFermionicSimulationGate((2 * sympy.Symbol('s'), 0, 0),
exponent=1),
)
eq.add_equality_group(
ofc.CubicFermionicSimulationGate((0, 0.7, 0), global_shift=2),
ofc.CubicFermionicSimulationGate((0, 0.35, 0),
global_shift=1,
exponent=2),
)
eq.add_equality_group(
ofc.CubicFermionicSimulationGate((1, 1, 1)),
ofc.CubicFermionicSimulationGate(((1 + 2 * np.pi), 1, 1)),
)
@pytest.mark.parametrize('exponent,control',
itertools.product([0, 1, -1, 0.25, -0.5, 0.1],
[0, 1, 2]))
def test_cubic_fermionic_simulation_gate_consistency_special(exponent, control):
weights = tuple(np.eye(1, 3, control)[0] * 0.5 * np.pi)
general_gate = ofc.CubicFermionicSimulationGate(weights, exponent=exponent)
general_unitary = cirq.unitary(general_gate)
indices = np.dot(list(itertools.product((0, 1), repeat=3)),
(2**np.roll(np.arange(3), -control))[::-1])
special_gate = cirq.ControlledGate(cirq.ISWAP**-exponent)
special_unitary = (
cirq.unitary(special_gate)[indices[:, np.newaxis], indices])
assert np.allclose(general_unitary, special_unitary)
@pytest.mark.parametrize(
'weights,exponent',
[(np.random.uniform(-5, 5, 3) + 1j * np.random.uniform(-5, 5, 3),
np.random.uniform(-5, 5)) for _ in range(5)])
def test_cubic_fermionic_simulation_gate_consistency_docstring(
weights, exponent):
generator = np.zeros((8, 8), dtype=np.complex128)
# w0 |110><101| + h.c.
generator[6, 5] = weights[0]
generator[5, 6] = weights[0].conjugate()
# w1 |110><011| + h.c.
generator[6, 3] = weights[1]
generator[3, 6] = weights[1].conjugate()
# w2 |101><011| + h.c.
generator[5, 3] = weights[2]
generator[3, 5] = weights[2].conjugate()
expected_unitary = la.expm(-1j * exponent * generator)
gate = ofc.CubicFermionicSimulationGate(weights, exponent=exponent)
actual_unitary = cirq.unitary(gate)
assert np.allclose(expected_unitary, actual_unitary)
def test_quartic_fermionic_simulation_consistency():
ofc.testing.assert_implements_consistent_protocols(
ofc.QuarticFermionicSimulationGate())
quartic_fermionic_simulation_simulator_test_cases = [
(ofc.QuarticFermionicSimulationGate(
(0, 0, 0)), 1., np.ones(16) / 4., np.ones(16) / 4., 5e-6),
(ofc.QuarticFermionicSimulationGate((0.2, -0.1, 0.7)), 0.,
np.array([1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) / 4.,
np.array([1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) / 4.,
5e-6),
(ofc.QuarticFermionicSimulationGate((0.2, -0.1, 0.7)), 0.3,
np.array([1, -1, -1, -1, -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]) / 4.,
np.array([
1, -1, -1, -np.exp(0.21j), -1, -np.exp(-0.03j),
np.exp(-0.06j), 1, 1,
np.exp(-0.06j),
np.exp(-0.03j), 1,
np.exp(0.21j), 1, 1, 1
]) / 4., 5e-6),
(ofc.QuarticFermionicSimulationGate((1. / 3, 0, 0)), 1.,
np.array([0, 0, 0, 0, 0, 0, 1., 0, 0, 1., 0, 0, 0, 0, 0, 0]) / np.sqrt(2),
np.array([0, 0, 0, 0, 0, 0, 1., 0, 0, 1., 0, 0, 0, 0, 0, 0]) / np.sqrt(2),
5e-6),
(ofc.QuarticFermionicSimulationGate((0, np.pi / 3, 0)), 1.,
np.array([1., 1., 0, 0, 0, 1., 0, 0, 0, 0., -1., 0, 0, 0, 0, 0]) / 2.,
np.array([
1., 1., 0, 0, 0, -np.exp(4j * np.pi / 3), 0, 0, 0, 0.,
-np.exp(1j * np.pi / 3), 0, 0, 0, 0, 0
]) / 2., 5e-6),
(ofc.QuarticFermionicSimulationGate((0, 0, -np.pi / 2)), 1.,
np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1., 0, 0,
0]), np.array([0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0]), 5e-6),
(ofc.QuarticFermionicSimulationGate((0, 0, -0.25 * np.pi)), 1.,
np.array([0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]),
np.array([0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1j, 0, 0, 0]) / np.sqrt(2),
5e-6),
(ofc.QuarticFermionicSimulationGate(
(-np.pi / 4, np.pi / 6, -np.pi / 2)), 1.,
np.array([0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0]) / np.sqrt(3),
np.array([
0, 0, 0, 1j, 0, -1j / 2., 1 / np.sqrt(2), 0, 0, 1j / np.sqrt(2),
np.sqrt(3) / 2, 0, 0, 0, 0, 0
]) / np.sqrt(3), 5e-6),
]
@pytest.mark.parametrize('gate, exponent, initial_state, correct_state, atol',
quartic_fermionic_simulation_simulator_test_cases)
def test_quartic_fermionic_simulation_on_simulator(gate, exponent,
initial_state, correct_state,
atol):
a, b, c, d = cirq.LineQubit.range(4)
circuit = cirq.Circuit(gate(a, b, c, d)**exponent)
result = circuit.final_wavefunction(initial_state)
cirq.testing.assert_allclose_up_to_global_phase(result,
correct_state,
atol=atol)
def test_quartic_fermionic_simulation_eq():
eq = cirq.testing.EqualsTester()
eq.add_equality_group(
ofc.QuarticFermionicSimulationGate((1.2, 0.4, -0.4), exponent=0.5),
ofc.QuarticFermionicSimulationGate((0.3, 0.1, -0.1), exponent=2),
ofc.QuarticFermionicSimulationGate((-0.6, -0.2, 0.2), exponent=-1),
ofc.QuarticFermionicSimulationGate((0.6, 0.2, 2 * np.pi - 0.2)),
)
eq.add_equality_group(
ofc.QuarticFermionicSimulationGate((-0.6, 0.0, 0.3), exponent=0.5))
eq.make_equality_group(lambda: ofc.QuarticFermionicSimulationGate(
(0.1, -0.3, 0.0), exponent=0.0))
eq.make_equality_group(lambda: ofc.QuarticFermionicSimulationGate(
(1., -1., 0.5), exponent=0.75))
def test_quadratic_fermionic_simulation_gate_text_diagram():
gate = ofc.QuadraticFermionicSimulationGate((1, 1))
a, b, c = cirq.LineQubit.range(3)
circuit = cirq.Circuit([gate(a, b), gate(b, c)])
assert super(type(gate), gate).wire_symbol(False) == type(gate).__name__
assert (super(type(gate), gate)._diagram_exponent(
cirq.CircuitDiagramInfoArgs.UNINFORMED_DEFAULT) == gate._exponent)
expected_text_diagram = """
0: ───↓↑(1, 1)──────────────
│
1: ───↓↑─────────↓↑(1, 1)───
│
2: ──────────────↓↑─────────
""".strip()
cirq.testing.assert_has_diagram(circuit, expected_text_diagram)
expected_text_diagram = """
0: ---a*a(1, 1)---------------
|
1: ---a*a---------a*a(1, 1)---
|
2: ---------------a*a---------
""".strip()
cirq.testing.assert_has_diagram(circuit,
expected_text_diagram,
use_unicode_characters=False)
def test_cubic_fermionic_simulation_gate_text_diagram():
gate = ofc.CubicFermionicSimulationGate((1, 1, 1))
qubits = cirq.LineQubit.range(5)
circuit = cirq.Circuit([gate(*qubits[:3]), gate(*qubits[2:5])])
assert super(type(gate), gate).wire_symbol(False) == type(gate).__name__
assert (super(type(gate), gate)._diagram_exponent(
cirq.CircuitDiagramInfoArgs.UNINFORMED_DEFAULT) == gate._exponent)
expected_text_diagram = """
0: ───↕↓↑(1, 1, 1)──────────────────
│
1: ───↕↓↑───────────────────────────
│
2: ───↕↓↑────────────↕↓↑(1, 1, 1)───
│
3: ──────────────────↕↓↑────────────
│
4: ──────────────────↕↓↑────────────
""".strip()
cirq.testing.assert_has_diagram(circuit, expected_text_diagram)
expected_text_diagram = """
0: ---na*a(1, 1, 1)-------------------
|
1: ---na*a----------------------------
|
2: ---na*a------------na*a(1, 1, 1)---
|
3: -------------------na*a------------
|
4: -------------------na*a------------
""".strip()
cirq.testing.assert_has_diagram(circuit,
expected_text_diagram,
use_unicode_characters=False)
test_weights = [1.0, 0.5, 0.25, 0.1, 0.0, -0.5]
@pytest.mark.parametrize('weights',
itertools.chain(
itertools.product(test_weights, repeat=3),
np.random.rand(10, 3)))
def test_quartic_fermionic_simulation_decompose(weights):
cirq.testing.assert_decompose_is_consistent_with_unitary(
ofc.QuarticFermionicSimulationGate(weights))
@pytest.mark.parametrize(
'weights,exponent',
[(np.random.uniform(-5, 5, 3) + 1j * np.random.uniform(-5, 5, 3),
np.random.uniform(-5, 5)) for _ in range(5)])
def test_quartic_fermionic_simulation_unitary(weights, exponent):
generator = np.zeros((1 << 4,) * 2, dtype=np.complex128)
# w0 |1001><0110| + h.c.
generator[9, 6] = weights[0]
generator[6, 9] = weights[0].conjugate()
# w1 |1010><0101| + h.c.
generator[10, 5] = weights[1]
generator[5, 10] = weights[1].conjugate()
# w2 |1100><0011| + h.c.
generator[12, 3] = weights[2]
generator[3, 12] = weights[2].conjugate()
expected_unitary = la.expm(-1j * exponent * generator)
gate = ofc.QuarticFermionicSimulationGate(weights, exponent=exponent)
actual_unitary = cirq.unitary(gate)
assert np.allclose(expected_unitary, actual_unitary)
def test_quartic_fermionic_simulation_gate_text_diagram():
gate = ofc.QuarticFermionicSimulationGate((1, 1, 1))
qubits = cirq.LineQubit.range(6)
circuit = cirq.Circuit([gate(*qubits[:4]), gate(*qubits[-4:])])
assert super(type(gate), gate).wire_symbol(False) == type(gate).__name__
for G in (gate, gate._with_exponent('e')):
assert (super(type(G), G)._diagram_exponent(
cirq.CircuitDiagramInfoArgs.UNINFORMED_DEFAULT) == G._exponent)
expected_text_diagram = """
0: ───⇊⇈(1, 1, 1)─────────────────
│
1: ───⇊⇈──────────────────────────
│
2: ───⇊⇈────────────⇊⇈(1, 1, 1)───
│ │
3: ───⇊⇈────────────⇊⇈────────────
│
4: ─────────────────⇊⇈────────────
│
5: ─────────────────⇊⇈────────────
""".strip()
cirq.testing.assert_has_diagram(circuit, expected_text_diagram)
expected_text_diagram = """
0: ---a*a*aa(1, 1, 1)---------------------
|
1: ---a*a*aa------------------------------
|
2: ---a*a*aa------------a*a*aa(1, 1, 1)---
| |
3: ---a*a*aa------------a*a*aa------------
|
4: ---------------------a*a*aa------------
|
5: ---------------------a*a*aa------------
""".strip()
cirq.testing.assert_has_diagram(circuit,
expected_text_diagram,
use_unicode_characters=False)
@pytest.mark.parametrize(
'weights,exponent',
[(np.random.uniform(-5, 5, 3) + 1j * np.random.uniform(-5, 5, 3),
np.random.uniform(-5, 5)) for _ in range(5)])
def test_quartic_fermionic_simulation_apply_unitary(weights, exponent):
gate = ofc.QuarticFermionicSimulationGate(weights, exponent=exponent)
cirq.testing.assert_has_consistent_apply_unitary(gate, atol=5e-6)
