kmers/
dot.rs

1use crate::Kmer;
2
3/// Compute the dot product between two k-mer frequency spectra.
4pub fn dot<Lhs, Rhs>(mut lhs: Lhs, mut rhs: Rhs) -> f64
5where
6    Lhs: Iterator<Item = (Kmer, usize)>,
7    Rhs: Iterator<Item = (Kmer, usize)>,
8{
9    let mut d = 0.0;
10
11    let mut lhs_next = lhs.next();
12    let mut rhs_next = rhs.next();
13    while let (Some(x), Some(y)) = (&lhs_next, &rhs_next) {
14        if x.0 < y.0 {
15            lhs_next = lhs.next();
16        } else if x.0 > y.0 {
17            rhs_next = rhs.next();
18        } else {
19            d += (x.1 * y.1) as f64;
20            lhs_next = lhs.next();
21            rhs_next = rhs.next();
22        }
23    }
24
25    d
26}
27
28#[cfg(test)]
29mod test {
30    use crate::frequency::{unique, frequency_vector_iter};
31
32    use super::*;
33
34    #[test]
35    fn test_1() {
36        let k = 11;
37
38        let lhs_seq = "CCACGACACATCATAGTCAAATTTCTAAAAATTAAGGACCACACAAACACCTTAAAAACAGAGAAAGAAAAATGACACATTACCTATACTGCAAAAACAATTCAAATGACAAAGTATTATTCATTGGAAACCACAGAGTTCAGAATGAAG";
39        let mut lhs_kmers = Kmer::make_many(k, lhs_seq);
40        lhs_kmers.sort();
41        lhs_kmers.dedup();
42
43        let rhs_seq = "CCACGACACATCACAGTCAAATTTCTAAAAATTAAGGACCACACAAACGCCTTAAAAACAGAGAAAGAAAAATGACATATTACCTATACTGCAAAAACAATTCAAATGACAAAGTATTATTCATTGGCAACCACAGAGTTCAGAATGAAG";
44        let mut rhs_kmers = Kmer::make_many(k, rhs_seq);
45        rhs_kmers.sort();
46        rhs_kmers.dedup();
47
48        let lhs_kmer_frequencies = unique(lhs_kmers.iter().map(|x| x.clone()));
49        let rhs_kmer_frequencies = unique(rhs_kmers.iter().map(|x| x.clone()));
50        let d = dot(lhs_kmer_frequencies, rhs_kmer_frequencies);
51        assert_eq!(d, 96.0);
52    }
53
54    #[test]
55    fn test_2() {
56        let k = 5;
57
58        let lhs_seq = "CCACGACACATCATAGTCAAATTTCTAAAAATTAAGGACCACACAAACACCTTAAAAACAGAGAAAGAAAAATGACACATTACCTATACTGCAAAAACAATTCAAATGACAAAGTATTATTCATTGGAAACCACAGAGTTCAGAATGAAG";
59        let lhs_frequencies = Kmer::frequency_vector(k, &lhs_seq);
60
61        let rhs_seq = "CCACGACACATCACAGTCAAATTTCTAAAAATTAAGGACCACACAAACGCCTTAAAAACAGAGAAAGAAAAATGACATATTACCTATACTGCAAAAACAATTCAAATGACAAAGTATTATTCATTGGCAACCACAGAGTTCAGAATGAAG";
62        let rhs_frequencies = Kmer::frequency_vector(k, &rhs_seq);
63
64        let lhs_kmer_frequencies = frequency_vector_iter(lhs_frequencies.iter());
65        let rhs_kmer_frequencies = frequency_vector_iter(rhs_frequencies.iter());
66        let d = dot(lhs_kmer_frequencies, rhs_kmer_frequencies);
67        assert_eq!(d, 188.0);
68    }
69}
kmers/dot.rs

kmers/
dot.rs
