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
extern crate phyloprob; extern crate phylofold; pub use phyloprob::*; pub use phylofold::*; pub use phylofold::{RnaId, RnaIdPair, Prob4dMat, SparseProbMat}; pub type Prob4dMatsWithRnaIdPairs = FxHashMap<RnaIdPair, Prob4dMat>; pub type ProbsWithRnaIds = Vec<Probs>; pub type ProbsWithPosPairs = FxHashMap<PosPair, Prob>; pub type ProbMatsWithRnaIds = Vec<SparseProbMat>; pub type Col = Vec<Char>; pub type Cols = Vec<Col>; pub type PosMaps = Vec<Pos>; pub type PosMapSets = Vec<PosMaps>; #[derive(Debug)] pub struct SeqAlign { pub cols: Cols, pub pos_map_sets: PosMapSets, } pub struct MeaCss { pub bpa_pos_pairs: PosPairs, pub ea: Mea, } pub type SeqId = String; pub type SeqIds = Vec<SeqId>; impl SeqAlign { pub fn new() -> SeqAlign { SeqAlign { cols: Cols::new(), pos_map_sets: PosMapSets::new(), } } } impl MeaCss { pub fn new() -> MeaCss { MeaCss { bpa_pos_pairs: PosPairs::new(), ea: 0., } } } pub const GAP: Char = '-' as Char; pub fn phyloalifold(mix_bpp_mat: &ProbMat, mix_upp_mat: &Probs, gamma: Prob, sa: &SeqAlign) -> MeaCss { let sa_len = sa.cols.len(); let mut mea_mat = vec![vec![0.; sa_len]; sa_len]; let sa_len = sa_len as Pos; for sub_sa_len in 1 .. sa_len + 1 { for i in 0 .. sa_len + 1 - sub_sa_len { let j = i + sub_sa_len - 1; let (long_i, long_j) = (i as usize, j as usize); if i == j { mea_mat[long_i][long_j] = mix_upp_mat[long_i]; continue; } let mut mea = mea_mat[long_i + 1][long_j] + mix_upp_mat[long_i]; let ea = mea_mat[long_i][long_j - 1] + mix_upp_mat[long_j]; if ea > mea { mea = ea; } let ea = mea_mat[long_i + 1][long_j - 1] + gamma * mix_bpp_mat[long_i][long_j]; if ea > mea { mea = ea; } for k in long_i .. long_j { let ea = mea_mat[long_i][k] + mea_mat[k + 1][long_j]; if ea > mea { mea = ea; } } mea_mat[long_i][long_j] = mea; } } let mut mea_css = MeaCss::new(); let mut pos_pair_stack = vec![(0, sa_len - 1)]; while pos_pair_stack.len() > 0 { let pos_pair = pos_pair_stack.pop().expect("Failed to pop an element of a vector."); let (i, j) = pos_pair; if j <= i {continue;} let (long_i, long_j) = (i as usize, j as usize); let mea = mea_mat[long_i][long_j]; if mea == mea_mat[long_i + 1][long_j] + mix_upp_mat[long_i] { pos_pair_stack.push((i + 1, j)); } else if mea == mea_mat[long_i][long_j - 1] + mix_upp_mat[long_j] { pos_pair_stack.push((i, j - 1)); } else if mea == mea_mat[long_i + 1][long_j - 1] + gamma * mix_bpp_mat[long_i][long_j] { pos_pair_stack.push((i + 1, j - 1)); mea_css.bpa_pos_pairs.push(pos_pair); } else { for k in i .. j { let long_k = k as usize; if mea == mea_mat[long_i][long_k] + mea_mat[long_k + 1][long_j] { pos_pair_stack.push((i, k)); pos_pair_stack.push((k + 1, j)); break; } } } } mea_css.ea = mea_mat[0][sa_len as usize - 1]; mea_css }