quantrs2_device/process_tomography/reconstruction/
ensemble.rs

1//! Ensemble methods reconstruction
2
3use scirs2_core::ndarray::{Array1, Array2, Array4};
4use scirs2_core::Complex64;
5
6use super::super::core::SciRS2ProcessTomographer;
7use super::super::results::{ExperimentalData, ReconstructionQuality};
8use super::utils::calculate_reconstruction_quality;
9use super::{bayesian, compressed_sensing, linear_inversion, maximum_likelihood};
10use crate::DeviceResult;
11
12/// Ensemble reconstruction combining multiple methods
13pub fn reconstruct_ensemble(
14    tomographer: &SciRS2ProcessTomographer,
15    experimental_data: &ExperimentalData,
16) -> DeviceResult<(Array4<Complex64>, ReconstructionQuality)> {
17    // Run multiple reconstruction methods
18    let (li_process, li_quality) =
19        linear_inversion::reconstruct_linear_inversion(tomographer, experimental_data)?;
20    let (mle_process, mle_quality) =
21        maximum_likelihood::reconstruct_maximum_likelihood(tomographer, experimental_data)?;
22    let (cs_process, cs_quality) =
23        compressed_sensing::reconstruct_compressed_sensing(tomographer, experimental_data)?;
24    let (bayes_process, bayes_quality) =
25        bayesian::reconstruct_bayesian(tomographer, experimental_data)?;
26
27    // Calculate weights based on reconstruction quality
28    let li_weight = calculate_weight(&li_quality);
29    let mle_weight = calculate_weight(&mle_quality);
30    let cs_weight = calculate_weight(&cs_quality);
31    let bayes_weight = calculate_weight(&bayes_quality);
32
33    let total_weight = li_weight + mle_weight + cs_weight + bayes_weight;
34
35    // Weighted combination
36    let dim = li_process.dim().0;
37    let mut ensemble_process = Array4::zeros((dim, dim, dim, dim));
38
39    for i in 0..dim {
40        for j in 0..dim {
41            for k in 0..dim {
42                for l in 0..dim {
43                    ensemble_process[[i, j, k, l]] = (li_process[[i, j, k, l]] * li_weight
44                        + mle_process[[i, j, k, l]] * mle_weight
45                        + cs_process[[i, j, k, l]] * cs_weight
46                        + bayes_process[[i, j, k, l]] * bayes_weight)
47                        / total_weight;
48                }
49            }
50        }
51    }
52
53    let log_likelihood = calculate_ensemble_log_likelihood(&ensemble_process, experimental_data)?;
54    let reconstruction_quality =
55        calculate_reconstruction_quality(&ensemble_process, experimental_data, log_likelihood);
56
57    Ok((ensemble_process, reconstruction_quality))
58}
59
60/// Calculate weight for ensemble based on reconstruction quality
61fn calculate_weight(quality: &ReconstructionQuality) -> f64 {
62    // Higher log-likelihood and better physical validity get higher weights
63    let ll_weight = (quality.log_likelihood + 100.0).max(0.1); // Offset to ensure positive
64    let physical_weight = quality.physical_validity.positivity_measure
65        * quality.physical_validity.trace_preservation_measure;
66    let stability_weight = 1.0 / (1.0 + quality.condition_number / 100.0);
67
68    ll_weight * physical_weight * stability_weight
69}
70
71/// Calculate ensemble log-likelihood
72fn calculate_ensemble_log_likelihood(
73    process_matrix: &Array4<Complex64>,
74    experimental_data: &ExperimentalData,
75) -> DeviceResult<f64> {
76    let mut log_likelihood = 0.0;
77
78    for (observed, &uncertainty) in experimental_data
79        .measurement_results
80        .iter()
81        .zip(experimental_data.measurement_uncertainties.iter())
82    {
83        let predicted = 0.5; // Placeholder
84        let diff = observed - predicted;
85        let variance = uncertainty * uncertainty;
86        log_likelihood -= 0.5 * (diff * diff / variance);
87    }
88
89    Ok(log_likelihood)
90}