poulpy/poulpy-hal/src/layouts/stats.rs
poulpy-hal/src/layouts/stats.rsRaw
1use rug::{
2 Float,
3 float::Round,
4 ops::{AddAssignRound, DivAssignRound, SubAssignRound},
5};
6
7use crate::layouts::{Backend, DataRef, VecZnx, VecZnxBig, VecZnxBigToRef, ZnxInfos};
8
9impl<D: DataRef> VecZnx<D> {
10 pub fn std(&self, base2k: usize, col: usize) -> f64 {
11 let prec: u32 = (self.size() * base2k) as u32;
12 let mut data: Vec<Float> = (0..self.n()).map(|_| Float::with_val(prec, 0)).collect();
13 self.decode_vec_float(base2k, col, &mut data);
14 // std = sqrt(sum((xi - avg)^2) / n)
15 let mut avg: Float = Float::with_val(prec, 0);
16 data.iter().for_each(|x| {
17 avg.add_assign_round(x, Round::Nearest);
18 });
19 avg.div_assign_round(Float::with_val(prec, data.len()), Round::Nearest);
20 data.iter_mut().for_each(|x| {
21 x.sub_assign_round(&avg, Round::Nearest);
22 });
23 let mut std: Float = Float::with_val(prec, 0);
24 data.iter().for_each(|x| std += x * x);
25 std.div_assign_round(Float::with_val(prec, data.len()), Round::Nearest);
26 std = std.sqrt();
27 std.to_f64()
28 }
29}
30
31impl<D: DataRef, B: Backend + Backend<ScalarBig = i64>> VecZnxBig<D, B> {
32 pub fn std(&self, base2k: usize, col: usize) -> f64 {
33 let self_ref: VecZnxBig<&[u8], B> = self.to_ref();
34 let znx: VecZnx<&[u8]> = VecZnx {
35 data: self_ref.data,
36 n: self_ref.n,
37 cols: self_ref.cols,
38 size: self_ref.size,
39 max_size: self_ref.max_size,
40 };
41 znx.std(base2k, col)
42 }
43}