bcore/process/
mod.rs

1use crate::api::Estimator;
2use crate::error::ApiError;
3use crate::Weight;
4use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis};
5
6pub fn spatial_average_concentration(
7    concentration_record: &ArrayView2<f64>,
8    full_volume: &ArrayView2<f64>,
9) -> Array1<f64> {
10    let weighted_concentration = concentration_record * full_volume;
11    let numerator = weighted_concentration.sum_axis(Axis(1));
12    let denominator = full_volume.sum_axis(Axis(1));
13    numerator / denominator
14}
15
16pub fn normalize_concentration(
17    concentration_record: &ArrayView2<f64>,
18    mean_concentration: &ArrayView1<f64>,
19) -> Array2<f64> {
20    concentration_record / mean_concentration
21}
22
23pub fn variance_concentration(
24    concentration_record: &ArrayView2<f64>,
25    full_volume: &ArrayView2<f64>,
26) -> Array1<f64> {
27    let mean_concentration = spatial_average_concentration(concentration_record, full_volume).insert_axis(Axis(1));    
28    ((concentration_record - mean_concentration).powi(2) * full_volume).sum_axis(Axis(1))
29}
30
31
32pub fn variance_concentration_from_mean(
33    concentration_record: &ArrayView2<f64>,
34    mean_concentration: &ArrayView1<f64>,
35    full_volume: &ArrayView2<f64>,
36) -> Array1<f64> {
37    ((concentration_record - mean_concentration).pow2() * full_volume).sum_axis(Axis(1))
38}
39
40// def normalize_concentration(
41//     raw_concentration: np.ndarray, volumes: np.ndarray
42// ) -> Tuple[np.ndarray, float, float]:
43//     vtot = np.sum(volumes, axis=1)
44//     mean_concentration = np.sum(raw_concentration * volumes, axis=1) / vtot
45//     mean_concentration = mean_concentration.reshape(-1, 1)
46//     variance = (
47//         np.sum(np.power(raw_concentration - mean_concentration, 2) * volumes, axis=1)
48//         / vtot
49//     )
50//     return raw_concentration / mean_concentration, mean_concentration, variance
51
52pub struct Histogram {
53    bins: Vec<f64>,
54    counts: Vec<f64>,
55    bin_counts: usize,
56}
57
58impl Histogram {
59    pub fn new(bin_counts: usize) -> Self {
60        Histogram {
61            bins: Vec::new(),
62            counts: Vec::new(),
63            bin_counts,
64        }
65    }
66
67    pub fn add(&mut self, values: Vec<f64>) {
68        if values.is_empty() {
69            return;
70        }
71
72        // Determine the range of the data
73        let min = values.iter().cloned().fold(f64::INFINITY, f64::min);
74        let max = values.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
75        let bin_width = (max - min) / self.bin_counts as f64;
76
77        // Initialize bins if not already done
78        if self.bins.is_empty() {
79            self.bins = (0..self.bin_counts)
80                .map(|i| min + i as f64 * bin_width)
81                .collect();
82            self.counts = vec![0.0; self.bin_counts];
83        }
84
85        // Add values to the appropriate bins
86        for value in values {
87            let bin_index = ((value - min) / bin_width).floor() as usize;
88            let bin_index = bin_index.min(self.bin_counts - 1); // Clamp to the last bin if necessary
89            self.counts[bin_index] += 1.0;
90        }
91    }
92
93    pub fn get_bins(&self) -> &[f64] {
94        &self.bins
95    }
96
97    pub fn get_counts(&self) -> &[f64] {
98        &self.counts
99    }
100}
101
102pub(crate) fn estimate(
103    etype: Estimator,
104    weight: &Weight,
105    rx: &Array1<f64>,
106) -> Result<f64, ApiError> {
107    let weighted_estimator = match weight {
108        Weight::Single(sw) => (rx * *sw).sum(),
109        Weight::Multiple(mw) => {
110            if mw.len() != rx.len() {
111                return Err(ApiError::ShapeError);
112            }
113            rx.iter().zip(mw).map(|(x, w)| x * w).sum()
114        }
115    };
116
117    if weighted_estimator == 0. {
118        return Ok(0.);
119    }
120    match etype {
121        // Estimator::MonteCarlo => self
122        //     .get_properties(key, i_export)
123        //     .map(|x| x.sum() / (x.len() as f64))
124        //     .unwrap_or(0.),
125        Estimator::MonteCarlo => {
126            let denum = match weight {
127                Weight::Single(sw) => (rx.dim() as f64) * *sw,
128                Weight::Multiple(mw) => mw.iter().sum(),
129            };
130
131            Ok(weighted_estimator / denum) //Normalise
132        }
133
134        Estimator::Weighted => Ok(weighted_estimator),
135    }
136}
bcore/process/mod.rs

bcore/process/
mod.rs
