wip rust debugging, iswap+zz experiment

.pre-commit-config.yaml

@@ -13,7 +13,7 @@ repos:
       - id: conventional-commits
 
   - repo: https://github.com/charliermarsh/ruff-pre-commit
-    rev: "v0.1.3"
+    rev: "v0.1.5"
     hooks:
       - id: ruff
         args: [--fix, --exit-non-zero-on-fix]
@@ -43,7 +43,7 @@ repos:
         name: isort (python)
 
   - repo: https://github.com/psf/black
-    rev: 23.10.1
+    rev: 23.11.0
     hooks:
       - id: black-jupyter
         language_version: python3.9

Makefile

@@ -7,7 +7,7 @@ init:
 	rm -rf .venv
 	$(PYTHON_VERSION) -m venv .venv
 	@$(PIP) install --upgrade pip
-	$(PIP) install -e .[dev] --quiet
+	$(PIP) install -e .[dev]
 	@$(PRE_COMMIT) install
 	@$(PRE_COMMIT) autoupdate
 	sudo chmod +x .git/hooks/pre-commit

rustlib/src/lib.rs

@@ -1,81 +1,45 @@
-use numpy::{PyArray2, ToPyArray};
 use pyo3::prelude::*;
 use ndarray::Array2;
+use ndarray::ArrayBase;
+use ndarray::OwnedRepr;
+use ndarray::Dim;
 use num_complex::Complex64;
 use rayon::prelude::*;
+use numpy::{PyArray2, ToPyArray};
 
-fn dag(matrix: &Array2<Complex64>) -> Array2<Complex64> {
-    matrix.mapv(|elem| elem.conj()).reversed_axes()
+fn dag(op: &Array2<Complex64>) -> Array2<Complex64> {
+    op.t().mapv(|elem| elem.conj())
 }
 
-
 #[pyfunction]
-fn apply_error_channel_rust(
+fn apply_operators_in_place(
     py: Python,
     state: &PyArray2<Complex64>,
     num_steps: usize,
-    kraus_operators: Vec<&PyArray2<Complex64>>,
+    operators: Vec<&PyArray2<Complex64>>
 ) -> PyResult<PyObject> {
     let mut new_state = state.to_owned_array();
-    let kraus_ops: Vec<Array2<Complex64>> = kraus_operators
+    let ops: Vec<Array2<Complex64>> = operators
         .iter()
         .map(|&op| op.to_owned_array())
         .collect();
 
-    for _ in 0..num_steps {
-        let temp_state: Array2<Complex64> = kraus_ops.par_iter()
-            .map(|e| {
-                let e_dag = dag(e);
-                e.dot(&new_state).dot(&e_dag)
-            })
-            .reduce_with(|a, b| a + b)
-            .unwrap_or_else(|| Array2::<Complex64>::zeros(new_state.dim()));
-
-        new_state = temp_state;
-    }
-
-    Ok(new_state.to_pyarray(py).to_owned().into())
-}
-
-
-#[pyfunction]
-fn apply_error_channel_with_unitary(
-    py: Python,
-    state: &PyArray2<Complex64>,
-    num_steps: usize,
-    kraus_operators: Vec<&PyArray2<Complex64>>,
-    fractional_unitary: Option<&PyArray2<Complex64>>,
-) -> PyResult<PyObject> {
-    let mut new_state = state.to_owned_array();
-    let kraus_ops: Vec<Array2<Complex64>> = kraus_operators
-        .iter()
-        .map(|&op| op.to_owned_array())
-        .collect();
+    let dags: Vec<Array2<Complex64>> = ops.iter().map(|op| dag(op)).collect();
 
     for _ in 0..num_steps {
-        // Apply the fractional unitary if it's provided
-        if let Some(frac_unitary) = fractional_unitary {
-            let frac_unitary_array = frac_unitary.to_owned_array();
-            new_state = frac_unitary_array.dot(&new_state).dot(&frac_unitary_array.t().mapv(|elem| elem.conj()));
-        }
+        let temp_states: Vec<Array2<Complex64>> = ops.par_iter().zip(dags.par_iter())
+            .map(|(op, op_dag)| op.dot(&new_state).dot(op_dag))
+            .collect();
 
-        // Apply the error channel
-        let mut temp_state = Array2::<Complex64>::zeros(new_state.dim());
-        for e in &kraus_ops {
-            let e_dag = dag(e);
-            temp_state += &e.dot(&new_state).dot(&e_dag);
-        }
-        new_state = temp_state;
+        let sum_state: Array2<Complex64> = temp_states.into_par_iter().reduce_with(|a, b| a + b).unwrap();
+        new_state = sum_state;
     }
 
     Ok(new_state.to_pyarray(py).to_owned().into())
 }
 
-
-
 #[pymodule]
 fn _lib(py: Python, m: &PyModule) -> PyResult<()> {
-    m.add_function(wrap_pyfunction!(apply_error_channel_rust, m)?)?;
-    m.add_function(wrap_pyfunction!(apply_error_channel_with_unitary, m)?)?;
+    m.add_function(wrap_pyfunction!(apply_operators_in_place, m)?)?;
     Ok(())
 }

src/notebooks/noisy_simulation/04_iswap_crosskerr.ipynb

@@ -0,0 +1,172 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": 334,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "from quantum_logical._lib import apply_operators_in_place as rust_apply\n",
+        "from quantum_logical.error_channel import AmplitudeDamping, PhaseDamping\n",
+        "import numpy as np\n",
+        "import matplotlib.pyplot as plt\n",
+        "from qutip.metrics import fidelity\n",
+        "from qutip.operators import sigmax, sigmay, sigmaz\n",
+        "from qutip import basis\n",
+        "from qutip import expect\n",
+        "from quantum_logical.error_channel import ErrorChannel\n",
+        "from qutip import Qobj, tensor\n",
+        "from scipy.linalg import fractional_matrix_power"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 335,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "from quantum_logical.interaction import ConversionGainInteraction\n",
+        "\n",
+        "gc, gg = np.pi / 2, 0\n",
+        "trotter_dt = 0.1\n",
+        "duration = 1.0"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 336,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "one_rho = basis(2, 1) * basis(2, 1).dag()  # |1><1|\n",
+        "initial_state = tensor([one_rho] * 2)\n",
+        "\n",
+        "\n",
+        "# |01><01|\n",
+        "initial_state = (\n",
+        "    tensor([basis(2, 0), basis(2, 1)]) * tensor([basis(2, 0), basis(2, 1)]).dag()\n",
+        ")\n",
+        "\n",
+        "# bell state |00> + |11>\n",
+        "initial_state = tensor([basis(2, 0), basis(2, 0)]) + tensor([basis(2, 1), basis(2, 1)])\n",
+        "initial_state = initial_state * initial_state.dag()\n",
+        "\n",
+        "# |00> + |01>\n",
+        "initial_state = tensor([basis(2, 0), basis(2, 0)]) + tensor([basis(2, 0), basis(2, 1)])\n",
+        "initial_state = initial_state * initial_state.dag()\n",
+        "\n",
+        "\n",
+        "initial_state /= initial_state.norm()\n",
+        "initial_state = initial_state.full()\n",
+        "initial_state = Qobj(initial_state)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 337,
+      "metadata": {},
+      "outputs": [
+        {
+          "data": {
+            "text/latex": [
+              "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False $ \\\\ \\left(\\begin{matrix}1.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & -1.0j & 0.0\\\\0.0 & -1.0j & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 1.0\\\\\\end{matrix}\\right)$"
+            ],
+            "text/plain": [
+              "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = False\n",
+              "Qobj data =\n",
+              "[[1.+0.j 0.+0.j 0.+0.j 0.+0.j]\n",
+              " [0.+0.j 0.+0.j 0.-1.j 0.+0.j]\n",
+              " [0.+0.j 0.-1.j 0.+0.j 0.+0.j]\n",
+              " [0.+0.j 0.+0.j 0.+0.j 1.+0.j]]"
+            ]
+          },
+          "execution_count": 337,
+          "metadata": {},
+          "output_type": "execute_result"
+        }
+      ],
+      "source": [
+        "H = ConversionGainInteraction(gc, gg)\n",
+        "iswap = H.construct_U(t=duration)\n",
+        "iswap"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 348,
+      "metadata": {},
+      "outputs": [
+        {
+          "data": {
+            "text/latex": [
+              "Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False $ \\\\ \\left(\\begin{matrix}(2.441\\times10^{+06}+9.604\\times10^{-12}j) & (2.430\\times10^{+06}+9.475\\times10^{+04}j) & (2.430\\times10^{+06}+9.475\\times10^{+04}j) & (2.441\\times10^{+06}+9.604\\times10^{-12}j)\\\\(2.430\\times10^{+06}-9.475\\times10^{+04}j) & (2.441\\times10^{+06}+1.029\\times10^{-11}j) & (2.441\\times10^{+06}+2.851\\times10^{-10}j) & (2.430\\times10^{+06}-9.475\\times10^{+04}j)\\\\(2.430\\times10^{+06}-9.475\\times10^{+04}j) & (2.441\\times10^{+06}-2.517\\times10^{-10}j) & (2.441\\times10^{+06}+1.201\\times10^{-11}j) & (2.430\\times10^{+06}-9.475\\times10^{+04}j)\\\\(2.441\\times10^{+06}+9.604\\times10^{-12}j) & (2.430\\times10^{+06}+9.475\\times10^{+04}j) & (2.430\\times10^{+06}+9.475\\times10^{+04}j) & (2.441\\times10^{+06}+9.604\\times10^{-12}j)\\\\\\end{matrix}\\right)$"
+            ],
+            "text/plain": [
+              "Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False\n",
+              "Qobj data =\n",
+              "[[2441406.5       +9.60357713e-12j 2430220.66960511+9.47506831e+04j\n",
+              "  2430220.66960511+9.47511831e+04j 2441406.        +9.60357713e-12j]\n",
+              " [2430220.66960511-9.47506831e+04j 2441406.        +1.02876332e-11j\n",
+              "  2441406.        +2.85089583e-10j 2430220.66960511-9.47506831e+04j]\n",
+              " [2430220.66960511-9.47511831e+04j 2441406.        -2.51666123e-10j\n",
+              "  2441406.5       +1.20071995e-11j 2430220.66960511-9.47506831e+04j]\n",
+              " [2441406.        +9.60357713e-12j 2430220.66960511+9.47506831e+04j\n",
+              "  2430220.66960511+9.47506831e+04j 2441406.        +9.60357713e-12j]]"
+            ]
+          },
+          "execution_count": 348,
+          "metadata": {},
+          "output_type": "execute_result"
+        }
+      ],
+      "source": [
+        "class MultiQubitErrorChannel(ErrorChannel):\n",
+        "    \"\"\"A class for error channels that can be applied to multi-qubit systems.\"\"\"\n",
+        "\n",
+        "    def __init__(self, chi, trotter_dt):\n",
+        "        \"\"\"\n",
+        "        Initializes with given arrays of T1 and T2 times for each qubit, and a trotter step size.\n",
+        "        \"\"\"\n",
+        "        self.N = 2\n",
+        "        self.chi = chi\n",
+        "        self.num_qubits = 2\n",
+        "        dims = 2**self.num_qubits\n",
+        "        self.trotter_dt = trotter_dt\n",
+        "        super().__init__(trotter_dt, dims)\n",
+        "\n",
+        "    def _init_kraus_operators(self):\n",
+        "        \"\"\"\n",
+        "        Creates the multi-qubit Kraus operators for the error channels.\n",
+        "        \"\"\"\n",
+        "        zz = tensor(sigmaz(), sigmaz())\n",
+        "        return (-1j * self.chi * zz * np.pi / 2 * self.trotter_dt).expm()\n",
+        "\n",
+        "\n",
+        "noise_channel = MultiQubitErrorChannel(0, trotter_dt)\n",
+        "new_state = noise_channel.apply_error_channel(initial_state, 1.0, unitary=iswap.full())\n",
+        "new_state"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": ".venv",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.9.16"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 2
+}