use scirs2_core::ndarray::{Array1, Array2, Array4};
use scirs2_core::Complex64;
use super::super::core::SciRS2ProcessTomographer;
use super::super::results::{ExperimentalData, ReconstructionQuality};
use super::utils::calculate_reconstruction_quality;
use crate::DeviceResult;
pub fn reconstruct_bayesian(
tomographer: &SciRS2ProcessTomographer,
experimental_data: &ExperimentalData,
) -> DeviceResult<(Array4<Complex64>, ReconstructionQuality)> {
let num_qubits = (experimental_data.input_states[0].dim().0 as f64).log2() as usize;
let dim = 1 << num_qubits;
let mut process_matrix = Array4::zeros((dim, dim, dim, dim));
for i in 0..dim {
process_matrix[[i, i, i, i]] = Complex64::new(1.0, 0.0);
}
let log_likelihood = calculate_bayesian_log_likelihood(&process_matrix, experimental_data)?;
let reconstruction_quality =
calculate_reconstruction_quality(&process_matrix, experimental_data, log_likelihood);
Ok((process_matrix, reconstruction_quality))
}
fn calculate_bayesian_log_likelihood(
process_matrix: &Array4<Complex64>,
experimental_data: &ExperimentalData,
) -> DeviceResult<f64> {
let mut log_likelihood = 0.0;
for (observed, &uncertainty) in experimental_data
.measurement_results
.iter()
.zip(experimental_data.measurement_uncertainties.iter())
{
let predicted = 0.5; let diff = observed - predicted;
let variance = uncertainty * uncertainty;
log_likelihood -= 0.5 * (diff * diff / variance);
}
let dim = process_matrix.dim().0;
for i in 0..dim {
for j in 0..dim {
for k in 0..dim {
for l in 0..dim {
if i == j && k == l && i == k {
let deviation = (process_matrix[[i, j, k, l]].re - 1.0).abs();
log_likelihood -= 0.5 * deviation * deviation;
} else {
let magnitude = process_matrix[[i, j, k, l]].norm();
log_likelihood -= 0.5 * magnitude * magnitude;
}
}
}
}
}
Ok(log_likelihood)
}