use std::collections::HashMap;
use crate::circuit::Circuit;
use crate::counts::Counts;
use crate::simulator_traits::Simulator;
use crate::sparse_mat::{kron, SMat};
use ndarray::{array, Array2};
use num::complex::Complex;
use num::Zero;
use rand::prelude::*;
use rayon::prelude::*;
use crate::instruction::{CliffordKind, Instruction, PauliKind};
use crate::sparse_quantum_state::SQS;
pub const MULTITHREADING_THRESHOLD: usize = 15;
#[derive(Clone, Debug)]
pub struct QuantumStateVector {
pub n: usize,
pub state: SQS,
pub reg: usize,
pub rng: ThreadRng,
}
pub fn kron_expt(u: &SMat, n: usize) -> SMat {
match n {
0 => SMat::from_dense(&array![[Complex::<f32>::from(1.0)]]),
1 => u.clone(),
_ => kron(&kron_expt(u, n - 1), u),
}
}
fn lift(t: &SMat, i: usize, n: usize) -> SMat {
let eye = SMat::eye(2);
let left = kron_expt(&eye, n - i - t.rows.ilog2() as usize);
let right = kron_expt(&eye, i);
kron(&left, &kron(t, &right))
}
fn ptt(p: &Vec<usize>) -> Vec<(usize, usize)> {
p.iter()
.enumerate()
.filter_map(|(mut src, &dest)| {
while src < dest {
src = p[src];
}
(src > dest).then_some((dest, src))
})
.collect()
}
fn tta(t: &Vec<(usize, usize)>) -> Vec<usize> {
t.iter()
.flat_map(|&(a, b)| {
if b.saturating_sub(a) <= 1 {
vec![a]
} else {
(a..b).chain((a..b - 1).rev()).collect()
}
})
.collect()
}
pub fn set_bit(x: usize, idx: usize, b: bool) -> usize {
let mask = !(1 << idx);
let flag = (b as usize) << idx;
x & mask | flag
}
impl QuantumStateVector {
pub fn new(n: usize) -> Self {
let rng = rand::thread_rng();
let state = SQS::new(2_u32.pow(n as u32));
state.set(0, Complex::new(1.0, 0.0));
Self {
n,
state,
reg: 0,
rng,
}
}
pub fn gate_apply_sq(&mut self, t: &Array2<Complex<f32>>, i: usize) {
self.state = self.state.dot(&lift(&SMat::from_dense(t), i, self.n));
}
pub fn gate_apply_nq(&mut self, t: &Array2<Complex<f32>>, i: &[usize]) {
let swap: Array2<Complex<f32>> = array![
[
Complex::<f32>::from(1.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0)
],
[
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(1.0),
Complex::<f32>::from(0.0)
],
[
Complex::<f32>::from(0.0),
Complex::<f32>::from(1.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0)
],
[
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(0.0),
Complex::<f32>::from(1.0)
]
];
let i_to_op = |i: &[usize]| {
i.iter().fold(kron_expt(&SMat::eye(2), self.n), |a, i| {
a.dot(&lift(&SMat::from_dense(&swap), *i, self.n))
})
};
let idk = i
.iter()
.copied()
.rev()
.chain(0..self.n)
.filter(|ind| !i.contains(ind))
.collect::<Vec<_>>();
let tr = tta(&ptt(&idk));
let tf = i_to_op(&tr);
let ft = i_to_op(&tr.into_iter().rev().collect::<Vec<_>>());
let all = tf.dot(&lift(&SMat::from_dense(&t), 0, self.n).dot(&ft));
self.state = self.state.dot(&all);
}
pub fn execute_instruction(&mut self, ins: &Instruction) {
match ins {
Instruction::Pauli { kind, target } => {
let ns = self.state.clone();
let op = |i| match kind {
PauliKind::X => {
let ndx = i ^ (1 << target);
ns.set(ndx, self.state.get(i));
}
PauliKind::Y => {
let phase_coef = if (i >> target) & 1 == 0 {
Complex::new(0.0, 1.0)
} else {
Complex::new(0.0, -1.0)
};
let ndx = i ^ (1 << target);
ns.set(ndx, self.state.get(i) * phase_coef);
}
PauliKind::Z => {
if (i >> target) & 1 == 1 {
ns.set(i, self.state.get(i) * Complex::new(-1.0, 0.0));
}
}
};
if self.n >= MULTITHREADING_THRESHOLD {
(0..self.state.n).par_bridge().for_each(|i| op(i));
} else {
(0..self.state.n).for_each(|i| op(i));
}
self.state = ns;
}
Instruction::Clifford { kind, target } => {
match kind {
CliffordKind::H => {
self.gate_apply_sq(kind.to_matrix(), target.to_single());
}
CliffordKind::S => {
self.gate_apply_sq(kind.to_matrix(), target.to_single());
}
CliffordKind::CNOT => {
let ns = self.state.clone();
let (target_qubit, control_qubit) = target.to_two();
for i in 0..self.state.n {
if (i >> control_qubit) & 1 == 1 {
let ndx = i ^ (1 << target_qubit);
ns.set(ndx, self.state.get(i));
}
}
self.state = ns;
}
}
}
Instruction::Gate { matrix, indices } => {
if indices.len() == 1 {
self.gate_apply_sq(matrix, indices[0]);
} else {
self.gate_apply_nq(matrix, indices);
}
}
Instruction::Measure { indices } => self.measure_qubits(indices, 0.0),
Instruction::MeasureNoisy { indices, p } => {
self.measure_qubits(indices, *p);
}
Instruction::MeasureState => self.measure(),
Instruction::Identity => {}
Instruction::ConditionalInstruction {
indices,
condition,
inst,
} => {
let input = indices
.iter()
.map(|x| ((self.reg >> *x) & 1) != 0)
.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 measure(&mut self) {
let mut rand = self.rng.gen_range(0.0..1.0_f32);
let mut k: i32 = -1;
for (i, theta) in self.state.iter().enumerate() {
rand -= theta.norm_sqr();
if rand < 0.0 {
k = i as i32;
break;
}
}
if k == -1 {
k = self.n as i32 - 1;
}
self.state.fill_zero();
self.state.set(0, 1.0.into());
self.reg = k as usize;
}
pub fn measure_qubits(&mut self, indices: &[usize], p: f32) {
let norm = self.state.iter().map(|x| x.norm_sqr()).sum::<f32>();
for qubit in indices {
let p_zero;
p_zero = self
.state
.iter()
.enumerate()
.par_bridge()
.map(|(idx, amplitude)| {
if (idx >> qubit) & 1 == 0 {
amplitude.norm_sqr()
} else {
0.0
}
})
.sum::<f32>()
/ norm;
let rf = self.rng.gen_range(0.0..1.0_f32);
let mut outcome = if rf < p_zero { 0 } else { 1 };
if p != 0.0 && self.rng.gen_range(0.0..1.0) < p {
outcome = !outcome;
}
self.reg = set_bit(self.reg, *qubit, rf > p_zero);
let prob = if outcome == 0 { p_zero } else { 1.0 - p_zero };
let scale = Complex::from(prob.sqrt().recip());
(0..self.n as u32).into_par_iter().for_each(|i| {
let bit = (i >> qubit) & 1;
let k = self.state.get(i);
self.state.set(
i,
if bit == outcome {
k * scale
} else {
Complex::zero()
},
);
});
}
}
fn get_counts(&self, circuit: &Circuit, n: usize) -> Counts {
let mut res = HashMap::new();
let nq = self.n;
for _ in 0..n {
let mut s = self.clone();
s.execute_circuit(circuit, 0);
let label = format!("{:0nq$b}", s.reg);
if let Some(i) = res.get_mut(&label) {
*i += 1
} else {
res.insert(label, 0);
}
}
Counts {
data: res,
shots: n,
}
}
}
impl Simulator for QuantumStateVector {
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
}
}