mod.rs 3.13 KB
Newer Older
afischer@lsv.uni-saarland.de's avatar
lib  
afischer@lsv.uni-saarland.de committed
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 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124

//// STIRLING APPROXIMATIONS
pub mod stirling_approximation{
    use std::f64;
    pub fn fac_stirling(n: f64) -> f64{
        let result: f64 = (2.0*f64::consts::PI*n).sqrt();
        let frac: f64 = (n/f64::consts::E).powf(n);
        result * frac / f64::consts::LN_2
    }


    pub fn log_fac_stirling(n: f64) -> f64{
        let result: f64 = n*n.ln() - n;
        result / f64::consts::LN_2
    }


    pub fn log_binomial_coefficient_stirling(n: usize, k: usize) -> f64{
        let res: f64;
        if n == k{
            res = 0.0;
        } else {
            if (n == 0) || (k == 0) {
                res = f64::INFINITY;
            } else {
                let n = n as f64;
                let k = k as f64;
                res = log_fac_stirling(n) - log_fac_stirling(k) - log_fac_stirling(n - k);
            }
        }

        res
    }


    pub fn gammaln_stirling(z: f64) -> f64{
        let result: f64 = z*z.ln() - z + 0.5*(2.0*f64::consts::PI/z).ln();
        result / f64::consts::LN_2
    }
}

//// DATA-TO-MODEL FOR DISTRIBUTIONS
pub mod d2m_codes{
    use basic_codes::stirling_approximation;
    pub fn d2m_nonzero_distribution_cost(counts: usize, bins: usize) -> f64{
        //println!("approximating d2m NONZERO distribution cost: {} counts in {} bins", counts, bins);
        if counts < bins{
            panic!("Number of counts in non-zero distribution must be >= number of bins!");
        }
        else if bins == 1
        {
            return 0.0;
        }
        stirling_approximation::log_binomial_coefficient_stirling(counts - 1, bins - 1)
    }

    pub fn d2m_distribution_cost(counts: usize, bins: usize) -> f64{
        //println!("approximating d2m distribution cost: {} counts in {} bins", counts, bins);
        if bins == 1
        {
            return 0.0;
        }
        stirling_approximation::log_binomial_coefficient_stirling(counts + bins - 1, bins - 1)

    }
}


//// uic

pub mod uic{
    use std;
    use atomic_f64::AtomicF64;

    const MAX_UIC_MEMOIZATION: usize = 1024*1024;

    lazy_static!{
        static ref UIC : [AtomicF64; 1024] = {
            let mut result : [AtomicF64;1024] = unsafe { std::mem::uninitialized() };
            for i in 0..MAX_UIC_MEMOIZATION {
                let uninit = std::mem::replace(result.get_mut(i).unwrap(), AtomicF64::new(0.0));
                std::mem::forget(uninit)
            };
            result
        };
    }

    use std::f64;
    const UIC_NORMALIZATION_CONSTANT: f64 = 1.5185673663648485;

    fn uic_proper(n:usize) -> f64
    {
        let mut res: f64 = 0.0;
        let mut n = n as f64;
        n = n.log(2.0);
        while n > 0.0{
            res = res + n;
            n = n.log(2.0);
        }

        res + UIC_NORMALIZATION_CONSTANT
    }

    pub fn uic_memoized(n: usize) -> f64 {
        if n < MAX_UIC_MEMOIZATION
        {
            let mut res: f64 = UIC[n-1].get();
            if res != 0.0
            {
                res
            }
            else
            {
                res = uic_proper(n);
                UIC[n-1].set(res);
                res
            }
        }
        else
        {
            uic_proper(n)
        }
    }
}