[chore]: unwrap project

9 files changed, 0 insertions, 3144 deletions

diff --git a/libpsi-core/src/core/circuit.rs b/libpsi-core/src/core/circuit.rs
deleted file mode 100644
index d0e4450..0000000
--- a/libpsi-core/src/core/circuit.rs
+++ /dev/null
@@ -1,481 +0,0 @@
-use super::{CustomGate, QuantumState, Runtime, RuntimeConfig};
-use crate::{format_amplitude, format_probability, Vector};
-use core::fmt;
-use std::sync::Arc;
-
-#[derive(Clone)]
-pub enum GateOp {
-    H(usize),
-    X(usize),
-    Y(usize),
-    Z(usize),
-    S(usize),
-    T(usize),
-    Sdg(usize),
-    Tdg(usize),
-    Sx(usize),
-    Sxdg(usize),
-    Rx(usize, f64),
-    Ry(usize, f64),
-    Rz(usize, f64),
-    P(usize, f64),
-    U1(usize, f64),
-    U2(usize, f64, f64),
-    U3(usize, f64, f64, f64),
-    CNOT(usize, usize),
-    CZ(usize, usize),
-    SWAP(usize, usize),
-    CRx(usize, usize, f64),
-    CRy(usize, usize, f64),
-    CRz(usize, usize, f64),
-    CP(usize, usize, f64),
-    CCNOT(usize, usize, usize),
-    CSWAP(usize, usize, usize),
-    Measure(usize, usize),
-    Custom(Arc<CustomGate>, Vec<usize>),
-}
-
-impl GateOp {
-    pub fn name(&self) -> &str {
-        match self {
-            GateOp::H(_) => "H",
-            GateOp::X(_) => "X",
-            GateOp::Y(_) => "Y",
-            GateOp::Z(_) => "Z",
-            GateOp::S(_) => "S",
-            GateOp::T(_) => "T",
-            GateOp::Sdg(_) => "S†",
-            GateOp::Tdg(_) => "T†",
-            GateOp::Sx(_) => "√X",
-            GateOp::Sxdg(_) => "√X†",
-            GateOp::Rx(_, _) => "Rx",
-            GateOp::Ry(_, _) => "Ry",
-            GateOp::Rz(_, _) => "Rz",
-            GateOp::P(_, _) => "P",
-            GateOp::U1(_, _) => "U1",
-            GateOp::U2(_, _, _) => "U2",
-            GateOp::U3(_, _, _, _) => "U3",
-            GateOp::CRx(_, _, _) => "CRx",
-            GateOp::CRy(_, _, _) => "CRy",
-            GateOp::CRz(_, _, _) => "CRz",
-            GateOp::CP(_, _, _) => "CP",
-            GateOp::CNOT(_, _) => "CNOT",
-            GateOp::CZ(_, _) => "CZ",
-            GateOp::SWAP(_, _) => "SWAP",
-            GateOp::CCNOT(_, _, _) => "CCNOT",
-            GateOp::CSWAP(_, _, _) => "CSWAP",
-            GateOp::Measure(_, _) => "M",
-            GateOp::Custom(gate, _) => &gate.name,
-        }
-    }
-
-    pub fn quantum_targets(&self) -> Vec<usize> {
-        match self {
-            GateOp::H(t)
-            | GateOp::X(t)
-            | GateOp::Y(t)
-            | GateOp::Z(t)
-            | GateOp::S(t)
-            | GateOp::T(t)
-            | GateOp::Sdg(t)
-            | GateOp::Tdg(t)
-            | GateOp::Sx(t)
-            | GateOp::Sxdg(t)
-            | GateOp::Rx(t, _)
-            | GateOp::Ry(t, _)
-            | GateOp::Rz(t, _)
-            | GateOp::P(t, _)
-            | GateOp::U1(t, _)
-            | GateOp::U2(t, _, _)
-            | GateOp::U3(t, _, _, _) => vec![*t],
-            GateOp::CNOT(c, t)
-            | GateOp::CZ(c, t)
-            | GateOp::SWAP(c, t)
-            | GateOp::CRx(c, t, _)
-            | GateOp::CRy(c, t, _)
-            | GateOp::CRz(c, t, _)
-            | GateOp::CP(c, t, _) => vec![*c, *t],
-            GateOp::CCNOT(c1, c2, t) | GateOp::CSWAP(c1, c2, t) => vec![*c1, *c2, *t],
-            GateOp::Measure(q, _) => vec![*q],
-            GateOp::Custom(_, targets) => targets.clone(),
-        }
-    }
-
-    pub fn classical_targets(&self) -> Vec<usize> {
-        match self {
-            GateOp::Measure(_, c) => vec![*c],
-            _ => vec![],
-        }
-    }
-
-    pub fn is_measurement(&self) -> bool {
-        matches!(self, GateOp::Measure(_, _))
-    }
-
-    pub fn is_custom(&self) -> bool {
-        matches!(self, GateOp::Custom(_, _))
-    }
-
-    pub fn is_non_clifford(&self) -> bool {
-        matches!(
-            self,
-            GateOp::T(_)
-                | GateOp::Tdg(_)
-                | GateOp::Sx(_)
-                | GateOp::Sxdg(_)
-                | GateOp::Rx(_, _)
-                | GateOp::Ry(_, _)
-                | GateOp::Rz(_, _)
-                | GateOp::P(_, _)
-                | GateOp::U1(_, _)
-                | GateOp::U2(_, _, _)
-                | GateOp::U3(_, _, _, _)
-                | GateOp::CRx(_, _, _)
-                | GateOp::CRy(_, _, _)
-                | GateOp::CRz(_, _, _)
-                | GateOp::CP(_, _, _)
-        )
-    }
-}
-
-pub struct QuantumCircuit {
-    num_qubits: usize,
-    num_classical: usize,
-    operations: Vec<GateOp>,
-    computed_state: Option<QuantumState>,
-}
-
-impl QuantumCircuit {
-    pub fn new(num_qubits: usize) -> QuantumCircuit {
-        QuantumCircuit {
-            num_qubits,
-            num_classical: 0,
-            operations: Vec::new(),
-            computed_state: None,
-        }
-    }
-
-    pub fn with_classical(num_qubits: usize, num_classical: usize) -> QuantumCircuit {
-        QuantumCircuit {
-            num_qubits,
-            num_classical,
-            operations: Vec::new(),
-            computed_state: None,
-        }
-    }
-
-    pub fn num_qubits(&self) -> usize {
-        self.num_qubits
-    }
-
-    pub fn num_classical(&self) -> usize {
-        self.num_classical
-    }
-
-    pub fn operations(&self) -> &[GateOp] {
-        &self.operations
-    }
-
-    pub fn is_computed(&self) -> bool {
-        self.computed_state.is_some()
-    }
-
-    pub fn compute(&mut self) -> &QuantumState {
-        self.compute_with(Runtime::default())
-    }
-
-    pub fn compute_with(&mut self, runtime: Runtime) -> &QuantumState {
-        if self.computed_state.is_some() {
-            return self.computed_state.as_ref().unwrap();
-        }
-
-        self.computed_state = Some(runtime.compute(self.num_qubits, &self.operations));
-        self.computed_state.as_ref().unwrap()
-    }
-
-    pub fn compute_with_config(&mut self, config: RuntimeConfig) -> &QuantumState {
-        if self.computed_state.is_some() {
-            return self.computed_state.as_ref().unwrap();
-        }
-
-        self.computed_state = Some(config.compute(self.num_qubits, &self.operations));
-        self.computed_state.as_ref().unwrap()
-    }
-
-    pub fn state(&mut self) -> &QuantumState {
-        self.compute()
-    }
-
-    pub fn state_with(&mut self, runtime: Runtime) -> &QuantumState {
-        self.compute_with(runtime)
-    }
-
-    pub fn state_with_config(&mut self, config: RuntimeConfig) -> &QuantumState {
-        self.compute_with_config(config)
-    }
-
-    pub fn h(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::H(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn x(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::X(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn y(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Y(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn z(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Z(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn s(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::S(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn t(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::T(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn sdg(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Sdg(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn tdg(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Tdg(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn sx(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Sx(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn sxdg(&mut self, target: usize) -> &mut Self {
-        self.operations.push(GateOp::Sxdg(target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn rx(&mut self, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::Rx(target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn ry(&mut self, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::Ry(target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn rz(&mut self, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::Rz(target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn p(&mut self, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::P(target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn u1(&mut self, target: usize, lambda: f64) -> &mut Self {
-        self.operations.push(GateOp::U1(target, lambda));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn u2(&mut self, target: usize, phi: f64, lambda: f64) -> &mut Self {
-        self.operations.push(GateOp::U2(target, phi, lambda));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn u3(&mut self, target: usize, theta: f64, phi: f64, lambda: f64) -> &mut Self {
-        self.operations.push(GateOp::U3(target, theta, phi, lambda));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn crx(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::CRx(control, target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn cry(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::CRy(control, target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn crz(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::CRz(control, target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn cp(&mut self, control: usize, target: usize, theta: f64) -> &mut Self {
-        self.operations.push(GateOp::CP(control, target, theta));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn cnot(&mut self, control: usize, target: usize) -> &mut Self {
-        self.operations.push(GateOp::CNOT(control, target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn cx(&mut self, control: usize, target: usize) -> &mut Self {
-        self.cnot(control, target)
-    }
-
-    pub fn cz(&mut self, control: usize, target: usize) -> &mut Self {
-        self.operations.push(GateOp::CZ(control, target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn swap(&mut self, qubit1: usize, qubit2: usize) -> &mut Self {
-        self.operations.push(GateOp::SWAP(qubit1, qubit2));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn ccnot(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self {
-        self.operations
-            .push(GateOp::CCNOT(control1, control2, target));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn toffoli(&mut self, control1: usize, control2: usize, target: usize) -> &mut Self {
-        self.ccnot(control1, control2, target)
-    }
-
-    pub fn cswap(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self {
-        self.operations
-            .push(GateOp::CSWAP(control, target1, target2));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn fredkin(&mut self, control: usize, target1: usize, target2: usize) -> &mut Self {
-        self.cswap(control, target1, target2)
-    }
-
-    pub fn measure(&mut self, qubit: usize, classical: usize) -> &mut Self {
-        if classical >= self.num_classical {
-            self.num_classical = classical + 1;
-        }
-        self.operations.push(GateOp::Measure(qubit, classical));
-        self
-    }
-
-    pub fn measure_all(&mut self) -> &mut Self {
-        for i in 0..self.num_qubits {
-            self.measure(i, i);
-        }
-        self
-    }
-
-    pub fn custom(&mut self, gate: &Arc<CustomGate>, targets: &[usize]) -> &mut Self {
-        self.operations
-            .push(GateOp::Custom(Arc::clone(gate), targets.to_vec()));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn apply_custom(&mut self, gate: CustomGate, targets: &[usize]) -> &mut Self {
-        self.operations
-            .push(GateOp::Custom(Arc::new(gate), targets.to_vec()));
-        self.computed_state = None;
-        self
-    }
-
-    pub fn reset(&mut self) -> &mut Self {
-        self.operations.clear();
-        self.computed_state = None;
-        self
-    }
-
-    pub fn probability(&mut self, state_index: usize) -> f64 {
-        self.compute();
-        let state = self.computed_state.as_ref().unwrap();
-        let amp = state.get(state_index);
-        amp.norm2()
-    }
-
-    pub fn probabilities(&mut self) -> Vec<f64> {
-        self.compute();
-        let n = 1 << self.num_qubits;
-        let state = self.computed_state.as_ref().unwrap();
-        (0..n).map(|i| state.get(i).norm2()).collect()
-    }
-
-    pub fn print_probabilities(&mut self) {
-        let probs = self.probabilities();
-        let n = self.num_qubits;
-        println!("Probabilities:");
-        for (i, p) in probs.iter().enumerate() {
-            if *p > 1e-10 {
-                let basis: String = format!("{:0width$b}", i, width = n);
-                println!("  |{}⟩: {}", basis, format_probability(*p));
-            }
-        }
-    }
-}
-
-impl fmt::Display for QuantumCircuit {
-    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
-        writeln!(
-            f,
-            "QuantumCircuit ({} qubits, {} classical)",
-            self.num_qubits, self.num_classical
-        )?;
-        writeln!(f, "Operations:")?;
-        for (i, op) in self.operations.iter().enumerate() {
-            match op {
-                GateOp::Measure(q, c) => writeln!(f, "  {}: {} q{} → c{}", i, op.name(), q, c)?,
-                GateOp::Custom(gate, targets) => {
-                    writeln!(f, "  {}: [{}] on {:?}", i, gate.name, targets)?
-                }
-                _ => writeln!(f, "  {}: {} on {:?}", i, op.name(), op.quantum_targets())?,
-            }
-        }
-        if let Some(state) = &self.computed_state {
-            writeln!(f, "State:")?;
-            let n = 1 << self.num_qubits;
-            for i in 0..n {
-                let amp = state.get(i);
-                if amp.real.abs() > 1e-10 || amp.imaginary.abs() > 1e-10 {
-                    let basis: String = format!("{:0width$b}", i, width = self.num_qubits);
-                    writeln!(f, "  |{}⟩: {}", basis, format_amplitude(&amp))?;
-                }
-            }
-        } else {
-            writeln!(f, "State: (not computed)")?;
-        }
-        Ok(())
-    }
-}
diff --git a/libpsi-core/src/core/classical_components.rs b/libpsi-core/src/core/classical_components.rs
deleted file mode 100644
index f6565af..0000000
--- a/libpsi-core/src/core/classical_components.rs
+++ /dev/null
@@ -1,63 +0,0 @@
-use core::ops;
-
-#[derive(Clone, Copy)]
-pub struct ClassicalBit<'a> {
-    state: bool,
-    name: &'a str,
-}
-
-#[derive(Clone)]
-pub struct ClassicalRegister<'a> {
-    bits: Vec<ClassicalBit<'a>>,
-    name: &'a str,
-}
-
-impl<'a> ClassicalBit<'a> {
-    pub fn new(name: &'a str, state: bool) -> ClassicalBit<'a> {
-        ClassicalBit { name, state }
-    }
-
-    pub fn get_name(&self) -> &'a str {
-        self.name
-    }
-
-    pub fn get_state(&self) -> bool {
-        self.state
-    }
-}
-
-impl<'a> ClassicalRegister<'a> {
-    pub fn new(name: &'a str, names: &'a [&'a str]) -> ClassicalRegister<'a> {
-        let mut bits: Vec<ClassicalBit<'a>> = Vec::new();
-        for i in 0..names.len() {
-            bits.push(ClassicalBit::new(names[i], false));
-        }
-        ClassicalRegister { name, bits }
-    }
-
-    pub fn set_bits(&mut self, bits: Vec<ClassicalBit<'a>>) {
-        self.bits = bits;
-    }
-
-    pub fn get_bits(&self) -> Vec<ClassicalBit<'a>> {
-        self.bits.clone()
-    }
-
-    pub fn get_name(&self) -> &'a str {
-        self.name
-    }
-}
-
-impl<'a> ops::Index<usize> for ClassicalRegister<'a> {
-    type Output = ClassicalBit<'a>;
-
-    fn index(&self, index: usize) -> &Self::Output {
-        &self.bits[index]
-    }
-}
-
-impl<'a> ops::IndexMut<usize> for ClassicalRegister<'a> {
-    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
-        &mut self.bits[index]
-    }
-}
diff --git a/libpsi-core/src/core/custom_gate.rs b/libpsi-core/src/core/custom_gate.rs
deleted file mode 100644
index 9e49ac9..0000000
--- a/libpsi-core/src/core/custom_gate.rs
+++ /dev/null
@@ -1,245 +0,0 @@
-use crate::{Complex, Matrix, QuantumGate};
-
-#[derive(Clone)]
-pub enum CustomGateDefinition {
-    Matrix(Matrix<Complex<f64>>),
-    Composite(Vec<(CompositeOp, Vec<usize>)>),
-}
-
-#[derive(Clone, Copy)]
-pub enum CompositeOp {
-    H,
-    X,
-    Y,
-    Z,
-    S,
-    T,
-    CNOT,
-    CZ,
-    SWAP,
-    CCNOT,
-    CSWAP,
-}
-
-#[derive(Clone)]
-pub struct CustomGate {
-    pub name: String,
-    pub num_qubits: usize,
-    pub definition: CustomGateDefinition,
-}
-
-impl CustomGate {
-    pub fn from_matrix(name: &str, matrix: Matrix<Complex<f64>>) -> Self {
-        let dim = matrix.rows;
-        let num_qubits = (dim as f64).log2() as usize;
-        assert_eq!(
-            1 << num_qubits,
-            dim,
-            "Matrix dimension must be a power of 2"
-        );
-        assert_eq!(matrix.rows, matrix.cols, "Matrix must be square");
-
-        CustomGate {
-            name: String::from(name),
-            num_qubits,
-            definition: CustomGateDefinition::Matrix(matrix),
-        }
-    }
-
-    pub fn from_composite(
-        name: &str,
-        num_qubits: usize,
-        ops: Vec<(CompositeOp, Vec<usize>)>,
-    ) -> Self {
-        CustomGate {
-            name: String::from(name),
-            num_qubits,
-            definition: CustomGateDefinition::Composite(ops),
-        }
-    }
-
-    pub fn to_quantum_gate(&self) -> QuantumGate<'static> {
-        match &self.definition {
-            CustomGateDefinition::Matrix(matrix) => {
-                let name: &'static str = Box::leak(self.name.clone().into_boxed_str());
-                QuantumGate {
-                    name,
-                    matrix: matrix.clone(),
-                    num_qubits: self.num_qubits,
-                }
-            }
-            CustomGateDefinition::Composite(ops) => {
-                let matrix = self.compute_composite_matrix(ops);
-                let name: &'static str = Box::leak(self.name.clone().into_boxed_str());
-                QuantumGate {
-                    name,
-                    matrix,
-                    num_qubits: self.num_qubits,
-                }
-            }
-        }
-    }
-
-    fn compute_composite_matrix(&self, ops: &[(CompositeOp, Vec<usize>)]) -> Matrix<Complex<f64>> {
-        use crate::gates::*;
-        use crate::Complex;
-
-        let dim = 1 << self.num_qubits;
-        let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]);
-        for i in 0..dim {
-            result.data[i * dim + i] = Complex::new(1.0, 0.0);
-        }
-
-        for (op, targets) in ops {
-            let gate: &QuantumGate = match op {
-                CompositeOp::H => &HADAMARD,
-                CompositeOp::X => &PAULI_X,
-                CompositeOp::Y => &PAULI_Y,
-                CompositeOp::Z => &PAULI_Z,
-                CompositeOp::S => &S_GATE,
-                CompositeOp::T => &T_GATE,
-                CompositeOp::CNOT => &CNOT,
-                CompositeOp::CZ => &CZ,
-                CompositeOp::SWAP => &SWAP,
-                CompositeOp::CCNOT => &TOFFOLI,
-                CompositeOp::CSWAP => &FREDKIN,
-            };
-
-            let full_gate = build_full_operator(&gate.matrix, targets, self.num_qubits);
-            result = matrix_multiply(&full_gate, &result);
-        }
-
-        result
-    }
-}
-
-fn build_full_operator(
-    gate_matrix: &Matrix<Complex<f64>>,
-    targets: &[usize],
-    total_qubits: usize,
-) -> Matrix<Complex<f64>> {
-    let dim = 1 << total_qubits;
-    let gate_dim = gate_matrix.rows;
-    let num_gate_qubits = targets.len();
-
-    let mut result = Matrix::new(dim, dim, vec![Complex::new(0.0, 0.0); dim * dim]);
-
-    for i in 0..dim {
-        for j in 0..dim {
-            let mut gate_i = 0usize;
-            let mut gate_j = 0usize;
-            let mut match_non_targets = true;
-
-            for q in 0..total_qubits {
-                let bit_i = (i >> (total_qubits - 1 - q)) & 1;
-                let bit_j = (j >> (total_qubits - 1 - q)) & 1;
-
-                if let Some(pos) = targets.iter().position(|&t| t == q) {
-                    gate_i |= bit_i << (num_gate_qubits - 1 - pos);
-                    gate_j |= bit_j << (num_gate_qubits - 1 - pos);
-                } else if bit_i != bit_j {
-                    match_non_targets = false;
-                    break;
-                }
-            }
-
-            if match_non_targets {
-                result.data[i * dim + j] = gate_matrix.data[gate_i * gate_dim + gate_j];
-            }
-        }
-    }
-
-    result
-}
-
-fn matrix_multiply(a: &Matrix<Complex<f64>>, b: &Matrix<Complex<f64>>) -> Matrix<Complex<f64>> {
-    let n = a.rows;
-    let mut result = Matrix::new(n, n, vec![Complex::new(0.0, 0.0); n * n]);
-
-    for i in 0..n {
-        for j in 0..n {
-            let mut sum = Complex::new(0.0, 0.0);
-            for k in 0..n {
-                sum = sum + a.data[i * n + k] * b.data[k * n + j];
-            }
-            result.data[i * n + j] = sum;
-        }
-    }
-
-    result
-}
-
-pub struct CustomGateBuilder {
-    name: String,
-    num_qubits: usize,
-    ops: Vec<(CompositeOp, Vec<usize>)>,
-}
-
-impl CustomGateBuilder {
-    pub fn new(name: &str, num_qubits: usize) -> Self {
-        CustomGateBuilder {
-            name: String::from(name),
-            num_qubits,
-            ops: Vec::new(),
-        }
-    }
-
-    pub fn h(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::H, vec![target]));
-        self
-    }
-
-    pub fn x(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::X, vec![target]));
-        self
-    }
-
-    pub fn y(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::Y, vec![target]));
-        self
-    }
-
-    pub fn z(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::Z, vec![target]));
-        self
-    }
-
-    pub fn s(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::S, vec![target]));
-        self
-    }
-
-    pub fn t(mut self, target: usize) -> Self {
-        self.ops.push((CompositeOp::T, vec![target]));
-        self
-    }
-
-    pub fn cnot(mut self, control: usize, target: usize) -> Self {
-        self.ops.push((CompositeOp::CNOT, vec![control, target]));
-        self
-    }
-
-    pub fn cz(mut self, control: usize, target: usize) -> Self {
-        self.ops.push((CompositeOp::CZ, vec![control, target]));
-        self
-    }
-
-    pub fn swap(mut self, a: usize, b: usize) -> Self {
-        self.ops.push((CompositeOp::SWAP, vec![a, b]));
-        self
-    }
-
-    pub fn ccnot(mut self, c1: usize, c2: usize, target: usize) -> Self {
-        self.ops.push((CompositeOp::CCNOT, vec![c1, c2, target]));
-        self
-    }
-
-    pub fn cswap(mut self, control: usize, t1: usize, t2: usize) -> Self {
-        self.ops.push((CompositeOp::CSWAP, vec![control, t1, t2]));
-        self
-    }
-
-    pub fn build(self) -> CustomGate {
-        CustomGate::from_composite(&self.name, self.num_qubits, self.ops)
-    }
-}
diff --git a/libpsi-core/src/core/gates.rs b/libpsi-core/src/core/gates.rs
deleted file mode 100644
index cee4ebc..0000000
--- a/libpsi-core/src/core/gates.rs
+++ /dev/null
@@ -1,253 +0,0 @@
-use crate::{complex, matrix, Complex, Matrix, QuantumGate};
-use std::f64::consts::FRAC_1_SQRT_2;
-
-pub fn rx_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let cos = (theta / 2.0).cos();
-    let sin = (theta / 2.0).sin();
-    matrix!(
-        [complex!(cos, 0.0), complex!(0.0, -sin)];
-        [complex!(0.0, -sin), complex!(cos, 0.0)]
-    )
-}
-
-pub fn ry_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let cos = (theta / 2.0).cos();
-    let sin = (theta / 2.0).sin();
-    matrix!(
-        [complex!(cos, 0.0), complex!(-sin, 0.0)];
-        [complex!(sin, 0.0), complex!(cos, 0.0)]
-    )
-}
-
-pub fn rz_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let half = theta / 2.0;
-    matrix!(
-        [complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(half.cos(), half.sin())]
-    )
-}
-
-pub fn p_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    matrix!(
-        [complex!(1.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())]
-    )
-}
-
-pub fn u1_matrix(lambda: f64) -> Matrix<Complex<f64>> {
-    p_matrix(lambda)
-}
-
-pub fn u2_matrix(phi: f64, lambda: f64) -> Matrix<Complex<f64>> {
-    let inv_sqrt2 = FRAC_1_SQRT_2;
-    matrix!(
-        [complex!(inv_sqrt2, 0.0), complex!(-inv_sqrt2 * lambda.cos(), -inv_sqrt2 * lambda.sin())];
-        [complex!(inv_sqrt2 * phi.cos(), inv_sqrt2 * phi.sin()), complex!((phi + lambda).cos() * inv_sqrt2, (phi + lambda).sin() * inv_sqrt2)]
-    )
-}
-
-pub fn u3_matrix(theta: f64, phi: f64, lambda: f64) -> Matrix<Complex<f64>> {
-    let cos = (theta / 2.0).cos();
-    let sin = (theta / 2.0).sin();
-    matrix!(
-        [complex!(cos, 0.0), complex!(-sin * lambda.cos(), -sin * lambda.sin())];
-        [complex!(sin * phi.cos(), sin * phi.sin()), complex!(cos * (phi + lambda).cos(), cos * (phi + lambda).sin())]
-    )
-}
-
-pub fn crx_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let cos = (theta / 2.0).cos();
-    let sin = (theta / 2.0).sin();
-    matrix!(
-        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(0.0, -sin)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, -sin), complex!(cos, 0.0)]
-    )
-}
-
-pub fn cry_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let cos = (theta / 2.0).cos();
-    let sin = (theta / 2.0).sin();
-    matrix!(
-        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(cos, 0.0), complex!(-sin, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(sin, 0.0), complex!(cos, 0.0)]
-    )
-}
-
-pub fn crz_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    let half = theta / 2.0;
-    matrix!(
-        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), -half.sin()), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(half.cos(), half.sin())]
-    )
-}
-
-pub fn cp_matrix(theta: f64) -> Matrix<Complex<f64>> {
-    matrix!(
-        [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
-        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(theta.cos(), theta.sin())]
-    )
-}
-
-#[rustfmt::skip]
-lazy_static::lazy_static! {
-    pub static ref HADAMARD: QuantumGate<'static> = QuantumGate {
-        name: "H",
-        matrix: matrix!([complex!(1.0, 0.0), complex!( 1.0, 0.0)];
-                        [complex!(1.0, 0.0), complex!(-1.0, 0.0)]) *
-                complex!(1.0/2.0_f64.sqrt(), 0.0),
-        num_qubits: 1,
-    };
-
-    pub static ref PAULI_X: QuantumGate<'static> = QuantumGate {
-        name: "X",
-        matrix: matrix!([complex!(0.0, 0.0), complex!(1.0, 0.0)];
-                        [complex!(1.0, 0.0), complex!(0.0, 0.0)]),
-        num_qubits: 1,
-    };
-
-    pub static ref PAULI_Y: QuantumGate<'static> = QuantumGate {
-        name: "Y", 
-        matrix: matrix!([complex!(0.0, 0.0), complex!(0.0, -1.0)];
-                        [complex!(0.0, 1.0), complex!(0.0,  0.0)]),
-        num_qubits: 1,
-    };
-
-    pub static ref PAULI_Z: QuantumGate<'static> = QuantumGate {
-        name: "Z", 
-        matrix: matrix!([complex!(1.0, 0.0), complex!( 0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(-1.0, 0.0)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref S_GATE: QuantumGate<'static> = QuantumGate {
-        name: "S",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 1.0)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref T_GATE: QuantumGate<'static> = QuantumGate {
-        name: "T",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, core::f64::consts::FRAC_1_SQRT_2)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref SDG_GATE: QuantumGate<'static> = QuantumGate {
-        name: "S†",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, -1.0)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref TDG_GATE: QuantumGate<'static> = QuantumGate {
-        name: "T†",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(core::f64::consts::FRAC_1_SQRT_2, -core::f64::consts::FRAC_1_SQRT_2)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref SX_GATE: QuantumGate<'static> = QuantumGate {
-        name: "√X",
-        matrix: matrix!([complex!(0.5, 0.5), complex!(0.5, -0.5)];
-                        [complex!(0.5, -0.5), complex!(0.5, 0.5)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref SXDG_GATE: QuantumGate<'static> = QuantumGate {
-        name: "√X†",
-        matrix: matrix!([complex!(0.5, -0.5), complex!(0.5, 0.5)];
-                        [complex!(0.5, 0.5), complex!(0.5, -0.5)]),
-        num_qubits: 1,
-    };
-    
-    pub static ref IDENTITY: QuantumGate<'static> = QuantumGate {
-        name: "I",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(1.0, 0.0)]),
-        num_qubits: 1,
-    };
-
-    pub static ref CNOT: QuantumGate<'static> = QuantumGate {
-        name: "CNOT", 
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]),
-        num_qubits: 2,
-    };
-    
-    pub static ref CZ: QuantumGate<'static> = QuantumGate {
-        name: "CZ",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!( 0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!( 0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(-1.0, 0.0)]),
-        num_qubits: 2,
-    };
-    
-    pub static ref SWAP: QuantumGate<'static> = QuantumGate {
-        name: "SWAP",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
-        num_qubits: 2,
-    };
-    
-    pub static ref ISWAP: QuantumGate<'static> = QuantumGate {
-        name: "iSWAP",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 1.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
-        num_qubits: 2,
-    };
-    
-    pub static ref SQRT_SWAP: QuantumGate<'static> = QuantumGate {
-        name: "√SWAP",
-        matrix: matrix!([complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.5, 0.5), complex!(0.5, -0.5), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.5, -0.5), complex!(0.5, 0.5), complex!(0.0, 0.0)];
-                        [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]),
-        num_qubits: 2,
-    };
-
-    pub static ref TOFFOLI: QuantumGate<'static> = QuantumGate {
-        name: "CCNOT",
-        matrix: matrix!(
-            [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)]
-        ),
-        num_qubits: 3,
-    };
-    
-    pub static ref FREDKIN: QuantumGate<'static> = QuantumGate {
-        name: "CSWAP",
-        matrix: matrix!(
-            [complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0)];
-            [complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(0.0, 0.0), complex!(1.0, 0.0)]
-        ),
-        num_qubits: 3,
-    };
-}
diff --git a/libpsi-core/src/core/kernel.rs b/libpsi-core/src/core/kernel.rs
deleted file mode 100644
index 7b66eea..0000000
--- a/libpsi-core/src/core/kernel.rs
+++ /dev/null
@@ -1,622 +0,0 @@
-use crate::maths::simd::{
-    apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel, SimdCapability,
-};
-use crate::{complex, Complex, Matrix};
-use rayon::prelude::*;
-use std::collections::HashSet;
-
-#[derive(Clone, Copy, Debug, PartialEq, Eq)]
-pub enum GateType {
-    Diagonal,
-    NonDiagonal,
-    Controlled,
-}
-
-#[derive(Clone)]
-pub struct Kernel {
-    pub matrix: Matrix<Complex<f64>>,
-    pub targets: Vec<usize>,
-    pub name: String,
-    pub gate_type: GateType,
-}
-
-impl Kernel {
-    pub fn new(name: &str, matrix: Matrix<Complex<f64>>, targets: Vec<usize>) -> Self {
-        let gate_type = Self::detect_gate_type(name, &matrix);
-        Self {
-            matrix,
-            targets,
-            name: name.to_string(),
-            gate_type,
-        }
-    }
-
-    fn detect_gate_type(name: &str, matrix: &Matrix<Complex<f64>>) -> GateType {
-        let diagonal_gates = [
-            "Z", "S", "T", "Sdg", "Tdg", "Rz", "P", "U1", "CZ", "CP", "CRz",
-        ];
-        if diagonal_gates.iter().any(|&g| name.starts_with(g)) {
-            return GateType::Diagonal;
-        }
-
-        let controlled_gates = [
-            "CNOT", "CZ", "SWAP", "CRx", "CRy", "CRz", "CP", "CCNOT", "CSWAP",
-        ];
-        if controlled_gates.iter().any(|&g| name.starts_with(g)) {
-            return GateType::Controlled;
-        }
-
-        if matrix.rows == 2 && matrix.cols == 2 {
-            let is_diag = matrix.data[1].real.abs() < 1e-10
-                && matrix.data[1].imaginary.abs() < 1e-10
-                && matrix.data[2].real.abs() < 1e-10
-                && matrix.data[2].imaginary.abs() < 1e-10;
-            if is_diag {
-                return GateType::Diagonal;
-            }
-        }
-
-        GateType::NonDiagonal
-    }
-
-    pub fn num_qubits(&self) -> usize {
-        self.targets.len()
-    }
-
-    pub fn target_set(&self) -> HashSet<usize> {
-        self.targets.iter().cloned().collect()
-    }
-
-    pub fn shares_qubits(&self, other: &Kernel) -> bool {
-        self.targets.iter().any(|t| other.targets.contains(t))
-    }
-
-    pub fn commutes_with(&self, other: &Kernel) -> bool {
-        if !self.shares_qubits(other) {
-            return true;
-        }
-
-        if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal {
-            if self.targets == other.targets {
-                return true;
-            }
-        }
-
-        false
-    }
-
-    pub fn can_fuse_with(&self, other: &Kernel) -> bool {
-        if self.targets.len() != 1 || other.targets.len() != 1 {
-            return false;
-        }
-        self.targets[0] == other.targets[0]
-    }
-
-    pub fn fuse(&self, other: &Kernel) -> Option<Kernel> {
-        if !self.can_fuse_with(other) {
-            return None;
-        }
-        let fused_matrix = other.matrix.dot(&self.matrix)?;
-        let new_type =
-            if self.gate_type == GateType::Diagonal && other.gate_type == GateType::Diagonal {
-                GateType::Diagonal
-            } else {
-                GateType::NonDiagonal
-            };
-        Some(Kernel {
-            matrix: fused_matrix,
-            targets: self.targets.clone(),
-            name: format!("{}+{}", self.name, other.name),
-            gate_type: new_type,
-        })
-    }
-}
-
-pub struct KernelBatch {
-    kernels: Vec<Kernel>,
-    num_qubits: usize,
-}
-
-impl KernelBatch {
-    pub fn new(num_qubits: usize) -> Self {
-        Self {
-            kernels: Vec::new(),
-            num_qubits,
-        }
-    }
-
-    pub fn add(&mut self, kernel: Kernel) {
-        self.kernels.push(kernel);
-    }
-
-    pub fn len(&self) -> usize {
-        self.kernels.len()
-    }
-
-    pub fn is_empty(&self) -> bool {
-        self.kernels.is_empty()
-    }
-
-    pub fn kernels(&self) -> &[Kernel] {
-        &self.kernels
-    }
-
-    pub fn optimize(&mut self) {
-        if self.kernels.len() < 2 {
-            return;
-        }
-
-        let mut optimized: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
-        let mut i = 0;
-
-        while i < self.kernels.len() {
-            let current = &self.kernels[i];
-
-            if i + 1 < self.kernels.len() {
-                let next = &self.kernels[i + 1];
-                if let Some(fused) = current.fuse(next) {
-                    optimized.push(fused);
-                    i += 2;
-                    continue;
-                }
-            }
-
-            optimized.push(current.clone());
-            i += 1;
-        }
-
-        self.kernels = optimized;
-    }
-
-    pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            *state = apply_kernel(state, kernel, self.num_qubits);
-        }
-    }
-
-    pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            *state = apply_kernel_parallel(state, kernel, self.num_qubits);
-        }
-    }
-
-    pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            if kernel.targets.len() == 1 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
-            } else {
-                *state = apply_kernel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            if kernel.targets.len() == 1 && self.num_qubits >= 10 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd_parallel(
-                    state,
-                    &gate,
-                    kernel.targets[0],
-                    self.num_qubits,
-                );
-            } else if kernel.targets.len() == 1 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
-            } else {
-                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn simd_capability(&self) -> SimdCapability {
-        SimdCapability::detect()
-    }
-}
-
-fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] {
-    [
-        [matrix.data[0], matrix.data[1]],
-        [matrix.data[2], matrix.data[3]],
-    ]
-}
-
-fn apply_kernel(state: &[Complex<f64>], kernel: &Kernel, num_qubits: usize) -> Vec<Complex<f64>> {
-    let dim = 1 << num_qubits;
-    let g = kernel.targets.len();
-    let gate_dim = 1 << g;
-
-    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
-    let mut non_target_mask: usize = (1 << num_qubits) - 1;
-    for &pos in &target_bits {
-        non_target_mask &= !(1 << pos);
-    }
-
-    let mut new_state = vec![complex!(0.0, 0.0); dim];
-
-    for i in 0..dim {
-        let mut target_idx = 0usize;
-        for (k, &pos) in target_bits.iter().enumerate() {
-            if (i >> pos) & 1 == 1 {
-                target_idx |= 1 << (g - 1 - k);
-            }
-        }
-
-        let mut sum = complex!(0.0, 0.0);
-
-        for j in 0..gate_dim {
-            let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
-            if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
-                continue;
-            }
-
-            let mut source_idx = i & non_target_mask;
-            for (k, &pos) in target_bits.iter().enumerate() {
-                if (j >> (g - 1 - k)) & 1 == 1 {
-                    source_idx |= 1 << pos;
-                }
-            }
-
-            sum = sum + gate_elem * state[source_idx];
-        }
-
-        new_state[i] = sum;
-    }
-
-    new_state
-}
-
-fn apply_kernel_parallel(
-    state: &[Complex<f64>],
-    kernel: &Kernel,
-    num_qubits: usize,
-) -> Vec<Complex<f64>> {
-    let dim = 1 << num_qubits;
-    let g = kernel.targets.len();
-    let gate_dim = 1 << g;
-
-    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
-    let mut non_target_mask: usize = (1 << num_qubits) - 1;
-    for &pos in &target_bits {
-        non_target_mask &= !(1 << pos);
-    }
-
-    (0..dim)
-        .into_par_iter()
-        .map(|i| {
-            let mut target_idx = 0usize;
-            for (k, &pos) in target_bits.iter().enumerate() {
-                if (i >> pos) & 1 == 1 {
-                    target_idx |= 1 << (g - 1 - k);
-                }
-            }
-
-            let mut sum = complex!(0.0, 0.0);
-
-            for j in 0..gate_dim {
-                let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
-                if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
-                    continue;
-                }
-
-                let mut source_idx = i & non_target_mask;
-                for (k, &pos) in target_bits.iter().enumerate() {
-                    if (j >> (g - 1 - k)) & 1 == 1 {
-                        source_idx |= 1 << pos;
-                    }
-                }
-
-                sum = sum + gate_elem * state[source_idx];
-            }
-
-            sum
-        })
-        .collect()
-}
-
-pub struct KernelBuilder {
-    num_qubits: usize,
-}
-
-impl KernelBuilder {
-    pub fn new(num_qubits: usize) -> Self {
-        Self { num_qubits }
-    }
-
-    pub fn num_qubits(&self) -> usize {
-        self.num_qubits
-    }
-}
-
-#[derive(Clone)]
-pub struct ExecutionLayer {
-    pub kernels: Vec<Kernel>,
-}
-
-impl ExecutionLayer {
-    pub fn new() -> Self {
-        Self {
-            kernels: Vec::new(),
-        }
-    }
-
-    pub fn can_add(&self, kernel: &Kernel) -> bool {
-        !self.kernels.iter().any(|k| k.shares_qubits(kernel))
-    }
-
-    pub fn add(&mut self, kernel: Kernel) {
-        self.kernels.push(kernel);
-    }
-
-    pub fn affected_qubits(&self) -> HashSet<usize> {
-        self.kernels
-            .iter()
-            .flat_map(|k| k.targets.iter().cloned())
-            .collect()
-    }
-}
-
-impl Default for ExecutionLayer {
-    fn default() -> Self {
-        Self::new()
-    }
-}
-
-pub struct StructureAwareKernelBatch {
-    kernels: Vec<Kernel>,
-    layers: Vec<ExecutionLayer>,
-    num_qubits: usize,
-    optimised: bool,
-}
-
-impl StructureAwareKernelBatch {
-    pub fn new(num_qubits: usize) -> Self {
-        Self {
-            kernels: Vec::new(),
-            layers: Vec::new(),
-            num_qubits,
-            optimised: false,
-        }
-    }
-
-    pub fn add(&mut self, kernel: Kernel) {
-        self.kernels.push(kernel);
-        self.optimised = false;
-    }
-
-    pub fn len(&self) -> usize {
-        self.kernels.len()
-    }
-
-    pub fn is_empty(&self) -> bool {
-        self.kernels.is_empty()
-    }
-
-    pub fn kernels(&self) -> &[Kernel] {
-        &self.kernels
-    }
-
-    pub fn layers(&self) -> &[ExecutionLayer] {
-        &self.layers
-    }
-
-    pub fn num_layers(&self) -> usize {
-        self.layers.len()
-    }
-
-    pub fn optimise(&mut self) {
-        if self.optimised || self.kernels.len() < 2 {
-            return;
-        }
-
-        self.reorder_commuting_gates();
-        self.multi_pass_fusion();
-        self.build_execution_layers();
-        self.optimised = true;
-    }
-
-    fn reorder_commuting_gates(&mut self) {
-        let mut changed = true;
-        let mut iterations = 0;
-        const MAX_ITERATIONS: usize = 100;
-
-        while changed && iterations < MAX_ITERATIONS {
-            changed = false;
-            iterations += 1;
-
-            for i in 0..self.kernels.len().saturating_sub(1) {
-                let current = &self.kernels[i];
-                let next = &self.kernels[i + 1];
-
-                if current.targets.len() == 1
-                    && next.targets.len() == 1
-                    && current.targets[0] != next.targets[0]
-                    && current.commutes_with(next)
-                {
-                    for j in (i + 2)..self.kernels.len() {
-                        let candidate = &self.kernels[j];
-
-                        if candidate.targets.len() == 1
-                            && candidate.targets[0] == current.targets[0]
-                        {
-                            let can_move = (i + 1..j).all(|k| {
-                                let between = &self.kernels[k];
-                                !between.shares_qubits(current) || current.commutes_with(between)
-                            });
-
-                            if can_move && current.can_fuse_with(candidate) {
-                                let kernel_to_move = self.kernels.remove(j);
-                                self.kernels.insert(i + 1, kernel_to_move);
-                                changed = true;
-                                break;
-                            }
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    fn multi_pass_fusion(&mut self) {
-        let mut changed = true;
-        let mut iterations = 0;
-        const MAX_ITERATIONS: usize = 50;
-
-        while changed && iterations < MAX_ITERATIONS {
-            changed = false;
-            iterations += 1;
-
-            let mut new_kernels: Vec<Kernel> = Vec::with_capacity(self.kernels.len());
-            let mut i = 0;
-
-            while i < self.kernels.len() {
-                if i + 1 < self.kernels.len() {
-                    let current = &self.kernels[i];
-                    let next = &self.kernels[i + 1];
-
-                    if let Some(fused) = current.fuse(next) {
-                        new_kernels.push(fused);
-                        i += 2;
-                        changed = true;
-                        continue;
-                    }
-                }
-
-                new_kernels.push(self.kernels[i].clone());
-                i += 1;
-            }
-
-            self.kernels = new_kernels;
-        }
-    }
-
-    fn build_execution_layers(&mut self) {
-        self.layers.clear();
-
-        for kernel in &self.kernels {
-            let mut placed = false;
-
-            for layer in &mut self.layers {
-                if layer.can_add(kernel) {
-                    layer.add(kernel.clone());
-                    placed = true;
-                    break;
-                }
-            }
-
-            if !placed {
-                let mut new_layer = ExecutionLayer::new();
-                new_layer.add(kernel.clone());
-                self.layers.push(new_layer);
-            }
-        }
-    }
-
-    pub fn execute(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            *state = apply_kernel(state, kernel, self.num_qubits);
-        }
-    }
-
-    pub fn execute_parallel(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            *state = apply_kernel_parallel(state, kernel, self.num_qubits);
-        }
-    }
-
-    pub fn execute_layered(&self, state: &mut Vec<Complex<f64>>) {
-        for layer in &self.layers {
-            for kernel in &layer.kernels {
-                *state = apply_kernel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn execute_layered_parallel(&self, state: &mut Vec<Complex<f64>>) {
-        for layer in &self.layers {
-            for kernel in &layer.kernels {
-                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn execute_simd(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            if kernel.targets.len() == 1 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
-            } else {
-                *state = apply_kernel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn execute_simd_parallel(&self, state: &mut Vec<Complex<f64>>) {
-        for kernel in &self.kernels {
-            if kernel.targets.len() == 1 && self.num_qubits >= 10 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd_parallel(
-                    state,
-                    &gate,
-                    kernel.targets[0],
-                    self.num_qubits,
-                );
-            } else if kernel.targets.len() == 1 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], self.num_qubits);
-            } else {
-                *state = apply_kernel_parallel(state, kernel, self.num_qubits);
-            }
-        }
-    }
-
-    pub fn stats(&self) -> KernelStats {
-        let single_qubit = self.kernels.iter().filter(|k| k.targets.len() == 1).count();
-        let two_qubit = self.kernels.iter().filter(|k| k.targets.len() == 2).count();
-        let multi_qubit = self.kernels.iter().filter(|k| k.targets.len() > 2).count();
-        let diagonal = self
-            .kernels
-            .iter()
-            .filter(|k| k.gate_type == GateType::Diagonal)
-            .count();
-
-        KernelStats {
-            total_kernels: self.kernels.len(),
-            single_qubit,
-            two_qubit,
-            multi_qubit,
-            diagonal,
-            execution_layers: self.layers.len(),
-        }
-    }
-}
-
-#[derive(Debug, Clone)]
-pub struct KernelStats {
-    pub total_kernels: usize,
-    pub single_qubit: usize,
-    pub two_qubit: usize,
-    pub multi_qubit: usize,
-    pub diagonal: usize,
-    pub execution_layers: usize,
-}
-
-impl std::fmt::Display for KernelStats {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(
-            f,
-            "Kernels: {} (1q: {}, 2q: {}, 3q+: {}, diag: {}), Layers: {}",
-            self.total_kernels,
-            self.single_qubit,
-            self.two_qubit,
-            self.multi_qubit,
-            self.diagonal,
-            self.execution_layers
-        )
-    }
-}
diff --git a/libpsi-core/src/core/mod.rs b/libpsi-core/src/core/mod.rs
deleted file mode 100644
index c38cf38..0000000
--- a/libpsi-core/src/core/mod.rs
+++ /dev/null
@@ -1,17 +0,0 @@
-pub mod circuit;
-pub mod classical_components;
-pub mod custom_gate;
-pub mod gates;
-pub mod kernel;
-pub mod noise;
-pub mod quantum_components;
-pub mod runtime;
-
-pub use circuit::*;
-pub use classical_components::*;
-pub use custom_gate::*;
-pub use gates::*;
-pub use kernel::*;
-pub use noise::*;
-pub use quantum_components::*;
-pub use runtime::*;
diff --git a/libpsi-core/src/core/noise.rs b/libpsi-core/src/core/noise.rs
deleted file mode 100644
index 6836c00..0000000
--- a/libpsi-core/src/core/noise.rs
+++ /dev/null
@@ -1,560 +0,0 @@
-use crate::{complex, Complex, Matrix};
-
-#[derive(Clone, Debug)]
-pub struct KrausOperator {
-    pub matrix: Matrix<Complex<f64>>,
-    pub name: String,
-}
-
-impl KrausOperator {
-    pub fn new(name: &str, matrix: Matrix<Complex<f64>>) -> Self {
-        Self {
-            matrix,
-            name: name.to_string(),
-        }
-    }
-}
-
-#[derive(Clone, Debug)]
-pub struct NoiseChannel {
-    pub name: String,
-    pub operators: Vec<KrausOperator>,
-    pub num_qubits: usize,
-}
-
-impl NoiseChannel {
-    pub fn new(name: &str, operators: Vec<KrausOperator>, num_qubits: usize) -> Self {
-        Self {
-            name: name.to_string(),
-            operators,
-            num_qubits,
-        }
-    }
-
-    pub fn depolarising(p: f64) -> Self {
-        let sqrt_1_p = (1.0 - p).sqrt();
-        let sqrt_p3 = (p / 3.0).sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_1_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(sqrt_p3, 0.0),
-                complex!(sqrt_p3, 0.0),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        let k2 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(0.0, -sqrt_p3),
-                complex!(0.0, sqrt_p3),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        let k3 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_p3, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(-sqrt_p3, 0.0),
-            ],
-        );
-
-        Self::new(
-            "Depolarising",
-            vec![
-                KrausOperator::new("K0", k0),
-                KrausOperator::new("K1(X)", k1),
-                KrausOperator::new("K2(Y)", k2),
-                KrausOperator::new("K3(Z)", k3),
-            ],
-            1,
-        )
-    }
-
-    pub fn amplitude_damping(gamma: f64) -> Self {
-        let sqrt_gamma = gamma.sqrt();
-        let sqrt_1_gamma = (1.0 - gamma).sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(1.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_gamma, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(sqrt_gamma, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        Self::new(
-            "AmplitudeDamping",
-            vec![
-                KrausOperator::new("K0", k0),
-                KrausOperator::new("K1", k1),
-            ],
-            1,
-        )
-    }
-
-    pub fn phase_damping(gamma: f64) -> Self {
-        let sqrt_gamma = gamma.sqrt();
-        let sqrt_1_gamma = (1.0 - gamma).sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(1.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_gamma, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_gamma, 0.0),
-            ],
-        );
-
-        Self::new(
-            "PhaseDamping",
-            vec![
-                KrausOperator::new("K0", k0),
-                KrausOperator::new("K1", k1),
-            ],
-            1,
-        )
-    }
-
-    pub fn bit_flip(p: f64) -> Self {
-        let sqrt_1_p = (1.0 - p).sqrt();
-        let sqrt_p = p.sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_1_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(sqrt_p, 0.0),
-                complex!(sqrt_p, 0.0),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        Self::new(
-            "BitFlip",
-            vec![
-                KrausOperator::new("K0(I)", k0),
-                KrausOperator::new("K1(X)", k1),
-            ],
-            1,
-        )
-    }
-
-    pub fn phase_flip(p: f64) -> Self {
-        let sqrt_1_p = (1.0 - p).sqrt();
-        let sqrt_p = p.sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_1_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(-sqrt_p, 0.0),
-            ],
-        );
-
-        Self::new(
-            "PhaseFlip",
-            vec![
-                KrausOperator::new("K0(I)", k0),
-                KrausOperator::new("K1(Z)", k1),
-            ],
-            1,
-        )
-    }
-
-    pub fn bit_phase_flip(p: f64) -> Self {
-        let sqrt_1_p = (1.0 - p).sqrt();
-        let sqrt_p = p.sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_1_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(0.0, -sqrt_p),
-                complex!(0.0, sqrt_p),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        Self::new(
-            "BitPhaseFlip",
-            vec![
-                KrausOperator::new("K0(I)", k0),
-                KrausOperator::new("K1(Y)", k1),
-            ],
-            1,
-        )
-    }
-
-    pub fn generalised_amplitude_damping(p: f64, gamma: f64) -> Self {
-        let sqrt_p = p.sqrt();
-        let sqrt_1_p = (1.0 - p).sqrt();
-        let sqrt_gamma = gamma.sqrt();
-        let sqrt_1_gamma = (1.0 - gamma).sqrt();
-
-        let k0 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_p, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_p * sqrt_1_gamma, 0.0),
-            ],
-        );
-
-        let k1 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(sqrt_p * sqrt_gamma, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        let k2 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(sqrt_1_p * sqrt_1_gamma, 0.0),
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p, 0.0),
-            ],
-        );
-
-        let k3 = Matrix::new(
-            2,
-            2,
-            vec![
-                complex!(0.0, 0.0),
-                complex!(0.0, 0.0),
-                complex!(sqrt_1_p * sqrt_gamma, 0.0),
-                complex!(0.0, 0.0),
-            ],
-        );
-
-        Self::new(
-            "GeneralisedAmplitudeDamping",
-            vec![
-                KrausOperator::new("K0", k0),
-                KrausOperator::new("K1", k1),
-                KrausOperator::new("K2", k2),
-                KrausOperator::new("K3", k3),
-            ],
-            1,
-        )
-    }
-}
-
-#[derive(Clone)]
-pub struct DensityMatrix {
-    pub data: Vec<Complex<f64>>,
-    pub dim: usize,
-    pub num_qubits: usize,
-}
-
-impl DensityMatrix {
-    pub fn new(num_qubits: usize) -> Self {
-        let dim = 1 << num_qubits;
-        let mut data = vec![complex!(0.0, 0.0); dim * dim];
-        data[0] = complex!(1.0, 0.0);
-        Self {
-            data,
-            dim,
-            num_qubits,
-        }
-    }
-
-    pub fn from_state_vector(state: &[Complex<f64>]) -> Self {
-        let dim = state.len();
-        let num_qubits = (dim as f64).log2() as usize;
-        let mut data = vec![complex!(0.0, 0.0); dim * dim];
-
-        for i in 0..dim {
-            for j in 0..dim {
-                data[i * dim + j] = state[i] * state[j].get_conjugate();
-            }
-        }
-
-        Self {
-            data,
-            dim,
-            num_qubits,
-        }
-    }
-
-    pub fn get(&self, row: usize, col: usize) -> Complex<f64> {
-        self.data[row * self.dim + col]
-    }
-
-    pub fn set(&mut self, row: usize, col: usize, value: Complex<f64>) {
-        self.data[row * self.dim + col] = value;
-    }
-
-    pub fn trace(&self) -> Complex<f64> {
-        let mut sum = complex!(0.0, 0.0);
-        for i in 0..self.dim {
-            sum = sum + self.get(i, i);
-        }
-        sum
-    }
-
-    pub fn purity(&self) -> f64 {
-        let mut sum = complex!(0.0, 0.0);
-        for i in 0..self.dim {
-            for j in 0..self.dim {
-                let rho_ij = self.get(i, j);
-                let rho_ji = self.get(j, i);
-                sum = sum + rho_ij * rho_ji;
-            }
-        }
-        sum.real
-    }
-
-    pub fn is_pure(&self, tolerance: f64) -> bool {
-        (self.purity() - 1.0).abs() < tolerance
-    }
-
-    pub fn probabilities(&self) -> Vec<f64> {
-        (0..self.dim).map(|i| self.get(i, i).real).collect()
-    }
-
-    pub fn apply_unitary(&mut self, gate: &Matrix<Complex<f64>>, targets: &[usize]) {
-        let g = targets.len();
-        let gate_dim = 1 << g;
-
-        let target_bits: Vec<usize> = targets
-            .iter()
-            .map(|&t| self.num_qubits - 1 - t)
-            .collect();
-
-        let mut non_target_mask: usize = (1 << self.num_qubits) - 1;
-        for &pos in &target_bits {
-            non_target_mask &= !(1 << pos);
-        }
-
-        let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim];
-
-        for i in 0..self.dim {
-            for j in 0..self.dim {
-                let mut sum = complex!(0.0, 0.0);
-
-                for k in 0..gate_dim {
-                    for l in 0..gate_dim {
-                        let mut src_i = i & non_target_mask;
-                        let mut src_j = j & non_target_mask;
-
-                        for (idx, &pos) in target_bits.iter().enumerate() {
-                            if (k >> (g - 1 - idx)) & 1 == 1 {
-                                src_i |= 1 << pos;
-                            }
-                            if (l >> (g - 1 - idx)) & 1 == 1 {
-                                src_j |= 1 << pos;
-                            }
-                        }
-
-                        let mut tgt_i = 0usize;
-                        let mut tgt_j = 0usize;
-                        for (idx, &pos) in target_bits.iter().enumerate() {
-                            if (i >> pos) & 1 == 1 {
-                                tgt_i |= 1 << (g - 1 - idx);
-                            }
-                            if (j >> pos) & 1 == 1 {
-                                tgt_j |= 1 << (g - 1 - idx);
-                            }
-                        }
-
-                        let u_ik = gate.data[tgt_i * gate_dim + k];
-                        let u_jl_dag = gate.data[tgt_j * gate_dim + l].get_conjugate();
-                        let rho_kl = self.get(src_i, src_j);
-
-                        sum = sum + u_ik * rho_kl * u_jl_dag;
-                    }
-                }
-
-                new_data[i * self.dim + j] = sum;
-            }
-        }
-
-        self.data = new_data;
-    }
-
-    pub fn apply_noise_channel(&mut self, channel: &NoiseChannel, target: usize) {
-        if channel.num_qubits != 1 {
-            panic!("Only single-qubit noise channels are currently supported");
-        }
-
-        let target_bit = self.num_qubits - 1 - target;
-        let mut new_data = vec![complex!(0.0, 0.0); self.dim * self.dim];
-
-        for kraus in &channel.operators {
-            let k = &kraus.matrix;
-
-            for i in 0..self.dim {
-                for j in 0..self.dim {
-                    let i_target = (i >> target_bit) & 1;
-                    let j_target = (j >> target_bit) & 1;
-
-                    for ki in 0..2 {
-                        for kj in 0..2 {
-                            let src_i = (i & !(1 << target_bit)) | (ki << target_bit);
-                            let src_j = (j & !(1 << target_bit)) | (kj << target_bit);
-
-                            let k_elem = k.data[i_target * 2 + ki];
-                            let k_dag_elem = k.data[j_target * 2 + kj].get_conjugate();
-                            let rho_elem = self.get(src_i, src_j);
-
-                            new_data[i * self.dim + j] =
-                                new_data[i * self.dim + j] + k_elem * rho_elem * k_dag_elem;
-                        }
-                    }
-                }
-            }
-        }
-
-        self.data = new_data;
-    }
-
-    pub fn measure_probability(&self, qubit: usize, outcome: usize) -> f64 {
-        let target_bit = self.num_qubits - 1 - qubit;
-        let mut prob = 0.0;
-
-        for i in 0..self.dim {
-            if (i >> target_bit) & 1 == outcome {
-                prob += self.get(i, i).real;
-            }
-        }
-
-        prob
-    }
-
-    pub fn fidelity_with_pure_state(&self, state: &[Complex<f64>]) -> f64 {
-        let mut sum = complex!(0.0, 0.0);
-
-        for i in 0..self.dim {
-            for j in 0..self.dim {
-                sum = sum + state[i].get_conjugate() * self.get(i, j) * state[j];
-            }
-        }
-
-        sum.real
-    }
-}
-
-impl std::fmt::Display for DensityMatrix {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        writeln!(f, "DensityMatrix ({} qubits, {}×{}):", self.num_qubits, self.dim, self.dim)?;
-        writeln!(f, "  Trace: {:.6}", self.trace().real)?;
-        writeln!(f, "  Purity: {:.6}", self.purity())?;
-        writeln!(f, "  Pure: {}", self.is_pure(1e-10))?;
-        writeln!(f, "  Probabilities: {:?}", self.probabilities())?;
-        Ok(())
-    }
-}
-
-impl std::fmt::Debug for DensityMatrix {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        writeln!(f, "DensityMatrix {}×{}:", self.dim, self.dim)?;
-        for i in 0..self.dim {
-            write!(f, "  [")?;
-            for j in 0..self.dim {
-                let val = self.get(i, j);
-                if j > 0 {
-                    write!(f, ", ")?;
-                }
-                write!(f, "{:.4}+{:.4}i", val.real, val.imaginary)?;
-            }
-            writeln!(f, "]")?;
-        }
-        Ok(())
-    }
-}
-
diff --git a/libpsi-core/src/core/quantum_components.rs b/libpsi-core/src/core/quantum_components.rs
deleted file mode 100644
index 75f6325..0000000
--- a/libpsi-core/src/core/quantum_components.rs
+++ /dev/null
@@ -1,318 +0,0 @@
-use crate::{column_vector, complex, ColumnVector, Complex, Float, Matrix, Vector, VectorMatrix};
-use core::{fmt, ops};
-
-#[macro_export]
-macro_rules! count {
-    () => { 0 };
-    ($head:expr $(,$tail:expr)*) => { 1 + count!($( $tail ),*) };
-}
-
-#[macro_export]
-macro_rules! qubit {
-    ($(($re:expr, $im:expr)),*) => {
-        {
-            let mut vector = Vec::new();
-            $(
-                vector.push(complex!($re, $im));
-            )*
-            QuantumBit::new(vector)
-        }
-    };
-}
-
-#[macro_export]
-macro_rules! quantum_register {
-    ($($bit:expr),*) => {
-        {
-            const N: usize = count!($($bit),*);
-            let mut bits: [QuantumBit; N] = [$($bit),*];
-            QuantumRegister::from(&mut bits)
-        }
-    };
-}
-
-pub type QuantumState = ColumnVector<Complex<f64>>;
-impl QuantumState {
-    pub fn state_0() -> QuantumState {
-        column_vector![complex!(1.0, 0.0), complex!(0.0, 0.0)]
-    }
-
-    pub fn state_1() -> QuantumState {
-        column_vector![complex!(0.0, 0.0), complex!(1.0, 0.0)]
-    }
-}
-
-fn identity_matrix<T: Float>(size: usize) -> Matrix<T> {
-    let mut data = vec![T::zero(); size * size];
-    for i in 0..size {
-        data[i * size + i] = T::one();
-    }
-    Matrix::new(size, size, data)
-}
-
-#[derive(Clone)]
-pub struct QuantumBit<'a> {
-    state: QuantumState,
-    name: &'a str,
-}
-
-#[derive(Clone)]
-pub struct QuantumRegister<'a> {
-    state_vector: QuantumState,
-    name: &'a str,
-    qubits: Vec<QuantumBit<'a>>,
-}
-
-#[derive(Clone)]
-pub struct QuantumGate<'a> {
-    pub name: &'a str,
-    pub matrix: Matrix<Complex<f64>>,
-    pub num_qubits: usize,
-}
-
-impl<'a> QuantumGate<'a> {
-    pub fn new(name: &'a str, matrix: Matrix<Complex<f64>>, num_qubits: usize) -> Self {
-        let expected_dim = 1 << num_qubits;
-        assert_eq!(
-            matrix.rows, expected_dim,
-            "Gate matrix rows must be 2^num_qubits"
-        );
-        assert_eq!(
-            matrix.cols, expected_dim,
-            "Gate matrix cols must be 2^num_qubits"
-        );
-        QuantumGate {
-            name,
-            matrix,
-            num_qubits,
-        }
-    }
-
-    pub fn from_matrix(name: &'a str, matrix: Matrix<Complex<f64>>) -> Self {
-        assert_eq!(matrix.rows, matrix.cols, "Gate matrix must be square");
-        let dim = matrix.rows;
-        assert!(
-            dim > 0 && (dim & (dim - 1)) == 0,
-            "Matrix dimension must be a power of 2"
-        );
-        let num_qubits = (dim as f64).log2() as usize;
-        QuantumGate {
-            name,
-            matrix,
-            num_qubits,
-        }
-    }
-}
-
-impl<'a> QuantumBit<'a> {
-    pub fn new(name: &'a str, state: QuantumState) -> QuantumBit<'a> {
-        QuantumBit { name, state }
-    }
-
-    pub fn get_state(&self) -> QuantumState {
-        self.state.clone()
-    }
-
-    pub fn get_name(&self) -> &'a str {
-        self.name
-    }
-}
-
-impl<'a> QuantumRegister<'a> {
-    pub fn new(name: &'a str, names: &[&'a str]) -> QuantumRegister<'a> {
-        let mut bits: Vec<QuantumBit<'a>> = Vec::new();
-        for i in 0..names.len() {
-            bits.push(QuantumBit::new(names[i], QuantumState::state_0()))
-        }
-
-        QuantumRegister::from(name, &mut bits)
-    }
-
-    pub fn from(name: &'a str, bits: &mut [QuantumBit<'a>]) -> QuantumRegister<'a> {
-        let mut register = QuantumRegister {
-            name,
-            qubits: bits.to_vec(),
-            state_vector: ColumnVector::new(vec![]),
-        };
-
-        register.update();
-        register
-    }
-
-    fn update(&mut self) {
-        let matrices: Vec<Matrix<Complex<f64>>> = self
-            .qubits
-            .iter()
-            .map(|qubit| qubit.state.to_matrix())
-            .collect();
-        let mut new_result = matrices[0].clone();
-        for matrix in &matrices[1..] {
-            new_result = new_result.kronecker(matrix);
-        }
-
-        self.state_vector = ColumnVector::from_matrix(&new_result);
-    }
-
-    pub fn get_bits(&self) -> Vec<QuantumBit<'_>> {
-        self.qubits.clone()
-    }
-
-    pub fn get_state(&self) -> QuantumState {
-        self.state_vector.clone()
-    }
-
-    pub fn get_name(&self) -> &'a str {
-        self.name
-    }
-
-    pub fn num_qubits(&self) -> usize {
-        self.qubits.len()
-    }
-
-    pub fn apply_gate(&mut self, gate: &QuantumGate, targets: &[usize]) {
-        let n = self.num_qubits();
-
-        assert_eq!(
-            gate.num_qubits,
-            targets.len(),
-            "Number of target qubits must match gate's qubit count"
-        );
-        for &t in targets {
-            assert!(
-                t < n,
-                "Target qubit index {} out of range for {}-qubit register",
-                t,
-                n
-            );
-        }
-
-        let mut sorted_targets = targets.to_vec();
-        sorted_targets.sort();
-        for i in 1..sorted_targets.len() {
-            assert_ne!(
-                sorted_targets[i],
-                sorted_targets[i - 1],
-                "Duplicate target qubit indices are not allowed"
-            );
-        }
-
-        let full_operator = self.build_full_operator(gate, targets);
-
-        self.state_vector = self
-            .state_vector
-            .mul_matrix(&full_operator)
-            .expect("Matrix multiplication failed during gate application");
-    }
-
-    fn build_full_operator(&self, gate: &QuantumGate, targets: &[usize]) -> Matrix<Complex<f64>> {
-        let n = self.num_qubits();
-        let g = gate.num_qubits;
-        let dim = 1 << n;
-
-        let mut contiguous = true;
-        for i in 1..targets.len() {
-            if targets[i] != targets[i - 1] + 1 {
-                contiguous = false;
-                break;
-            }
-        }
-
-        if contiguous && g == n {
-            return gate.matrix.clone();
-        }
-
-        if contiguous {
-            return self.build_contiguous_operator(gate, targets[0]);
-        }
-
-        let mut result = Matrix::new(dim, dim, vec![complex!(0.0, 0.0); dim * dim]);
-
-        for col in 0..dim {
-            for row in 0..dim {
-                let mut target_row_bits = 0usize;
-                let mut target_col_bits = 0usize;
-
-                for (i, &t) in targets.iter().enumerate() {
-                    let qubit_pos = n - 1 - t;
-                    if (row >> qubit_pos) & 1 == 1 {
-                        target_row_bits |= 1 << (g - 1 - i);
-                    }
-                    if (col >> qubit_pos) & 1 == 1 {
-                        target_col_bits |= 1 << (g - 1 - i);
-                    }
-                }
-
-                let mut non_target_match = true;
-                for q in 0..n {
-                    if !targets.contains(&q) {
-                        let qubit_pos = n - 1 - q;
-                        if ((row >> qubit_pos) & 1) != ((col >> qubit_pos) & 1) {
-                            non_target_match = false;
-                            break;
-                        }
-                    }
-                }
-
-                if non_target_match {
-                    result.set(row, col, gate.matrix.get(target_row_bits, target_col_bits));
-                }
-            }
-        }
-
-        result
-    }
-
-    fn build_contiguous_operator(
-        &self,
-        gate: &QuantumGate,
-        start_idx: usize,
-    ) -> Matrix<Complex<f64>> {
-        let n = self.num_qubits();
-        let g = gate.num_qubits;
-
-        let mut result: Option<Matrix<Complex<f64>>> = None;
-
-        for i in 0..n {
-            let part: Matrix<Complex<f64>> = if i == start_idx {
-                gate.matrix.clone()
-            } else if i > start_idx && i < start_idx + g {
-                continue;
-            } else {
-                identity_matrix(2)
-            };
-
-            result = Some(match result {
-                None => part,
-                Some(r) => r.kronecker(&part),
-            });
-        }
-
-        result.unwrap_or_else(|| identity_matrix(1 << n))
-    }
-
-    pub fn apply_gates(&mut self, operations: &[(&QuantumGate, &[usize])]) {
-        for (gate, targets) in operations {
-            self.apply_gate(gate, targets);
-        }
-    }
-}
-
-impl<'a> ops::Index<usize> for QuantumRegister<'a> {
-    type Output = QuantumBit<'a>;
-
-    fn index(&self, index: usize) -> &Self::Output {
-        &self.qubits[index]
-    }
-}
-
-impl<'a> ops::IndexMut<usize> for QuantumRegister<'a> {
-    fn index_mut(&mut self, index: usize) -> &mut Self::Output {
-        &mut self.qubits[index]
-    }
-}
-
-impl<'a> fmt::Display for QuantumGate<'a> {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "{}", self.name)
-    }
-}
diff --git a/libpsi-core/src/core/runtime.rs b/libpsi-core/src/core/runtime.rs
deleted file mode 100644
index 382fa2f..0000000
--- a/libpsi-core/src/core/runtime.rs
+++ /dev/null
@@ -1,585 +0,0 @@
-use super::{
-    GateOp, Kernel, KernelBatch, QuantumGate, QuantumRegister, QuantumState,
-    StructureAwareKernelBatch,
-};
-use crate::gates::{
-    cp_matrix, crx_matrix, cry_matrix, crz_matrix, p_matrix, rx_matrix, ry_matrix, rz_matrix,
-    u1_matrix, u2_matrix, u3_matrix, CNOT, CZ, FREDKIN, HADAMARD, PAULI_X, PAULI_Y, PAULI_Z,
-    SDG_GATE, SWAP, SXDG_GATE, SX_GATE, S_GATE, TDG_GATE, TOFFOLI, T_GATE,
-};
-use crate::maths::simd::{apply_single_qubit_gate_simd, apply_single_qubit_gate_simd_parallel};
-use crate::maths::vector::Vector;
-use crate::{complex, Complex, Matrix};
-use rayon::prelude::*;
-
-const PARALLEL_THRESHOLD: usize = 8;
-
-#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
-pub struct RuntimeConfig {
-    pub parallel: bool,
-    pub simd: bool,
-    pub batched: bool,
-    pub structure_aware: bool,
-    pub parallel_threshold: usize,
-}
-
-impl RuntimeConfig {
-    pub fn new() -> Self {
-        Self {
-            parallel: false,
-            simd: false,
-            batched: false,
-            structure_aware: false,
-            parallel_threshold: PARALLEL_THRESHOLD,
-        }
-    }
-
-    pub fn parallel(mut self) -> Self {
-        self.parallel = true;
-        self
-    }
-
-    pub fn simd(mut self) -> Self {
-        self.simd = true;
-        self
-    }
-
-    pub fn batched(mut self) -> Self {
-        self.batched = true;
-        self
-    }
-
-    pub fn structure_aware(mut self) -> Self {
-        self.structure_aware = true;
-        self
-    }
-
-    pub fn with_threshold(mut self, threshold: usize) -> Self {
-        self.parallel_threshold = threshold;
-        self
-    }
-
-    pub fn optimal() -> Self {
-        Self::new().structure_aware().simd().parallel()
-    }
-
-    pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
-        let dim = 1 << num_qubits;
-        let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
-        state[0] = complex!(1.0, 0.0);
-
-        let use_parallel = self.parallel && num_qubits >= self.parallel_threshold;
-
-        if self.structure_aware {
-            let mut batch = Runtime::build_structure_aware_batch(num_qubits, operations);
-            batch.optimise();
-            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
-        } else if self.batched {
-            let mut batch = Runtime::build_kernel_batch(num_qubits, operations);
-            batch.optimize();
-            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
-        } else {
-            let batch = Runtime::build_kernel_batch(num_qubits, operations);
-            self.execute_kernels(&mut state, batch.kernels(), num_qubits, use_parallel);
-        }
-
-        QuantumState::new(state)
-    }
-
-    fn execute_kernels(
-        &self,
-        state: &mut Vec<Complex<f64>>,
-        kernels: &[Kernel],
-        num_qubits: usize,
-        use_parallel: bool,
-    ) {
-        for kernel in kernels {
-            if self.simd && kernel.targets.len() == 1 {
-                let gate = matrix_to_2x2(&kernel.matrix);
-                if use_parallel {
-                    apply_single_qubit_gate_simd_parallel(
-                        state,
-                        &gate,
-                        kernel.targets[0],
-                        num_qubits,
-                    );
-                } else {
-                    apply_single_qubit_gate_simd(state, &gate, kernel.targets[0], num_qubits);
-                }
-            } else if use_parallel {
-                *state = apply_gate_parallel(state, &kernel.matrix, &kernel.targets, num_qubits);
-            } else {
-                *state = apply_kernel_direct(state, kernel, num_qubits);
-            }
-        }
-    }
-}
-
-impl std::fmt::Display for RuntimeConfig {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        let mut features = Vec::new();
-        if self.structure_aware {
-            features.push("structure-aware");
-        }
-        if self.batched && !self.structure_aware {
-            features.push("batched");
-        }
-        if self.simd {
-            features.push("SIMD");
-        }
-        if self.parallel {
-            features.push("parallel");
-        }
-        if features.is_empty() {
-            features.push("basic");
-        }
-        write!(f, "Runtime[{}]", features.join("+"))
-    }
-}
-
-#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
-pub enum Runtime {
-    #[default]
-    BasicRT,
-    BasicRTMT,
-    BatchedRT,
-    BatchedRTMT,
-    SimdRT,
-    SimdRTMT,
-    StructureAwareRT,
-    StructureAwareMT,
-    WFEvolution,
-    WFEvolutionMT,
-    GPUAccelerated,
-    Custom(RuntimeConfig),
-}
-
-impl Runtime {
-    pub fn custom() -> RuntimeConfig {
-        RuntimeConfig::new()
-    }
-
-    pub fn optimal() -> RuntimeConfig {
-        RuntimeConfig::optimal()
-    }
-
-    pub fn to_config(&self) -> RuntimeConfig {
-        match self {
-            Runtime::BasicRT => RuntimeConfig::new(),
-            Runtime::BasicRTMT => RuntimeConfig::new().parallel(),
-            Runtime::BatchedRT => RuntimeConfig::new().batched(),
-            Runtime::BatchedRTMT => RuntimeConfig::new().batched().parallel(),
-            Runtime::SimdRT => RuntimeConfig::new().batched().simd(),
-            Runtime::SimdRTMT => RuntimeConfig::new().batched().simd().parallel(),
-            Runtime::StructureAwareRT => RuntimeConfig::new().structure_aware().simd(),
-            Runtime::StructureAwareMT => RuntimeConfig::new().structure_aware().simd().parallel(),
-            Runtime::Custom(config) => *config,
-            _ => RuntimeConfig::new(),
-        }
-    }
-
-    pub fn compute(&self, num_qubits: usize, operations: &[GateOp]) -> QuantumState {
-        match self {
-            Runtime::BasicRT => Self::compute_basic(num_qubits, operations),
-            Runtime::BasicRTMT => Self::compute_basic_mt(num_qubits, operations),
-            Runtime::Custom(config) => config.compute(num_qubits, operations),
-            Runtime::WFEvolution => {
-                unimplemented!("WFEvolution (Schrödinger equation) runtime not yet implemented")
-            }
-            Runtime::WFEvolutionMT => {
-                unimplemented!(
-                    "WFEvolutionMT (multi-threaded Schrödinger) runtime not yet implemented"
-                )
-            }
-            Runtime::GPUAccelerated => {
-                unimplemented!("GPUAccelerated runtime not yet implemented")
-            }
-            _ => self.to_config().compute(num_qubits, operations),
-        }
-    }
-
-    pub fn build_kernel_batch(num_qubits: usize, operations: &[GateOp]) -> KernelBatch {
-        let mut batch = KernelBatch::new(num_qubits);
-
-        for op in operations {
-            if let Some(kernel) = Self::op_to_kernel(op) {
-                batch.add(kernel);
-            }
-        }
-
-        batch
-    }
-
-    fn op_to_kernel(op: &GateOp) -> Option<Kernel> {
-        let (matrix, targets, name): (Matrix<Complex<f64>>, Vec<usize>, &str) = match op {
-            GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t], "H"),
-            GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t], "X"),
-            GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t], "Y"),
-            GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t], "Z"),
-            GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t], "S"),
-            GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t], "T"),
-            GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t], "Sdg"),
-            GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t], "Tdg"),
-            GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t], "Sx"),
-            GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t], "Sxdg"),
-            GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t], "Rx"),
-            GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t], "Ry"),
-            GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t], "Rz"),
-            GateOp::P(t, theta) => (p_matrix(*theta), vec![*t], "P"),
-            GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t], "U1"),
-            GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t], "U2"),
-            GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t], "U3"),
-            GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t], "CNOT"),
-            GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t], "CZ"),
-            GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b], "SWAP"),
-            GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t], "CRx"),
-            GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t], "CRy"),
-            GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t], "CRz"),
-            GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t], "CP"),
-            GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t], "CCNOT"),
-            GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2], "CSWAP"),
-            GateOp::Measure(_, _) => return None,
-            GateOp::Custom(gate, tgts) => {
-                let qg = gate.to_quantum_gate();
-                (qg.matrix, tgts.clone(), "Custom")
-            }
-        };
-
-        Some(Kernel::new(name, matrix, targets))
-    }
-
-    pub fn build_structure_aware_batch(
-        num_qubits: usize,
-        operations: &[GateOp],
-    ) -> StructureAwareKernelBatch {
-        let mut batch = StructureAwareKernelBatch::new(num_qubits);
-
-        for op in operations {
-            if let Some(kernel) = Self::op_to_kernel(op) {
-                batch.add(kernel);
-            }
-        }
-
-        batch
-    }
-
-    fn compute_basic(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
-        let names: Vec<String> = (0..num_qubits).map(|i| format!("q{}", i)).collect();
-        let leaked_names: &'static [String] = Box::leak(names.into_boxed_slice());
-        let name_refs: Vec<&'static str> = leaked_names.iter().map(|s| s.as_str()).collect();
-
-        let mut register = QuantumRegister::new(
-            Box::leak(Box::new("circuit".to_string())).as_str(),
-            &name_refs,
-        );
-
-        for op in operations {
-            match op {
-                // Clifford gates
-                GateOp::H(t) => register.apply_gate(&HADAMARD, &[*t]),
-                GateOp::X(t) => register.apply_gate(&PAULI_X, &[*t]),
-                GateOp::Y(t) => register.apply_gate(&PAULI_Y, &[*t]),
-                GateOp::Z(t) => register.apply_gate(&PAULI_Z, &[*t]),
-                GateOp::S(t) => register.apply_gate(&S_GATE, &[*t]),
-                GateOp::CNOT(c, t) => register.apply_gate(&CNOT, &[*c, *t]),
-                GateOp::CZ(c, t) => register.apply_gate(&CZ, &[*c, *t]),
-                GateOp::SWAP(a, b) => register.apply_gate(&SWAP, &[*a, *b]),
-                GateOp::CCNOT(c1, c2, t) => register.apply_gate(&TOFFOLI, &[*c1, *c2, *t]),
-                GateOp::CSWAP(c, t1, t2) => register.apply_gate(&FREDKIN, &[*c, *t1, *t2]),
-
-                // Non-Clifford fixed gates
-                GateOp::T(t) => register.apply_gate(&T_GATE, &[*t]),
-                GateOp::Sdg(t) => register.apply_gate(&SDG_GATE, &[*t]),
-                GateOp::Tdg(t) => register.apply_gate(&TDG_GATE, &[*t]),
-                GateOp::Sx(t) => register.apply_gate(&SX_GATE, &[*t]),
-                GateOp::Sxdg(t) => register.apply_gate(&SXDG_GATE, &[*t]),
-
-                // Parametric single-qubit gates (non-Clifford for most angles)
-                GateOp::Rx(t, theta) => {
-                    let gate = QuantumGate {
-                        name: "Rx",
-                        matrix: rx_matrix(*theta),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::Ry(t, theta) => {
-                    let gate = QuantumGate {
-                        name: "Ry",
-                        matrix: ry_matrix(*theta),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::Rz(t, theta) => {
-                    let gate = QuantumGate {
-                        name: "Rz",
-                        matrix: rz_matrix(*theta),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::P(t, theta) => {
-                    let gate = QuantumGate {
-                        name: "P",
-                        matrix: p_matrix(*theta),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::U1(t, lambda) => {
-                    let gate = QuantumGate {
-                        name: "U1",
-                        matrix: u1_matrix(*lambda),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::U2(t, phi, lambda) => {
-                    let gate = QuantumGate {
-                        name: "U2",
-                        matrix: u2_matrix(*phi, *lambda),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-                GateOp::U3(t, theta, phi, lambda) => {
-                    let gate = QuantumGate {
-                        name: "U3",
-                        matrix: u3_matrix(*theta, *phi, *lambda),
-                        num_qubits: 1,
-                    };
-                    register.apply_gate(&gate, &[*t]);
-                }
-
-                // Controlled parametric gates
-                GateOp::CRx(c, t, theta) => {
-                    let gate = QuantumGate {
-                        name: "CRx",
-                        matrix: crx_matrix(*theta),
-                        num_qubits: 2,
-                    };
-                    register.apply_gate(&gate, &[*c, *t]);
-                }
-                GateOp::CRy(c, t, theta) => {
-                    let gate = QuantumGate {
-                        name: "CRy",
-                        matrix: cry_matrix(*theta),
-                        num_qubits: 2,
-                    };
-                    register.apply_gate(&gate, &[*c, *t]);
-                }
-                GateOp::CRz(c, t, theta) => {
-                    let gate = QuantumGate {
-                        name: "CRz",
-                        matrix: crz_matrix(*theta),
-                        num_qubits: 2,
-                    };
-                    register.apply_gate(&gate, &[*c, *t]);
-                }
-                GateOp::CP(c, t, theta) => {
-                    let gate = QuantumGate {
-                        name: "CP",
-                        matrix: cp_matrix(*theta),
-                        num_qubits: 2,
-                    };
-                    register.apply_gate(&gate, &[*c, *t]);
-                }
-
-                // Measurement and custom gates
-                GateOp::Measure(_, _) => {}
-                GateOp::Custom(gate, targets) => {
-                    let quantum_gate = gate.to_quantum_gate();
-                    register.apply_gate(&quantum_gate, targets);
-                }
-            }
-        }
-
-        register.get_state()
-    }
-
-    fn compute_basic_mt(num_qubits: usize, operations: &[GateOp]) -> QuantumState {
-        // For small circuits, fall back to single-threaded (overhead not worth it)
-        if num_qubits < PARALLEL_THRESHOLD {
-            return Self::compute_basic(num_qubits, operations);
-        }
-
-        let dim = 1 << num_qubits;
-
-        // Initialize state to |0...0⟩
-        let mut state: Vec<Complex<f64>> = vec![complex!(0.0, 0.0); dim];
-        state[0] = complex!(1.0, 0.0);
-
-        for op in operations {
-            let (gate_matrix, targets): (Matrix<Complex<f64>>, Vec<usize>) = match op {
-                // Clifford gates
-                GateOp::H(t) => (HADAMARD.matrix.clone(), vec![*t]),
-                GateOp::X(t) => (PAULI_X.matrix.clone(), vec![*t]),
-                GateOp::Y(t) => (PAULI_Y.matrix.clone(), vec![*t]),
-                GateOp::Z(t) => (PAULI_Z.matrix.clone(), vec![*t]),
-                GateOp::S(t) => (S_GATE.matrix.clone(), vec![*t]),
-                GateOp::CNOT(c, t) => (CNOT.matrix.clone(), vec![*c, *t]),
-                GateOp::CZ(c, t) => (CZ.matrix.clone(), vec![*c, *t]),
-                GateOp::SWAP(a, b) => (SWAP.matrix.clone(), vec![*a, *b]),
-                GateOp::CCNOT(c1, c2, t) => (TOFFOLI.matrix.clone(), vec![*c1, *c2, *t]),
-                GateOp::CSWAP(c, t1, t2) => (FREDKIN.matrix.clone(), vec![*c, *t1, *t2]),
-
-                // Non-Clifford fixed gates
-                GateOp::T(t) => (T_GATE.matrix.clone(), vec![*t]),
-                GateOp::Sdg(t) => (SDG_GATE.matrix.clone(), vec![*t]),
-                GateOp::Tdg(t) => (TDG_GATE.matrix.clone(), vec![*t]),
-                GateOp::Sx(t) => (SX_GATE.matrix.clone(), vec![*t]),
-                GateOp::Sxdg(t) => (SXDG_GATE.matrix.clone(), vec![*t]),
-
-                // Parametric single-qubit gates
-                GateOp::Rx(t, theta) => (rx_matrix(*theta), vec![*t]),
-                GateOp::Ry(t, theta) => (ry_matrix(*theta), vec![*t]),
-                GateOp::Rz(t, theta) => (rz_matrix(*theta), vec![*t]),
-                GateOp::P(t, theta) => (p_matrix(*theta), vec![*t]),
-                GateOp::U1(t, lambda) => (u1_matrix(*lambda), vec![*t]),
-                GateOp::U2(t, phi, lambda) => (u2_matrix(*phi, *lambda), vec![*t]),
-                GateOp::U3(t, theta, phi, lambda) => (u3_matrix(*theta, *phi, *lambda), vec![*t]),
-
-                // Controlled parametric gates
-                GateOp::CRx(c, t, theta) => (crx_matrix(*theta), vec![*c, *t]),
-                GateOp::CRy(c, t, theta) => (cry_matrix(*theta), vec![*c, *t]),
-                GateOp::CRz(c, t, theta) => (crz_matrix(*theta), vec![*c, *t]),
-                GateOp::CP(c, t, theta) => (cp_matrix(*theta), vec![*c, *t]),
-
-                // Measurement (skip) and custom gates
-                GateOp::Measure(_, _) => continue,
-                GateOp::Custom(custom_gate, tgts) => {
-                    let quantum_gate = custom_gate.to_quantum_gate();
-                    state = apply_gate_parallel(&state, &quantum_gate.matrix, tgts, num_qubits);
-                    continue;
-                }
-            };
-
-            state = apply_gate_parallel(&state, &gate_matrix, &targets, num_qubits);
-        }
-
-        QuantumState::new(state)
-    }
-}
-
-/// Apply a gate to the state vector in parallel using sparse application
-/// This is O(2^n * 2^g) instead of O(2^2n) for full matrix multiplication
-fn apply_gate_parallel(
-    state: &[Complex<f64>],
-    gate_matrix: &Matrix<Complex<f64>>,
-    targets: &[usize],
-    num_qubits: usize,
-) -> Vec<Complex<f64>> {
-    let dim = 1 << num_qubits;
-    let g = targets.len();
-    let gate_dim = 1 << g;
-
-    // Convert target qubit indices to bit positions (from MSB)
-    let target_bits: Vec<usize> = targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
-    // Create a mask for non-target qubits
-    let mut non_target_mask: usize = (1 << num_qubits) - 1;
-    for &pos in &target_bits {
-        non_target_mask &= !(1 << pos);
-    }
-
-    // Parallel computation of new state
-    let new_state: Vec<Complex<f64>> = (0..dim)
-        .into_par_iter()
-        .map(|i| {
-            // Extract the target qubit bits from index i
-            let mut target_idx = 0usize;
-            for (k, &pos) in target_bits.iter().enumerate() {
-                if (i >> pos) & 1 == 1 {
-                    target_idx |= 1 << (g - 1 - k);
-                }
-            }
-
-            // Compute the contribution to state[i]
-            let mut sum = complex!(0.0, 0.0);
-
-            // For each possible input state that could contribute
-            for j in 0..gate_dim {
-                // Get the gate matrix element
-                let gate_elem = gate_matrix.data[target_idx * gate_dim + j];
-
-                // Skip if zero (sparse optimization)
-                if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
-                    continue;
-                }
-
-                // Compute the source index by replacing target bits in i with bits from j
-                let mut source_idx = i & non_target_mask;
-                for (k, &pos) in target_bits.iter().enumerate() {
-                    if (j >> (g - 1 - k)) & 1 == 1 {
-                        source_idx |= 1 << pos;
-                    }
-                }
-
-                sum = sum + gate_elem * state[source_idx];
-            }
-
-            sum
-        })
-        .collect();
-
-    new_state
-}
-
-fn matrix_to_2x2(matrix: &Matrix<Complex<f64>>) -> [[Complex<f64>; 2]; 2] {
-    [
-        [matrix.data[0], matrix.data[1]],
-        [matrix.data[2], matrix.data[3]],
-    ]
-}
-
-fn apply_kernel_direct(
-    state: &[Complex<f64>],
-    kernel: &Kernel,
-    num_qubits: usize,
-) -> Vec<Complex<f64>> {
-    let dim = 1 << num_qubits;
-    let g = kernel.targets.len();
-    let gate_dim = 1 << g;
-
-    let target_bits: Vec<usize> = kernel.targets.iter().map(|&t| num_qubits - 1 - t).collect();
-
-    let mut non_target_mask: usize = (1 << num_qubits) - 1;
-    for &pos in &target_bits {
-        non_target_mask &= !(1 << pos);
-    }
-
-    let mut new_state = vec![complex!(0.0, 0.0); dim];
-
-    for i in 0..dim {
-        let mut target_idx = 0usize;
-        for (k, &pos) in target_bits.iter().enumerate() {
-            if (i >> pos) & 1 == 1 {
-                target_idx |= 1 << (g - 1 - k);
-            }
-        }
-
-        let mut sum = complex!(0.0, 0.0);
-
-        for j in 0..gate_dim {
-            let gate_elem = kernel.matrix.data[target_idx * gate_dim + j];
-
-            if gate_elem.real.abs() < 1e-15 && gate_elem.imaginary.abs() < 1e-15 {
-                continue;
-            }
-
-            let mut source_idx = i & non_target_mask;
-            for (k, &pos) in target_bits.iter().enumerate() {
-                if (j >> (g - 1 - k)) & 1 == 1 {
-                    source_idx |= 1 << pos;
-                }
-            }
-
-            sum = sum + gate_elem * state[source_idx];
-        }
-
-        new_state[i] = sum;
-    }
-
-    new_state
-}