-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheuclideanmeasure.cpp
60 lines (52 loc) · 1.66 KB
/
euclideanmeasure.cpp
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
/*!
* @brief
*
* Copyright (C) 2018 Kenneth Ingham
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "euclideanmeasure.h"
long double
euclideanmeasure::compare(const FastaRecord& a, const FastaRecord& b)
{
long double dist = 0.0;
kmerset* ksa = get_counts(a);
kmerset* ksb = get_counts(b);
// a -> b
for (auto kmi : *ksa) {
int t;
kmer_t km = kmi.first;
if (ksb->count(km) > 0)
t = ksa->at(km) - ksb->at(km);
else
t = ksa->at(km);
dist += t*t;
}
// b -> a, but only for the kmers that are in b and not in a
for (auto kmi : *ksb) {
int t;
kmer_t km = kmi.first;
if (ksa->count(km) == 0) {
t = ksb->at(km);
dist += t*t;
}
}
// mapped into [0,1]
//return dist == 0 ? 0 : 1.0 - 1.0/sqrt(dist);
// Actual Euclidean distance
//return sqrt(dist);
// Save time and do not do the sqrt, since the relative distances
// remain the same, even if the absolute distances are different.
return dist;
};