1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

use{
    crate::*,
    std::{convert::*, f64::consts::LOG10_E}
};

#[derive(Debug)]
/// Possible errors that can occur during gluing together
/// WangLandau intervals or Entropic Sampling intervals
pub enum GlueErrors{
    /// `original_hist.borders_clone()` failed
    BorderCreation(HistErrors),
    /// Nothing to be glued, glue interval list was empty
    EmptyList,
    /// Binary search failed - PartialOrd::partial_cmp returned None
    BinarySearch,
    /// # Glue interval and intervals to be glued do not match
    /// * Likely `original_hist` is to small
    OutOfBounds,
    /// The intervals need to overlap, otherwise no gluing can occur
    NoOverlap,
}

impl From<HistErrors> for GlueErrors{
    fn from(e: HistErrors) -> Self {
        GlueErrors::BorderCreation(e)
    }
}

/// # Normalize log10 probability density
/// * input: Slice containing log10 of (non normalized) probability density
/// * afterwards, it will be normalized, i.e., sum_i 10^log10_density\[i\] ≈ 1
pub fn norm_log10_sum_to_1(log10_density: &mut[f64]){
    // prevent errors due to small or very large numbers
    subtract_max(log10_density);

    // calculate actual sum in non log space
    let sum = log10_density.iter()
        .fold(0.0, |acc, &val| {
            if val.is_finite(){
               acc +  10_f64.powf(val)
            } else {
                acc
            }
        }  
    );
    
    let sum = sum.log10();
    log10_density.iter_mut()
        .for_each(|val| *val -= sum);
}



pub(crate) fn height_correction(log10_vec: &mut [Vec<f64>], z_vec: &[f64]){
    log10_vec.iter_mut()
        .skip(1)
        .zip(z_vec.iter())
        .for_each( |(vec, &z)|
            vec.iter_mut()
                .for_each(
                    |val| 
                    {
                        *val += z;
                    }
                )
        );
}





/// subtracts maximum, if it is finite
pub(crate) fn subtract_max(list: &mut[f64]) -> f64
{
    let max = list
        .iter()
        .copied()
        .fold(f64::NAN, f64::max);

    if max.is_finite() {
        list.iter_mut()
            .for_each(|val| *val -= max);
    }
    max
}

pub(crate) fn ln_to_log10(slice: &mut [f64])
{
    slice.iter_mut()
            .for_each(|val| *val *= LOG10_E);
}

pub(crate) fn log10_to_ln(slice: &mut [f64])
{
    slice.iter_mut()
        .for_each(|val| *val /= LOG10_E);
}