demes-forward 0.6.0

rust API for iterating demes models forwards in time
Documentation 
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
125
126
127
128
129
130
131
132
133
134
135
136
137
138

pub fn update_ancestry_proportions(
    sources: &[usize],
    source_proportions: &[f64],
    ancestry_proportions: &mut [f64],
) {
    assert_eq!(sources.len(), source_proportions.len());
    let sum = source_proportions.iter().fold(0.0, |a, b| a + b);
    ancestry_proportions.iter_mut().for_each(|a| *a *= 1. - sum);
    sources
        .iter()
        .zip(source_proportions.iter())
        .for_each(|(source, proportion)| ancestry_proportions[*source] += proportion);
}

pub fn ancestry_proportions_from_graph(
    graph: &crate::ForwardGraph,
    child_deme: usize,
) -> Option<Vec<f64>> {
    graph.offspring_deme_sizes()?;

    let mut rv = vec![0.0; graph.offspring_deme_sizes().unwrap().len()];

    let deme = graph.demes_graph().get_deme(child_deme).unwrap();
    let bwtime = graph
        .time_to_backward(graph.last_time_updated().unwrap())
        .unwrap()
        .unwrap();
    if bwtime > deme.start_time() || bwtime < deme.end_time() {
        return Some(rv);
    }

    if !deme.ancestor_indexes().is_empty() && bwtime == deme.start_time() {
        for (a, p) in deme
            .ancestor_indexes()
            .iter()
            .zip(deme.proportions().iter())
        {
            rv[*a] = f64::from(*p);
        }
    } else {
        rv[child_deme] = 1.0;
    }

    let mut sources: Vec<usize> = vec![];
    let mut source_proportions: Vec<f64> = vec![];

    let bwtime: f64 = graph
        .time_to_backward(graph.last_time_updated().unwrap())
        .unwrap()
        .unwrap()
        .into();
    // NOTE: we subract 1 because:
    // * "last time updated" refers to the birth time of the current parents.
    // * this fn is trying to get the ancestry proportions of the next generation of children.
    // * the children are 1 generation closer to the present, so we subtract 1!
    let bwtime = demes::Time::try_from(bwtime - 1.0).unwrap();
    for p in graph
        .demes_graph()
        .pulses()
        .iter()
        .filter(|p| p.time() == f64::from(bwtime))
    {
        sources.clear();
        source_proportions.clear();
        let dest = graph
            .demes_graph()
            .demes()
            .iter()
            .position(|d| d.name() == p.dest())
            .unwrap();
        if dest == child_deme {
            for (s, d) in p.sources().iter().zip(p.proportions().iter()) {
                let source = graph
                    .demes_graph()
                    .demes()
                    .iter()
                    .position(|d| d.name() == s)
                    .unwrap();
                sources.push(source);
                source_proportions.push(f64::from(*d));
            }
            update_ancestry_proportions(&sources, &source_proportions, &mut rv);
        }
    }

    sources.clear();
    source_proportions.clear();

    for m in graph
        .demes_graph()
        .migrations()
        .iter()
        .filter(|m| bwtime >= m.end_time() && bwtime < m.start_time())
    {
        let d = graph
            .demes_graph()
            .demes()
            .iter()
            .position(|deme| deme.name() == m.dest())
            .unwrap();
        if d == child_deme {
            let s = graph
                .demes_graph()
                .demes()
                .iter()
                .position(|deme| deme.name() == m.source())
                .unwrap();
            sources.push(s);
            source_proportions.push(f64::from(m.rate()));
        }
    }
    update_ancestry_proportions(&sources, &source_proportions, &mut rv);

    Some(rv)
}

pub fn test_model_duration(graph: &mut crate::ForwardGraph) {
    for time in graph.time_iterator() {
        graph.update_state(time).unwrap();
        // assert!(graph.parental_demes().is_some(), "{}", time);
        assert!(
            graph.num_extant_parental_demes() > 0,
            "{:?} {:?}",
            time,
            graph.end_time()
        );
        assert!(graph
            .parental_deme_sizes()
            .unwrap()
            .iter()
            .any(|size| size > &0.0));
        if time == graph.end_time() - 1.0.into() {
            assert!(graph.offspring_deme_sizes().is_none(), "time = {time:?}");
        } else {
            assert!(graph.offspring_deme_sizes().is_some(), "time = {time:?}");
        }
    }
}