Skip to main content

poulpy_hal/layouts/
stats.rs

1use dashu_float::{FBig, round::mode::HalfEven};
2
3use crate::layouts::{Backend, HostDataRef, VecZnx, VecZnxBig};
4
5/// Summary statistics (max absolute value and standard deviation) of a
6/// polynomial vector's decoded floating-point coefficients.
7pub struct Stats {
8    max: f64,
9    std: f64,
10}
11
12impl Stats {
13    /// Returns the maximum absolute coefficient value.
14    pub fn max(&self) -> f64 {
15        self.max
16    }
17
18    /// Returns the standard deviation of the coefficients.
19    pub fn std(&self) -> f64 {
20        self.std
21    }
22}
23
24impl<D: HostDataRef> VecZnx<D> {
25    /// Computes [`Stats`] (max absolute value and standard deviation) for
26    /// column `col` by decoding all limbs into arbitrary-precision floats.
27    pub fn stats(&self, base2k: usize, col: usize) -> Stats {
28        let mut data: Vec<FBig<HalfEven>> = (0..self.n()).map(|_| FBig::ZERO).collect();
29        self.decode_vec_float(base2k, col, &mut data);
30
31        // std = sqrt(sum((xi - avg)^2) / n)
32        let mut avg: FBig<HalfEven> = FBig::ZERO;
33        let mut max: FBig<HalfEven> = FBig::ZERO;
34
35        data.iter().for_each(|x| {
36            avg = avg.clone() + x.clone();
37            let abs_x = if x < &FBig::<HalfEven>::ZERO { -x.clone() } else { x.clone() };
38            if abs_x > max {
39                max = abs_x;
40            }
41        });
42        avg /= FBig::from(data.len() as i64);
43        data.iter_mut().for_each(|x| {
44            *x = x.clone() - avg.clone();
45        });
46        let mut variance: FBig<HalfEven> = FBig::ZERO;
47        data.iter().for_each(|x| {
48            variance = variance.clone() + x.clone() * x.clone();
49        });
50        variance /= FBig::from(data.len() as i64);
51
52        // Final output is f64; to_f64() always succeeds (returns nearest f64).
53        // f64::try_from(FBig) fails with LossOfPrecision for nearly all values.
54        Stats {
55            std: variance.to_f64().value().sqrt(),
56            max: max.to_f64().value(),
57        }
58    }
59}
60
61impl<D: HostDataRef, B: Backend + Backend<ScalarBig = i64>> VecZnxBig<D, B> {
62    pub fn stats(&self, base2k: usize, col: usize) -> Stats {
63        let shape = self.shape();
64        let znx: VecZnx<&[u8]> =
65            VecZnx::from_data_with_max_size(self.data.as_ref(), shape.n(), shape.cols(), shape.size(), shape.max_size());
66        znx.stats(base2k, col)
67    }
68}