use std::{collections::HashMap, ops::BitXorAssign};
use bitvec::vec::BitVec;
use pauli::Pauli;
use rand::{rngs::ThreadRng, thread_rng, Rng};
use rayon::prelude::*;
use crate::{
circuit::Circuit,
counts::Counts,
instruction::{CliffordKind, Instruction, PauliKind},
simulator_traits::Simulator,
};
pub mod pauli;
#[derive(Clone)]
pub struct Stabilizer {
pub stabilizer_pauli: Vec<pauli::Pauli>,
pub destabilizer_pauli: Vec<pauli::Pauli>,
pub stabilizer_phase: BitVec<u64>,
pub destabilizer_phase: BitVec<u64>,
pub n: usize,
pub rng: ThreadRng,
pub reg: BitVec,
}
fn set_bit(value: &mut u64, bit: usize, state: bool) {
*value = (*value & !(1 << bit)) | ((state as u64) << bit);
}
impl Stabilizer {
pub fn new(n: usize) -> Self {
let mut stabilizer_pauli: Vec<Pauli> = (0..n).map(|_| Pauli::new(n)).collect();
let mut destabilizer_pauli: Vec<Pauli> = (0..n).map(|_| Pauli::new(n)).collect();
for i in 0..n {
destabilizer_pauli[i].x.resize(n, false);
destabilizer_pauli[i].z.resize(n, false);
destabilizer_pauli[i].x.set(i, true);
stabilizer_pauli[i].x.resize(n, false);
stabilizer_pauli[i].z.resize(n, false);
stabilizer_pauli[i].z.set(i, true);
}
let mut stabilizer_phase: BitVec<u64> = BitVec::new();
let mut destabilizer_phase: BitVec<u64> = BitVec::new();
stabilizer_phase.resize(n, false);
destabilizer_phase.resize(n, false);
let mut reg = BitVec::new();
reg.resize(n, false);
Self {
stabilizer_pauli,
destabilizer_pauli,
stabilizer_phase,
destabilizer_phase,
n,
rng: thread_rng(),
reg,
}
}
pub fn anticommutes_row(&self, r: usize, s: usize) -> bool {
let rz = &self.stabilizer_pauli[r].z;
let rx = &self.stabilizer_pauli[r].x;
let sz = &self.stabilizer_pauli[s].z;
let sx = &self.stabilizer_pauli[s].x;
(0..self.n)
.into_iter()
.map(|i| (rx[i] as u8 * sz[i] as u8 + rz[i] as u8 * sx[i] as u8) % 2)
.sum::<u8>()
% 2
== 1
}
pub fn anticommutes(&self, i: usize) -> Option<usize> {
self.stabilizer_pauli
.iter()
.enumerate()
.filter(|(_index, r)| *r.x.get(i).unwrap())
.map(|(index, _x)| index)
.min()
}
pub fn rowsum(&mut self, h: usize, i: usize) {
let g = |x1: i8, z1: i8, x2: i8, z2: i8| {
if x1 == z1 && z1 == 0 {
return 0;
}
if x1 == z1 && z1 == 1 {
return z2 - x2;
}
if x1 == 1 && z1 == 0 {
return z2 * (2 * x2 - 1);
}
if x1 == 0 && z1 == 1 {
return x2 * (1 - 2 * z2);
}
panic!();
};
let h_phase = if h < self.n {
self.destabilizer_phase[h] as i8
} else {
self.stabilizer_phase[h - self.n] as i8
};
let h_row: &mut Pauli = if h < self.n {
unsafe { &mut *(&mut self.destabilizer_pauli[h] as *mut _) }
} else {
unsafe { &mut *(&mut self.stabilizer_pauli[h - self.n] as *mut _) }
};
let sum = 2 * h_phase
+ 2 * self.stabilizer_phase[i] as i8
+ (0..self.n)
.map(|j| {
g(
self.stabilizer_pauli[i].x[j] as i8,
self.stabilizer_pauli[i].z[j] as i8,
h_row.x[j] as i8,
h_row.z[j] as i8,
)
})
.sum::<i8>();
let cond = sum % 2 != 0;
if h < self.n {
self.destabilizer_phase.set(h, cond);
} else {
self.stabilizer_phase.set(h - self.n, cond);
}
(0..self.n).for_each(|j| {
let fml = self.stabilizer_pauli[i].x[j];
h_row.x.get_mut(j).unwrap().bitxor_assign(fml)
});
(0..self.n).for_each(|j| {
let fml = self.stabilizer_pauli[i].z[j];
h_row.z.get_mut(j).unwrap().bitxor_assign(fml)
});
}
pub fn measure_sq(&mut self, target: usize) -> bool {
let p = self.anticommutes(target);
if let Some(p) = p {
for i in 0..self.n {
if i == p || !self.stabilizer_pauli[i].x[target] {
continue;
}
self.rowsum(p, i);
self.destabilizer_pauli[p] = self.stabilizer_pauli[p].clone();
}
self.stabilizer_pauli[p].x.fill(false);
self.stabilizer_pauli[p].z.fill(false);
self.stabilizer_pauli[p].z.set(target, true);
let outcome = self.rng.gen_bool(0.5);
self.stabilizer_phase.set(p, outcome);
outcome
} else {
let mut accum = false;
for i in 0..self.n {
if self.anticommutes_row(target, i) {
accum ^= self.destabilizer_phase[i];
}
}
accum
}
}
pub fn print_tableu(&self) {
println!();
for (index, i) in self.destabilizer_pauli.iter().enumerate() {
print!(
"{:?} | {:?} | {}",
i.x.iter()
.map(|x| if *x { 1 } else { 0 })
.collect::<Vec<u8>>(),
i.z.iter()
.map(|x| if *x { 1 } else { 0 })
.collect::<Vec<u8>>(),
self.destabilizer_phase[index] as u8
);
println!();
}
println!();
for (index, i) in self.stabilizer_pauli.iter().enumerate() {
print!(
"{:?} | {:?} | {}",
i.x.iter()
.map(|x| if *x { 1 } else { 0 })
.collect::<Vec<u8>>(),
i.z.iter()
.map(|x| if *x { 1 } else { 0 })
.collect::<Vec<u8>>(),
self.stabilizer_phase[index] as u8
);
println!();
}
}
pub fn execute_instruction(&mut self, ins: &Instruction) {
match ins {
Instruction::Pauli { kind, target } => match kind {
PauliKind::X => {
self.destabilizer_pauli
.chunks(64)
.chain(self.stabilizer_pauli.chunks(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let z_mask = rows.iter().enumerate().fold(0_u64, |acc, (i, x)| {
let mut j = acc;
set_bit(&mut j, i, x.z[*target]);
j
});
*phases ^= z_mask;
});
}
PauliKind::Z => {
self.destabilizer_pauli
.chunks(64)
.chain(self.stabilizer_pauli.chunks(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let x_mask = rows.iter().enumerate().fold(0_u64, |acc, (i, x)| {
let mut j = acc;
set_bit(&mut j, i, x.x[*target]);
j
});
*phases ^= x_mask;
});
}
PauliKind::Y => {
self.destabilizer_pauli
.chunks(64)
.chain(self.stabilizer_pauli.chunks(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let xz_mask = rows.iter().enumerate().fold(0_u64, |acc, (i, r)| {
let mut j = acc;
set_bit(&mut j, i, r.z[*target] ^ r.x[*target]);
j
});
*phases ^= xz_mask;
});
}
},
Instruction::Clifford { kind, target } => match kind {
CliffordKind::H => {
let qubit = target.to_single();
self.destabilizer_pauli
.chunks_mut(64)
.chain(self.stabilizer_pauli.chunks_mut(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let xz_mask = rows.iter_mut().enumerate().fold(0_u64, |acc, (i, r)| {
let mut j = acc;
set_bit(&mut j, i, r.z[qubit] & r.x[qubit]);
let temp = r.x[qubit];
let ttemp = r.z[qubit];
r.x.set(qubit, ttemp);
r.z.set(qubit, temp);
j
});
*phases ^= xz_mask;
});
}
CliffordKind::S => {
let qubit = target.to_single();
self.destabilizer_pauli
.chunks_mut(64)
.chain(self.stabilizer_pauli.chunks_mut(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let xz_mask = rows.iter_mut().enumerate().fold(0_u64, |acc, (i, r)| {
let mut j = acc;
set_bit(&mut j, i, r.z[qubit] & r.x[qubit]);
let temp = r.x[qubit];
r.z.get_mut(qubit).unwrap().bitxor_assign(temp);
j
});
*phases ^= xz_mask;
});
}
CliffordKind::CNOT => {
let (a, b) = target.to_two();
self.destabilizer_pauli
.chunks_mut(64)
.chain(self.stabilizer_pauli.chunks_mut(64))
.zip(
self.destabilizer_phase
.as_raw_mut_slice()
.iter_mut()
.chain(self.stabilizer_phase.as_raw_mut_slice().iter_mut()),
)
.for_each(|(rows, phases)| {
let xz_mask = rows.iter_mut().enumerate().fold(0_u64, |acc, (i, r)| {
let mut j = acc;
set_bit(&mut j, i, r.x[a] & r.z[b] & (r.x[b] ^ r.z[a] ^ true));
let xr = r.x[b] ^ r.z[a];
r.x.get_mut(b).unwrap().bitxor_assign(xr);
let zr = r.z[a] ^ r.z[b];
r.z.get_mut(a).unwrap().bitxor_assign(zr);
j
});
*phases ^= xz_mask;
});
}
},
Instruction::Gate { .. } => {
panic!("stabilizer simulator only supports pauli and clifford gates")
}
Instruction::StateVecMat { .. } => {
panic!("StateVecMat instructions are only supported on the statevector simulator");
}
Instruction::Identity => {}
Instruction::Measure { indices } => {
for i in indices {
let r = self.measure_sq(*i);
self.reg.set(*i, r);
}
}
Instruction::MeasureNoisy { indices, p } => {
for i in indices {
let mut r = self.measure_sq(*i);
if self.rng.gen_range(0.0..1.0) < *p {
r = !r;
self.execute_instruction(&Instruction::Pauli {
kind: PauliKind::X,
target: *i,
});
}
self.reg.set(*i, r);
}
}
Instruction::MeasureState => {
for i in 0..self.n {
let r = self.measure_sq(i);
self.reg.set(i, r);
}
}
Instruction::ConditionalInstruction {
indices,
condition,
inst,
} => {
let input = indices.iter().map(|x| self.reg[*x]).collect();
if condition(input) {
for ins in inst {
self.execute_instruction(&ins);
}
}
}
}
}
fn execute_circuit(&mut self, circuit: &Circuit, index: usize) {
if circuit.n > self.n {
panic!(
"number of qubits required for circuit is less than number of qubits in machine."
);
}
if index == circuit.ins.len() {
return;
}
self.execute_instruction(&circuit.ins[index]);
self.execute_circuit(circuit, index + 1);
}
pub fn get_reg_string(&self) -> String {
let mut s = String::new();
for i in 0..self.n {
s.push(if *self.reg.get(i).unwrap() { '1' } else { '0' });
}
s
}
fn get_counts(&self, circuit: &Circuit, n: usize) -> Counts {
let mut res = HashMap::new();
for _ in 0..n {
let mut s = self.clone();
s.execute_circuit(circuit, 0);
let label = s.get_reg_string();
if let Some(i) = res.get_mut(&label) {
*i += 1
} else {
res.insert(label, 0);
}
}
Counts {
data: res,
shots: n,
}
}
}
impl Simulator for Stabilizer {
fn execute_instruction(&mut self, ins: &Instruction) {
self.execute_instruction(ins);
}
fn execute_circuit(&mut self, circuit: &Circuit) {
self.execute_circuit(circuit, 0);
}
fn get_counts(&self, circuit: &Circuit, n: usize) -> Counts {
self.get_counts(circuit, n)
}
fn get_random(&mut self) -> &mut ThreadRng {
&mut self.rng
}
}