-
Notifications
You must be signed in to change notification settings - Fork 55
/
Copy pathmh.hpp
164 lines (140 loc) · 4.45 KB
/
mh.hpp
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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
#include<vector>
#include<cstring>
#include<math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include "util.hpp"
void sample_cons_params(struct node nodes[],struct config conf,gsl_rng *rand,int tp);
double multi_param_post(struct node nodes[], struct datum data[], int old,struct config conf);
double param_post(struct node nodes[], struct datum data[], int old,struct config conf, int tp);
void update_params(struct node nodes[], struct config conf);
void get_pi(struct node nodes[], double pi[], struct config conf, int old, int tp);
void load_ssm_data(char fname[], struct datum data[], struct config conf);
void load_cnv_data(char fname[], struct datum data[], struct config conf);
void load_data_states(char fname[], struct datum data[], struct node nodes[], struct config conf);
void load_tree(char fname[], struct node nodes[], struct config conf);
void write_params(char fname[], struct node nodes[], struct config conf);
void mh_loop(struct node nodes[], struct datum data[], char* fname, struct config conf);
struct config{
int MH_ITR;
float MH_STD;
int N_SSM_DATA; // no. of data points
int N_CNV_DATA; // no. of data points
int NNODES; // no. of nodes in the tree
int TREE_HEIGHT;
int NTPS; // no. of samples / time points
};
struct state{
struct node* nd;
int nr,nv;
};
struct node{
int id;
vector<double> param,pi;
vector<double> param1,pi1; // dummy
int ndata;
vector<int> dids;
int nchild;
vector<int> cids; // children ids
int ht;
};
struct datum{
int id;
vector<int> a,d;
double mu_r,mu_v;
vector<double> log_bin_norm_const;//log_bin_coeff(d,a);
struct datum* cnv; // for SSM datum, this is a pointer to its CNV datum
//int cnv;// just an indicator for cnv or ssm datum
// this is used to compute the binomial parameter
vector <struct state> states1, states2, states3, states4; // maternal and paternal state
double log_ll1111(vector<double> phi, int old){
double llh = 0.0;
for(int tp=0; tp<phi.size();tp++)
llh+=log_complete_ll(phi[tp],mu_r,mu_v,old,tp);
return llh;
}
double log_ll(double phi, int old, int tp){
return log_complete_ll(phi,mu_r,mu_v,old,tp);
}
double log_complete_ll(double phi, double mu_r, double mu_v, int old, int tp){
double nr=0;
double nv=0;
double mu = 0;
double llh = 0;
if(cnv==NULL){ // cnv data
mu = (1 - phi) * mu_r + phi * mu_v;
llh = log_binomial_likelihood(a[tp], d[tp], mu) + log_bin_norm_const[tp];
}
else{ // ssm data
double ll[4]; // maternal and paternal
nr=0;
nv=0;
for(int i=0;i<states1.size();i++){
if(old==0){
nr += (states1[i].nd->pi1[tp])*states1[i].nr;
nv += (states1[i].nd->pi1[tp])*states1[i].nv;
}else{
nr += (states1[i].nd->pi[tp])*states1[i].nr;
nv += (states1[i].nd->pi[tp])*states1[i].nv;
}
}
if(nr+nv>0){
mu = (nv*(1-mu_r) + nr*mu_r)/(nr+nv);
ll[0] = log_binomial_likelihood(a[tp], d[tp], mu) + log(0.25) + log_bin_norm_const[tp];
}else{
ll[0]=log(pow(10,-99));
}
// repetitive...
nr=nv=0;
for(int i=0;i<states2.size();i++){
if(old==0){
nr += (states2[i].nd->pi1[tp])*states2[i].nr;
nv += (states2[i].nd->pi1[tp])*states2[i].nv;
}else{
nr += (states2[i].nd->pi[tp])*states2[i].nr;
nv += (states2[i].nd->pi[tp])*states2[i].nv;
}
}
if(nr+nv>0){
mu = (nv*(1-mu_r) + nr*mu_r)/(nr+nv);
ll[1] = log_binomial_likelihood(a[tp], d[tp], mu) + log(0.25) + log_bin_norm_const[tp];
}else{
ll[1]=log(pow(10,-99));
}
nr=nv=0;
for(int i=0;i<states3.size();i++){
if(old==0){
nr += (states3[i].nd->pi1[tp])*states3[i].nr;
nv += (states3[i].nd->pi1[tp])*states3[i].nv;
}else{
nr += (states3[i].nd->pi[tp])*states3[i].nr;
nv += (states3[i].nd->pi[tp])*states3[i].nv;
}
}
if(nr+nv>0){
mu = (nv*(1-mu_r) + nr*mu_r)/(nr+nv);
ll[2] = log_binomial_likelihood(a[tp], d[tp], mu) + log(0.25) + log_bin_norm_const[tp];
}else{
ll[2]=log(pow(10,-99));
}
nr=nv=0;
for(int i=0;i<states4.size();i++){
if(old==0){
nr += (states4[i].nd->pi1[tp])*states4[i].nr;
nv += (states4[i].nd->pi1[tp])*states4[i].nv;
}else{
nr += (states4[i].nd->pi[tp])*states4[i].nr;
nv += (states4[i].nd->pi[tp])*states4[i].nv;
}
}
if(nr+nv>0){
mu = (nv*(1-mu_r) + nr*mu_r)/(nr+nv);
ll[3] = log_binomial_likelihood(a[tp], d[tp], mu) + log(0.25) + log_bin_norm_const[tp];
}else{
ll[3]=log(pow(10,-99));
}
llh = logsumexp(ll,4);
}
return llh;
}
};