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

// use std::f64;
//

// use ordered_float::NotNaN;
// use bio::stats::{LogProb, Prob};
// use bio::stats::combinatorics::combinations;
//
// use model::{Variant, DiscreteAlleleFreqs, AlleleFreq};
// use priors::{Model, InfiniteSitesNeutralVariationModel, PairModel};

use itertools::Itertools;

use model::{AlleleFreq, DiscreteAlleleFreqs};

pub fn allele_freqs(ploidy: u32) -> DiscreteAlleleFreqs {
    (0..ploidy + 1).map(|m| AlleleFreq(m as f64 / ploidy as f64)).collect_vec()
}


/*
pub struct NormalNormalModel {
    inner: InfiniteSitesNeutralVariationModel,
    ploidy: u32,
    allele_freqs: DiscreteAlleleFreqs
}


impl NormalNormalModel {
    pub fn new(ploidy: u32, heterozygosity: Prob) -> Self {
        NormalNormalModel {
            inner: InfiniteSitesNeutralVariationModel::new(2, ploidy, heterozygosity),
            ploidy: ploidy,
            allele_freqs: allele_freqs(ploidy)
        }
    }
}


impl PairModel<DiscreteAlleleFreqs, DiscreteAlleleFreqs> for NormalNormalModel {
    fn prior_prob(&self, af_first: AlleleFreq, af_second: AlleleFreq, _: &Variant) -> LogProb {
        let sample_prior = |af: AlleleFreq| {
            let m = *af * self.ploidy as f64;
            let valid_genotypes = combinations(self.ploidy, m);
            let all_genotypes = combinations()
        };
    }

    fn joint_prob<L, O>(
        &self,
        af_first: &DiscreteAlleleFreqs,
        af_second: &DiscreteAlleleFreqs,
        likelihood_first: &L,
        likelihood_second: &O,
        _: &Variant
    ) -> LogProb where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        let p_second = LogProb::ln_sum_exp(&af_second.iter().map(|&af_second| {
            likelihood_second(af_second, AlleleFreq(0.0))
        }).collect_vec());

        let prob = LogProb::ln_sum_exp(&af_first.iter().map(|&af_first| {
            let p_first = likelihood_first(af_first, AlleleFreq(0.0));
            let prob = p_first + p_second;

            prob
        }).collect_vec());

        prob
    }

    fn marginal_prob<L, O>(
        &self,
        likelihood_first: &L,
        likelihood_second: &O,
        variant: &Variant
    ) -> LogProb where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        let p = self.joint_prob(
            self.allele_freqs().0,
            self.allele_freqs().1,
            likelihood_first,
            likelihood_second,
            variant
        );

        p
    }

    fn map<L, O>(
        &self,
        likelihood_first: &L,
        likelihood_second: &O,
        _: &Variant
    ) -> (AlleleFreq, AlleleFreq) where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        fn calc_map<L: Fn(AlleleFreq, AlleleFreq) -> LogProb>(likelihood: &L, afs: &DiscreteAlleleFreqs) -> AlleleFreq {
            let (_, map) = afs.iter().minmax_by_key(|&af| {
                let p = likelihood(*af, AlleleFreq(0.0));
                NotNaN::new(*p).expect("probability is NaN")
            }).into_option().expect("prior has empty allele frequency spectrum");
            *map
        }

        let map_first = calc_map(likelihood_first, self.allele_freqs().0);
        let map_second = calc_map(likelihood_second, self.allele_freqs().1);

        (map_first, map_second)
    }

    fn allele_freqs(&self) -> (&DiscreteAlleleFreqs, &DiscreteAlleleFreqs) {
        (&self.allele_freqs, &self.allele_freqs)
    }
}*/