-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathExperimentDGE.hpp
More file actions
134 lines (107 loc) · 3.76 KB
/
Copy pathExperimentDGE.hpp
File metadata and controls
134 lines (107 loc) · 3.76 KB
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
#ifndef EXPERIMENT_DGE_HPP
#define EXPERIMENT_DGE_HPP
#include "Experiment.hpp"
template <typename Type>
template <typename Fnc>
Matrix Experiment<Type>::local_dge(Fnc fnc,
const features_t &global_features) const {
auto spot_weights = marginalize_spots();
Matrix m(G, T);
#pragma omp parallel for if (DO_PARALLEL)
for (size_t g = 0; g < G; ++g)
for (size_t t = 0; t < T; ++t)
m(g, t) = local_dge_sub(fnc, global_features, g, t, spot_weights[t]);
return m;
}
template <typename Type>
template <typename Fnc>
Float Experiment<Type>::local_dge_sub(Fnc fnc,
const features_t &global_features,
size_t g, size_t t, Float theta_,
Float p) const {
const Float eps = 1e-6;
Float q = 1 - p;
Float lp = log(p);
Float lq = log(q);
double mi = 0;
size_t x = 0;
double cumsum = 0;
while (cumsum < 2 - 2 * eps) {
Float l1 = log_negative_binomial(x, global_features.prior.r(g, t),
baseline_phi(g) * phi(g, t) * theta_,
global_features.prior.p(g, t));
Float l2 = log_negative_binomial(x, global_features.prior.r(g, t),
fnc(baseline_phi(g), phi(g, t)) * theta_,
global_features.prior.p(g, t));
cumsum += exp(l1);
cumsum += exp(l2);
l1 += lp;
l2 += lq;
Float lz = logSumExp(l1, l2);
Float p1 = exp(l1);
Float p2 = exp(l2);
mi += p1 * (l1 - lz - lp) + p2 * (l2 - lz - lq);
x++;
}
// LOG(debug) << "t = " << t << " g = " << g << " x = " << x;
return mi / log(2.0);
}
template <typename Type>
Matrix Experiment<Type>::pairwise_dge(const features_t &global_features) const {
size_t T_ = 0;
for (size_t t1 = 0; t1 < T; ++t1)
for (size_t t2 = t1 + 1; t2 < T; ++t2)
T_++;
Matrix m(G, T_);
size_t t_ = 0;
for (size_t t1 = 0; t1 < T; ++t1)
for (size_t t2 = t1 + 1; t2 < T; ++t2) {
auto v = pairwise_dge_sub(global_features, t1, t2);
for (size_t g = 0; g < G; ++g)
m(g, t_) = v(g);
t_++;
}
return m;
}
template <typename Type>
Vector Experiment<Type>::pairwise_dge_sub(const features_t &global_features,
size_t t1, size_t t2) const {
LOG(verbose) << "Performing DGE for factor " << t1 << " and factor " << t2;
Vector v(G);
#pragma omp parallel for if (DO_PARALLEL)
for (size_t g = 0; g < G; ++g)
v(g) = pairwise_dge_sub(global_features, t1, t2, g);
return v;
}
template <typename Type>
Float Experiment<Type>::pairwise_dge_sub(const features_t &global_features,
size_t t1, size_t t2, size_t g,
Float theta_, Float p) const {
const Float eps = 1e-6;
Float q = 1 - p;
Float lp = log(p);
Float lq = log(q);
double mi = 0;
size_t x = 0;
double cumsum = 0;
while (cumsum < 2 - 2 * eps) {
Float l1 = log_negative_binomial(x, global_features.prior.r(g, t1),
baseline_phi(g) * phi(g, t1) * theta_,
global_features.prior.p(g, t1));
Float l2 = log_negative_binomial(x, global_features.prior.r(g, t2),
baseline_phi(g) * phi(g, t2) * theta_,
global_features.prior.p(g, t2));
cumsum += exp(l1);
cumsum += exp(l2);
l1 += lp;
l2 += lq;
Float lz = logSumExp(l1, l2);
Float p1 = exp(l1);
Float p2 = exp(l2);
mi += p1 * (l1 - lz - lp) + p2 * (l2 - lz - lq);
x++;
}
// LOG(debug) << "t1 = " << t1 << " t2 = " << t2 << " g = " << g << " x = " << x;
return mi / log(2.0);
}
#endif