1 // -*-mode:c++; c-style:k&r; c-basic-offset:4;-*-
2 //
3 // Copyright 2014, Julian Catchen <jcatchen@uoregon.edu>
4 //
5 // This file is part of Stacks.
6 //
7 // Stacks is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // Stacks is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with Stacks. If not, see <http://www.gnu.org/licenses/>.
19 //
20
21 #ifndef __SMOOTHING_UTILS_H__
22 #define __SMOOTHING_UTILS_H__
23
24 #include <cmath>
25 #include <vector>
26
27 extern double sigma;
28
29 inline
30 double *
calc_weights()31 calc_weights()
32 {
33 //
34 // Calculate weights for window smoothing operations.
35 //
36 // For each genomic region centered on a nucleotide position c, the contribution of the population
37 // genetic statistic at position p to the region average was weighted by the Gaussian function:
38 // exp( (-1 * (p - c)^2) / (2 * sigma^2))
39 //
40 int limit = 3 * sigma;
41 double *weights = new double[limit + 1];
42
43 for (int i = 0; i <= limit; i++)
44 weights[i] = exp((-1 * pow(i, 2)) / (2 * pow(sigma, 2)));
45
46 return weights;
47 }
48
49 template<class StatT>
50 int
determine_window_limits(vector<StatT * > & sites,uint center_bp,uint & pos_l,uint & pos_u)51 determine_window_limits(vector<StatT *> &sites, uint center_bp, uint &pos_l, uint &pos_u)
52 {
53 int limit = 3 * sigma;
54 int limit_l = center_bp - limit > 0 ? center_bp - limit : 0;
55 int limit_u = center_bp + limit;
56
57 while (pos_l < sites.size()) {
58 if (sites[pos_l] == NULL) {
59 pos_l++;
60 } else {
61 if (sites[pos_l]->bp < limit_l)
62 pos_l++;
63 else
64 break;
65 }
66 }
67 while (pos_u < sites.size()) {
68 if (sites[pos_u] == NULL) {
69 pos_u++;
70 } else {
71 if (sites[pos_u]->bp < limit_u)
72 pos_u++;
73 else
74 break;
75 }
76 }
77 return 0;
78 }
79
80 #endif // __SMOOTHING_UTILS_H__
81