diff options
| -rw-r--r-- | README.md | 28 | ||||
| -rw-r--r-- | libpsi-core/src/core/mod.rs | 2 | ||||
| -rw-r--r-- | libpsi-core/src/core/noise.rs | 560 | ||||
| -rw-r--r-- | libpsi-core/src/lib.rs | 1 | ||||
| -rw-r--r-- | tester/src/main.rs | 8 | ||||
| -rw-r--r-- | tester/src/noise.rs | 172 |
6 files changed, 771 insertions, 0 deletions
@@ -61,6 +61,34 @@ Automatic detection and use of platform-specific SIMD instructions: - Multi-pass fusion until convergence - Execution layer grouping for parallelism +### Noise Channels (Density Matrix) + +Realistic quantum noise simulation using Kraus operators: + +| Channel | Description | +|---------|-------------| +| `depolarising(p)` | Random Pauli error with probability $p$ | +| `amplitude_damping(γ)` | Energy decay ($T_1$ relaxation) | +| `phase_damping(γ)` | Phase decoherence ($T_2$ dephasing) | +| `bit_flip(p)` | $X$ error with probability $p$ | +| `phase_flip(p)` | $Z$ error with probability $p$ | +| `bit_phase_flip(p)` | $Y$ error with probability $p$ | + +```rust +use libpsi_core::{DensityMatrix, NoiseChannel}; + +// Create density matrix from circuit state +let dm = DensityMatrix::from_state_vector(&state_vec); + +// Apply noise +let noise = NoiseChannel::depolarising(0.05); +dm.apply_noise_channel(&noise, 0); // Apply to qubit 0 + +// Check properties +println!("Purity: {}", dm.purity()); // 1.0 = pure, <1.0 = mixed +println!("Fidelity: {}", dm.fidelity_with_pure_state(&ideal_state)); +``` + ## Project Structure - **`libpsi-core`**: Core quantum simulation library diff --git a/libpsi-core/src/core/mod.rs b/libpsi-core/src/core/mod.rs index 343e319..c38cf38 100644 --- a/libpsi-core/src/core/mod.rs +++ b/libpsi-core/src/core/mod.rs @@ -3,6 +3,7 @@ pub mod classical_components; pub mod custom_gate; pub mod gates; pub mod kernel; +pub mod noise; pub mod quantum_components; pub mod runtime; @@ -11,5 +12,6 @@ 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 new file mode 100644 index 0000000..6836c00 --- /dev/null +++ b/libpsi-core/src/core/noise.rs @@ -0,0 +1,560 @@ +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/lib.rs b/libpsi-core/src/lib.rs index 96c77a0..5876918 100644 --- a/libpsi-core/src/lib.rs +++ b/libpsi-core/src/lib.rs @@ -13,5 +13,6 @@ pub use core::classical_components::*; pub use core::custom_gate::*; pub use core::gates; pub use core::kernel::*; +pub use core::noise::*; pub use core::quantum_components::*; pub use core::runtime::*; diff --git a/tester/src/main.rs b/tester/src/main.rs index 854472c..8a9a41c 100644 --- a/tester/src/main.rs +++ b/tester/src/main.rs @@ -3,6 +3,7 @@ mod clifford; mod common; mod custom_gates; mod kernels; +mod noise; mod non_clifford; mod simd; @@ -25,6 +26,7 @@ fn print_usage() { println!(" custom Run custom gate tests only"); println!(" kernels Run kernel batching tests only"); println!(" simd Run SIMD acceleration tests only"); + println!(" noise Run noise channel tests only"); println!(" bench Run benchmark tests only"); println!(" help Show this help message"); println!(); @@ -34,6 +36,7 @@ fn print_usage() { println!(" tester non-clifford # Run only rotation/parametric gate tests"); println!(" tester kernels # Run only kernel batching tests"); println!(" tester simd # Run only SIMD tests"); + println!(" tester noise # Run only noise channel tests"); println!(" tester custom bench # Run custom gates and benchmarks"); } @@ -58,6 +61,7 @@ fn main() { let run_custom = run_all || args.iter().any(|a| a == "custom"); let run_kernels = run_all || args.iter().any(|a| a == "kernels"); let run_simd = run_all || args.iter().any(|a| a == "simd"); + let run_noise = run_all || args.iter().any(|a| a == "noise"); let run_bench = run_all || args.iter().any(|a| a == "bench"); if run_clifford { @@ -80,6 +84,10 @@ fn main() { simd::run_all(&mut results); } + if run_noise { + noise::run_all(&mut results); + } + if run_bench { benchmarks::run_all(&mut results); } diff --git a/tester/src/noise.rs b/tester/src/noise.rs new file mode 100644 index 0000000..54edfe0 --- /dev/null +++ b/tester/src/noise.rs @@ -0,0 +1,172 @@ +use crate::common::{print_section, BenchmarkResult}; +use libpsi_core::{ + complex, DensityMatrix, NoiseChannel, QuantumCircuit, Runtime, Vector, +}; +use std::time::Instant; + +pub fn run_all(results: &mut Vec<BenchmarkResult>) { + println!("═══════════════════════════════════════════════════════════════"); + println!(" NOISE CHANNEL TESTS"); + println!("═══════════════════════════════════════════════════════════════\n"); + + test_density_matrix_basics(results); + test_noise_channels(results); + test_noisy_circuit(results); +} + +pub fn test_density_matrix_basics(results: &mut Vec<BenchmarkResult>) { + print_section("Density Matrix Basics"); + + let dm = DensityMatrix::new(2); + println!("Initial |00⟩ state:"); + println!("{}", dm); + + let mut circuit = QuantumCircuit::new(2); + circuit.h(0).cnot(0, 1); + circuit.compute_with(Runtime::BasicRT); + let state = circuit.state(); + + let state_vec: Vec<_> = (0..state.size()) + .map(|i| state.get(i)) + .collect(); + + let dm_bell = DensityMatrix::from_state_vector(&state_vec); + println!("Bell state |Φ+⟩:"); + println!("{}", dm_bell); + println!("Full matrix:"); + println!("{:?}", dm_bell); + + let is_pure = dm_bell.is_pure(1e-10); + println!("Purity check: {}\n", if is_pure { "✓ Pure" } else { "✗ Mixed" }); + + results.push(BenchmarkResult { + name: "DM: Bell state".to_string(), + basic_time: std::time::Duration::from_micros(0), + mt_time: std::time::Duration::from_micros(0), + results_match: is_pure, + }); +} + +pub fn test_noise_channels(results: &mut Vec<BenchmarkResult>) { + print_section("Noise Channel Effects"); + + let channels: Vec<(&str, NoiseChannel)> = vec![ + ("Depolarising (p=0.1)", NoiseChannel::depolarising(0.1)), + ("Amplitude Damping (γ=0.2)", NoiseChannel::amplitude_damping(0.2)), + ("Phase Damping (γ=0.2)", NoiseChannel::phase_damping(0.2)), + ("Bit Flip (p=0.1)", NoiseChannel::bit_flip(0.1)), + ("Phase Flip (p=0.1)", NoiseChannel::phase_flip(0.1)), + ("Bit-Phase Flip (p=0.1)", NoiseChannel::bit_phase_flip(0.1)), + ]; + + let plus_state = vec![ + complex!(1.0 / 2.0_f64.sqrt(), 0.0), + complex!(1.0 / 2.0_f64.sqrt(), 0.0), + ]; + + println!("Starting with |+⟩ state: (|0⟩ + |1⟩)/√2\n"); + + for (name, channel) in channels { + let mut dm = DensityMatrix::from_state_vector(&plus_state); + let initial_purity = dm.purity(); + + let start = Instant::now(); + dm.apply_noise_channel(&channel, 0); + let elapsed = start.elapsed(); + + let final_purity = dm.purity(); + let fidelity = dm.fidelity_with_pure_state(&plus_state); + + println!("{:30}", name); + println!(" Purity: {:.4} → {:.4}", initial_purity, final_purity); + println!(" Fidelity with |+⟩: {:.4}", fidelity); + println!(" Probabilities: {:?}", dm.probabilities()); + println!(" Time: {:.2}μs\n", elapsed.as_secs_f64() * 1_000_000.0); + + let purity_decreased = final_purity <= initial_purity + 1e-10; + + results.push(BenchmarkResult { + name: format!("Noise: {}", name), + basic_time: elapsed, + mt_time: elapsed, + results_match: purity_decreased, + }); + } +} + +pub fn test_noisy_circuit(results: &mut Vec<BenchmarkResult>) { + print_section("Noisy Circuit Simulation"); + + let mut circuit = QuantumCircuit::new(2); + circuit.h(0).cnot(0, 1); + circuit.compute_with(Runtime::BasicRT); + let state = circuit.state(); + let state_vec: Vec<_> = (0..state.size()).map(|i| state.get(i)).collect(); + + let mut dm = DensityMatrix::from_state_vector(&state_vec); + println!("Bell state before noise:"); + println!("{}", dm); + + let depol = NoiseChannel::depolarising(0.05); + + let start = Instant::now(); + dm.apply_noise_channel(&depol, 0); + dm.apply_noise_channel(&depol, 1); + let elapsed = start.elapsed(); + + println!("Bell state after 5% depolarising on both qubits:"); + println!("{}", dm); + + let fidelity = dm.fidelity_with_pure_state(&state_vec); + println!("Fidelity with ideal Bell state: {:.4}", fidelity); + println!("Time: {:.2}μs\n", elapsed.as_secs_f64() * 1_000_000.0); + + let mut dm2 = DensityMatrix::from_state_vector(&state_vec); + let amp_damp = NoiseChannel::amplitude_damping(0.1); + + dm2.apply_noise_channel(&_damp, 0); + dm2.apply_noise_channel(&_damp, 1); + + println!("Bell state after 10% amplitude damping on both qubits:"); + println!("{}", dm2); + println!("Probabilities show decay towards |00⟩: {:?}", dm2.probabilities()); + + results.push(BenchmarkResult { + name: "Noisy Bell circuit".to_string(), + basic_time: elapsed, + mt_time: elapsed, + results_match: fidelity > 0.8 && fidelity < 1.0, + }); + + println!(); + print_section("T1/T2 Relaxation Simulation"); + + let one_state = vec![complex!(0.0, 0.0), complex!(1.0, 0.0)]; + let mut dm_t1 = DensityMatrix::from_state_vector(&one_state); + + println!("Simulating T1 decay of |1⟩ state:"); + println!(" Initial: P(0)={:.4}, P(1)={:.4}", dm_t1.probabilities()[0], dm_t1.probabilities()[1]); + + let t1_channel = NoiseChannel::amplitude_damping(0.3); + for step in 1..=5 { + dm_t1.apply_noise_channel(&t1_channel, 0); + println!( + " Step {}: P(0)={:.4}, P(1)={:.4}, Purity={:.4}", + step, + dm_t1.probabilities()[0], + dm_t1.probabilities()[1], + dm_t1.purity() + ); + } + + let decayed = dm_t1.probabilities()[0] > 0.8; + println!(" Decay complete: {}\n", if decayed { "✓" } else { "✗" }); + + results.push(BenchmarkResult { + name: "T1 decay simulation".to_string(), + basic_time: std::time::Duration::from_micros(0), + mt_time: std::time::Duration::from_micros(0), + results_match: decayed, + }); +} + |