Qtcorgi simulators

doing_a_test.py

@@ -0,0 +1,117 @@
+# import jax
+import pennylane as qml
+import numpy as np
+# import scipy
+#
+# n_wires = 2
+# num_qscripts = 2
+# qscripts = []
+# for i in range(num_qscripts):
+#     unitary = scipy.stats.unitary_group(dim=3**n_wires, seed=(42 + i)).rvs()
+#     op = qml.QutritUnitary(unitary, wires=range(n_wires))
+#     qs = qml.tape.QuantumScript([op], [qml.expval(qml.GellMann(0, 3))])
+#     qscripts.append(qs)
+#
+# dev = qml.devices.DefaultQutritMixed()
+# program, execution_config = dev.preprocess()
+# new_batch, post_processing_fn = program(qscripts)
+# results = dev.execute(new_batch, execution_config=execution_config)
+# print(post_processing_fn(results))
+#
+# @jax.jit
+# def f(x):
+#     qs = qml.tape.QuantumScript([qml.TRX(x, 0)], [qml.expval(qml.GellMann(0, 3))])
+#     program, execution_config = dev.preprocess()
+#     new_batch, post_processing_fn = program([qs])
+#     results = dev.execute(new_batch, execution_config=execution_config)
+#     return post_processing_fn(results)[0]
+#
+# jax.config.update("jax_enable_x64", True)
+# print(f(jax.numpy.array(1.2)))
+# print(jax.grad(f)(jax.numpy.array(1.2)))
+
+# 1/2, 1/3, 1/6
+
+# 1/2, 1/3, 1/6
+#
+# 1/2+1/3-2/6
+#
+# 2/6-3/6
+#
+#
+# gellMann_8_coeffs = np.array([1/np.sqrt(3), -1, 1/np.sqrt(3), 1, 1/np.sqrt(3), -1]) / 2
+# gellMann_8_obs = [qml.GellMann(0, i) for i in [1, 2, 4, 5, 6, 7]]
+# H1 = qml.Hamiltonian(gellMann_8_coeffs, gellMann_8_obs).matrix()
+#
+# G8 = qml.GellMann.compute_matrix(8)
+# Had = qml.THadamard.compute_matrix()
+# aHad = np.conj(Had).T
+#
+# H2 = np.round(aHad@G8@Had, 5)
+#
+#
+# # print(np.allclose(H1, H2))
+# # print(H1)
+# # print(H2)
+# # print(np.round(H2*np.sqrt(3), 5))
+# obs = aHad@G8@Had
+#
+# diag_gates = qml.THermitian(obs, 0).diagonalizing_gates()#[0].matrix()
+# print(len(diag_gates))
+# diag_gates = diag_gates[0].matrix()
+#
+# print(np.round((np.conj(diag_gates).T)@obs@diag_gates, 5))
+
+
+def setup_state(nr_wires):
+    """Sets up a basic state used for testing."""
+    setup_unitary = np.array(
+        [
+            [1 / np.sqrt(2), 1 / np.sqrt(3), 1 / np.sqrt(6)],
+            [np.sqrt(2 / 29), np.sqrt(3 / 29), -2 * np.sqrt(6 / 29)],
+            [-5 / np.sqrt(58), 7 / np.sqrt(87), 1 / np.sqrt(174)],
+        ]
+    ).T
+    qml.QutritUnitary(setup_unitary, wires=0)
+    qml.QutritUnitary(setup_unitary, wires=1)
+    if nr_wires == 3:
+        qml.TAdd(wires=(0, 2))
+
+
+dev = qml.device(
+        "default.qutrit.mixed",
+        wires=2,
+        damping_measurement_gammas=(0.2, 0.1, 0.4),
+        trit_flip_measurement_probs=(0.1, 0.2, 0.5),
+    )
+# Create matricies for the observables with diagonalizing matrix :math:`THadamard^\dag`
+inv_sqrt_3_i = 1j / np.sqrt(3)
+non_commuting_obs_one = np.array(
+    [
+        [0, -1 + inv_sqrt_3_i, -1 - inv_sqrt_3_i],
+        [-1 - inv_sqrt_3_i, 0, -1 + inv_sqrt_3_i],
+        [-1 + inv_sqrt_3_i, -1 - inv_sqrt_3_i, 0],
+    ]
+)
+non_commuting_obs_one /= 2
+
+@qml.qnode(dev)
+def circuit():
+    setup_state(2)
+
+    qml.THadamard(wires=0)
+    qml.THadamard(wires=1, subspace=(0, 1))
+
+    return qml.expval(qml.THermitian(non_commuting_obs_one, 0))
+
+
+
+
+
+print(my_test())
+
+
+
+
+
+

my_testing.py

@@ -0,0 +1,97 @@
+import numpy as np
+import pennylane as qml
+#
+# U = np.array([[1/np.sqrt(2), 1/np.sqrt(3), 1/np.sqrt(6)], [np.sqrt(2/29), np.sqrt(3/29), -2 * np.sqrt(6/29)], [-5/np.sqrt(58), 7/np.sqrt(87), 1/np.sqrt(174)]])
+#
+#
+# print(np.linalg.norm(U[1]))
+
+# inv_sqrt_3 = 1 / np.sqrt(3)
+# gellMann_8_coeffs = np.array([inv_sqrt_3, -1, inv_sqrt_3, 1, inv_sqrt_3, -1]) / 2
+# gellMann_8_obs = [qml.GellMann(0, i) for i in [1, 2, 4, 5, 6, 7]]
+
+# obs = qml.expval(qml.Hamiltonian(gellMann_8_coeffs, gellMann_8_obs))
+#
+# obs.matrix()
+#
+# ham = qml.Hamiltonian()
+
+# from scipy.stats import unitary_group
+# X = qml.PauliX.compute_matrix()
+# Y = qml.PauliX.compute_matrix()
+# Z = qml.PauliX.compute_matrix()
+# I = np.eye(2)
+
+
+
+#gellMann_8_obs = [qml.Hermitian(sum([np.random.rand()*m for m in [X, Y, Z, I]]), wires=0)]
+# for obs in gellMann_8_obs:
+#     diagm = obs.diagonalizing_gates()[0].matrix()
+#     diagm_adj = np.conj(diagm).T
+#     obsm = obs.matrix()
+#
+#     print(np.round(diagm@obsm@diagm_adj, 5), "\n")
+#     print(np.round(diagm_adj @ diagm, 5), "\n===============================================================\n")
+had = qml.THadamard.compute_matrix()
+ahad = np.conj(had).T
+G8 = qml.GellMann.compute_matrix(8)
+G3 = qml.GellMann.compute_matrix(3)
+
+
+# print(np.round(had@G8@ahad, 5))
+#
+# print()
+# print(np.round(had@G3@ahad, 5))
+
+inv_sqrt_3 = 1 / np.sqrt(3)
+inv_sqrt_3_i = inv_sqrt_3 * 1j
+#
+gellMann_3_equivalent = (
+            np.array(
+                [[0, 1+inv_sqrt_3_i, 1-inv_sqrt_3_i],
+                 [1-inv_sqrt_3_i, 0, 1 + inv_sqrt_3_i],
+                 [1+inv_sqrt_3_i, 1-inv_sqrt_3_i, 0]]
+            )
+            / 2
+        )
+gellMann_8_equivalent = (
+                np.array(
+                    [[0, (inv_sqrt_3 - 1j), (inv_sqrt_3 + 1j)],
+                     [inv_sqrt_3 + 1j, 0, inv_sqrt_3 - 1j],
+                     [inv_sqrt_3 - 1j, inv_sqrt_3 + 1j, 0]]
+                )
+                / 2
+        )
+
+dg = qml.THermitian(gellMann_8_equivalent, 0).diagonalizing_gates()[0].matrix()
+# dga = np.conj(dg).T
+# print(np.round(dg@gellMann_8_equivalent@dga, 5))
+#
+# print(np.abs((dg@had@(np.array([1/2,1/3,1/6])**(1/2))))**2)
+#print(qml.GellMann(0, 1).diagonalizing_gates()[0].matrix())
+
+# print(np.round(dg@had, 4))
+obs = np.diag([1, 2, 3])
+print(np.round(had@obs@ahad, 4))
+
+obs = np.diag([-2, -1, 1])
+
+non_commuting_obs_two = np.array(
+            [
+                [-2/3, -2/3 + inv_sqrt_3_i, -2/3 - inv_sqrt_3_i],
+                [-2/3 - inv_sqrt_3_i, -2/3, -2/3 + inv_sqrt_3_i],
+                [-2/3 + inv_sqrt_3_i, -2/3 - inv_sqrt_3_i, -2/3],
+            ]
+        )
+
+print(np.round(had@obs@ahad, 4))
+
+print(np.allclose(had@obs@ahad, non_commuting_obs_two))
+# print(np.allclose(had@G8@ahad, gellMann_8_equivalent))
+# #print(had@G8@ahad)
+# print(np.allclose(ahad, dg))
+import jax
+print(jax.numpy.array([jax.numpy.nan, 1., 2.]))
+
+
+

pennylane/devices/default_mixed.py

@@ -48,6 +48,7 @@
 from pennylane.wires import Wires
 
 from .._version import __version__
+from .qtcorgi_helper.qtcorgi_simulator import get_qubit_final_state_from_initial
 
 logger = logging.getLogger(__name__)
 logger.addHandler(logging.NullHandler())
@@ -781,9 +782,17 @@ def apply(self, operations, rotations=None, **kwargs):
                     f"on a {self.short_name} device."
                 )
 
-        for operation in operations:
-            self._apply_operation(operation)
-
+        prep = False
+        if len(operations) > 0:
+            if (
+                isinstance(operations[0], StatePrep)
+                or isinstance(operations[0], BasisState)
+                or isinstance(operations[0], QubitDensityMatrix)
+            ):
+                prep = True
+                self._apply_operation(operation)
+
+        self._state = get_qubit_final_state_from_initial(operations[prep:], self._state)
         # store the pre-rotated state
         self._pre_rotated_state = self._state