mapgd  0.4
A program for the Maximum-likelihood analysis of population genomic data.
 All Data Structures Functions Variables Friends Groups Pages
newton-method-ld.h
1 /*This code was automatically generated by extras/automatic_src/make_tangent_ld.py*/
2 
3 #ifndef NEWTON_LD
4 #define NEWTON_LD
5 
6 #include "population.h"
7 #include "typedef.h"
8 #include "constants.h"
9 
10 
11 float_t H0 (const Population &P1, const Population &P2, const float_t &D);
12 
13 float_t J00 (const Population &P1, const Population &P2, const float_t &D);
14 
15 inline float_t lnL_NR (const Population &P1, const Population &P2, const float_t &D) {
16  float_t sum=0;
17  float_t Q1=P1.m;
18  float_t Q2=P2.m;
19  float_t f=(P1.f+P2.f)/2.;
20  int p_size=P1.likelihoods.size();
21  for (int x=0; x<p_size; ++x) {
22  float_t MM1=P1.likelihoods[x].MM;
23  float_t Mm1=P1.likelihoods[x].Mm;
24  float_t mm1=P1.likelihoods[x].mm;
25  float_t MM2=P2.likelihoods[x].MM;
26  float_t Mm2=P2.likelihoods[x].Mm;
27  float_t mm2=P2.likelihoods[x].mm;
28  if (P1.likelihoods[x].N > 0 and P2.likelihoods[x].N > 0){
29  sum+=log(MM1*MM2*(D + (Q1 - 1.0)*(Q2 - 1.0))*(D - f*(D + (Q1 - 1.0)*(Q2 - 1.0) - 1) + (Q1 - 1.0)*(Q2 - 1.0)) + 2*MM1*Mm2*(D + Q2*(Q1 - 1.0))*(D + (Q1 - 1.0)*(Q2 - 1.0))*(f - 1) + MM1*mm2*(D + Q2*(Q1 - 1.0))*(D + Q2*(Q1 - 1.0) - f*(D + Q2*(Q1 - 1.0) + 1)) + 2*MM2*Mm1*(D + Q1*(Q2 - 1.0))*(D + (Q1 - 1.0)*(Q2 - 1.0))*(f - 1) + MM2*mm1*(D + Q1*(Q2 - 1.0))*(D + Q1*(Q2 - 1.0) - f*(D + Q1*(Q2 - 1.0) + 1)) - 2*Mm1*Mm2*(f - 1)*((D + Q1*Q2)*(D + (Q1 - 1.0)*(Q2 - 1.0)) + (D + Q1*(Q2 - 1.0))*(D + Q2*(Q1 - 1.0))) + 2*Mm1*mm2*(D + Q1*Q2)*(D + Q2*(Q1 - 1.0))*(f - 1) + 2*Mm2*mm1*(D + Q1*Q2)*(D + Q1*(Q2 - 1.0))*(f - 1) + mm1*mm2*(D + Q1*Q2)*(D + Q1*Q2 - f*(D + Q1*Q2 - 1)))
30  ;}
31  }
32  return sum;
33 }
34 
35 inline void set_c27 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
36  const float_t Q1=std::get<0>(d).m;
37  const float_t Q2=std::get<1>(d).m;
38  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
39  const float_t D=std::get<2>(d);
40  con.c[27]= -2*D + 2*Q1*(-Q2 + 1.0) ;
41 }
42 
43 inline void set_c26 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
44  const float_t Q1=std::get<0>(d).m;
45  const float_t Q2=std::get<1>(d).m;
46  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
47  const float_t D=std::get<2>(d);
48  con.c[26]= 2*D + 2*(-Q1 + 1.0)*(-Q2 + 1.0) ;
49 }
50 
51 inline void set_c25 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
52  const float_t Q1=std::get<0>(d).m;
53  const float_t Q2=std::get<1>(d).m;
54  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
55  const float_t D=std::get<2>(d);
56  con.c[25]= -D + Q1*(-Q2 + 1.0) ;
57 }
58 
59 inline void set_c23 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
60  const float_t Q1=std::get<0>(d).m;
61  const float_t Q2=std::get<1>(d).m;
62  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
63  const float_t D=std::get<2>(d);
64  con.c[23]= D + (-Q1 + 1.0)*(-Q2 + 1.0) ;
65 }
66 
67 inline void set_c11 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
68  const float_t Q1=std::get<0>(d).m;
69  const float_t Q2=std::get<1>(d).m;
70  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
71  const float_t D=std::get<2>(d);
72  con.c[11]= D + Q1*(Q2 - 1.0) ;
73 }
74 
75 inline void set_c10 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
76  const float_t Q1=std::get<0>(d).m;
77  const float_t Q2=std::get<1>(d).m;
78  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
79  const float_t D=std::get<2>(d);
80  con.c[10]= D + Q2*(Q1 - 1.0) ;
81 }
82 
83 inline void set_c6 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
84  const float_t Q1=std::get<0>(d).m;
85  const float_t Q2=std::get<1>(d).m;
86  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
87  const float_t D=std::get<2>(d);
88  con.c[6]= D + Q1*Q2 ;
89 }
90 
91 inline void set_c0 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
92  const float_t Q1=std::get<0>(d).m;
93  const float_t Q2=std::get<1>(d).m;
94  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
95  const float_t D=std::get<2>(d);
96  con.c[0]= D + (Q1 - 1.0)*(Q2 - 1.0) ;
97 }
98 
99 inline void set_c33 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
100  const float_t Q1=std::get<0>(d).m;
101  const float_t Q2=std::get<1>(d).m;
102  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
103  const float_t D=std::get<2>(d);
104  con.c[33]= 8*D + 2*Q1*Q2 - 2*Q1*(-Q2 + 1.0) - 2*Q2*(-Q1 + 1.0) + 2*(-Q1 + 1.0)*(-Q2 + 1.0) ;
105 }
106 
107 inline void set_c32 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
108  const float_t Q1=std::get<0>(d).m;
109  const float_t Q2=std::get<1>(d).m;
110  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
111  const float_t D=std::get<2>(d);
112  con.c[32]= 2*D + 2*Q1*Q2 - con.c[6]*f + f*(-D - Q1*Q2 + 1) ;
113 }
114 
115 inline void set_c31 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
116  const float_t Q1=std::get<0>(d).m;
117  const float_t Q2=std::get<1>(d).m;
118  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
119  const float_t D=std::get<2>(d);
120  con.c[31]= 2*D - 2*Q1*(-Q2 + 1.0) + con.c[25]*f - f*(D - Q1*(-Q2 + 1.0) + 1) ;
121 }
122 
123 inline void set_c30 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
124  const float_t Q1=std::get<0>(d).m;
125  const float_t Q2=std::get<1>(d).m;
126  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
127  const float_t D=std::get<2>(d);
128  con.c[30]= 2*D - 2*Q2*(-Q1 + 1.0) + f*(-D + Q2*(-Q1 + 1.0)) - f*(D - Q2*(-Q1 + 1.0) + 1) ;
129 }
130 
131 inline void set_c29 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
132  const float_t Q1=std::get<0>(d).m;
133  const float_t Q2=std::get<1>(d).m;
134  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
135  const float_t D=std::get<2>(d);
136  con.c[29]= 2*D - con.c[23]*f + f*(-D - (-Q1 + 1.0)*(-Q2 + 1.0) + 1) + 2*(-Q1 + 1.0)*(-Q2 + 1.0) ;
137 }
138 
139 inline void set_c28 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
140  const float_t Q1=std::get<0>(d).m;
141  const float_t Q2=std::get<1>(d).m;
142  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
143  const float_t D=std::get<2>(d);
144  con.c[28]= -2*f + 2 ;
145 }
146 
147 inline void set_c24 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
148  const float_t Q1=std::get<0>(d).m;
149  const float_t Q2=std::get<1>(d).m;
150  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
151  const float_t D=std::get<2>(d);
152  con.c[24]= -2*D + 2*Q2*(-Q1 + 1.0) ;
153 }
154 
155 inline void set_c22 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
156  const float_t Q1=std::get<0>(d).m;
157  const float_t Q2=std::get<1>(d).m;
158  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
159  const float_t D=std::get<2>(d);
160  con.c[22]= con.c[26]*con.c[6] + con.c[27]*(-D + Q2*(-Q1 + 1.0)) ;
161 }
162 
163 inline void set_c21 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
164  const float_t Q1=std::get<0>(d).m;
165  const float_t Q2=std::get<1>(d).m;
166  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
167  const float_t D=std::get<2>(d);
168  con.c[21]= -f + 1 ;
169 }
170 
171 inline void set_c20 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
172  const float_t Q1=std::get<0>(d).m;
173  const float_t Q2=std::get<1>(d).m;
174  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
175  const float_t D=std::get<2>(d);
176  con.c[20]= pow(con.c[6], 2) + con.c[6]*f*(-D - Q1*Q2 + 1) ;
177 }
178 
179 inline void set_c19 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
180  const float_t Q1=std::get<0>(d).m;
181  const float_t Q2=std::get<1>(d).m;
182  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
183  const float_t D=std::get<2>(d);
184  con.c[19]= pow(con.c[25], 2) + con.c[25]*f*(D - Q1*(-Q2 + 1.0) + 1) ;
185 }
186 
187 inline void set_c18 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
188  const float_t Q1=std::get<0>(d).m;
189  const float_t Q2=std::get<1>(d).m;
190  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
191  const float_t D=std::get<2>(d);
192  con.c[18]= f*(-D + Q2*(-Q1 + 1.0))*(D - Q2*(-Q1 + 1.0) + 1) + pow(-D + Q2*(-Q1 + 1.0), 2) ;
193 }
194 
195 inline void set_c17 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
196  const float_t Q1=std::get<0>(d).m;
197  const float_t Q2=std::get<1>(d).m;
198  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
199  const float_t D=std::get<2>(d);
200  con.c[17]= pow(con.c[23], 2) + con.c[23]*f*(-D - (-Q1 + 1.0)*(-Q2 + 1.0) + 1) ;
201 }
202 
203 inline void set_c16 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
204  const float_t Q1=std::get<0>(d).m;
205  const float_t Q2=std::get<1>(d).m;
206  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
207  const float_t D=std::get<2>(d);
208  con.c[16]= 4*D + Q1*Q2 + Q1*(Q2 - 1.0) + Q2*(Q1 - 1.0) + (Q1 - 1.0)*(Q2 - 1.0) ;
209 }
210 
211 inline void set_c15 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
212  const float_t Q1=std::get<0>(d).m;
213  const float_t Q2=std::get<1>(d).m;
214  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
215  const float_t D=std::get<2>(d);
216  con.c[15]= 2*D + 2*Q1*Q2 - con.c[6]*f - f*(D + Q1*Q2 - 1) ;
217 }
218 
219 inline void set_c14 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
220  const float_t Q1=std::get<0>(d).m;
221  const float_t Q2=std::get<1>(d).m;
222  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
223  const float_t D=std::get<2>(d);
224  con.c[14]= 2*D + 2*Q1*(Q2 - 1.0) - con.c[11]*f - f*(D + Q1*(Q2 - 1.0) + 1) ;
225 }
226 
227 inline void set_c13 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
228  const float_t Q1=std::get<0>(d).m;
229  const float_t Q2=std::get<1>(d).m;
230  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
231  const float_t D=std::get<2>(d);
232  con.c[13]= 2*D + 2*Q2*(Q1 - 1.0) - con.c[10]*f - f*(D + Q2*(Q1 - 1.0) + 1) ;
233 }
234 
235 inline void set_c12 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
236  const float_t Q1=std::get<0>(d).m;
237  const float_t Q2=std::get<1>(d).m;
238  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
239  const float_t D=std::get<2>(d);
240  con.c[12]= 2*D - con.c[0]*f - f*(D + (Q1 - 1.0)*(Q2 - 1.0) - 1) + 2*(Q1 - 1.0)*(Q2 - 1.0) ;
241 }
242 
243 inline void set_c9 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
244  const float_t Q1=std::get<0>(d).m;
245  const float_t Q2=std::get<1>(d).m;
246  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
247  const float_t D=std::get<2>(d);
248  con.c[9]= con.c[0]*con.c[6] + con.c[10]*con.c[11] ;
249 }
250 
251 inline void set_c8 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
252  const float_t Q1=std::get<0>(d).m;
253  const float_t Q2=std::get<1>(d).m;
254  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
255  const float_t D=std::get<2>(d);
256  con.c[8]= f - 1 ;
257 }
258 
259 inline void set_c7 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
260  const float_t Q1=std::get<0>(d).m;
261  const float_t Q2=std::get<1>(d).m;
262  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
263  const float_t D=std::get<2>(d);
264  con.c[7]= -D - Q1*Q2 + f*(D + Q1*Q2 - 1) ;
265 }
266 
267 inline void set_c5 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
268  const float_t Q1=std::get<0>(d).m;
269  const float_t Q2=std::get<1>(d).m;
270  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
271  const float_t D=std::get<2>(d);
272  con.c[5]= -D + Q1*(-Q2 + 1.0) + f*(D + Q1*(Q2 - 1.0) + 1) ;
273 }
274 
275 inline void set_c4 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
276  const float_t Q1=std::get<0>(d).m;
277  const float_t Q2=std::get<1>(d).m;
278  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
279  const float_t D=std::get<2>(d);
280  con.c[4]= D - Q1*(-Q2 + 1.0) ;
281 }
282 
283 inline void set_c3 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
284  const float_t Q1=std::get<0>(d).m;
285  const float_t Q2=std::get<1>(d).m;
286  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
287  const float_t D=std::get<2>(d);
288  con.c[3]= -D + Q2*(-Q1 + 1.0) + f*(D + Q2*(Q1 - 1.0) + 1) ;
289 }
290 
291 inline void set_c2 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
292  const float_t Q1=std::get<0>(d).m;
293  const float_t Q2=std::get<1>(d).m;
294  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
295  const float_t D=std::get<2>(d);
296  con.c[2]= D - Q2*(-Q1 + 1.0) ;
297 }
298 
299 inline void set_c1 (Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &con, const std::tuple<const Population &, const Population &, const float_t &> &d) {
300  const float_t Q1=std::get<0>(d).m;
301  const float_t Q2=std::get<1>(d).m;
302  const float_t f=(std::get<0>(d).f+std::get<1>(d).f)/2.;
303  const float_t D=std::get<2>(d);
304  con.c[1]= -D + f*(D + (-Q1 + 1.0)*(-Q2 + 1.0) - 1) - (Q1 - 1.0)*(Q2 - 1.0) ;
305 }
306 
307 static void (*cfn_ld[34])(Constants<float, std::tuple<const Population &, const Population &, const float_t &> > &, const std::tuple<const Population &, const Population &, const float_t &> &)={&set_c27 , &set_c26
308 , &set_c25
309 , &set_c23
310 , &set_c11
311 , &set_c10
312 , &set_c6
313 , &set_c0
314 , &set_c33
315 , &set_c32
316 , &set_c31
317 , &set_c30
318 , &set_c29
319 , &set_c28
320 , &set_c24
321 , &set_c22
322 , &set_c21
323 , &set_c20
324 , &set_c19
325 , &set_c18
326 , &set_c17
327 , &set_c16
328 , &set_c15
329 , &set_c14
330 , &set_c13
331 , &set_c12
332 , &set_c9
333 , &set_c8
334 , &set_c7
335 , &set_c5
336 , &set_c4
337 , &set_c3
338 , &set_c2
339 , &set_c1
340 };
341 
342 #define LD_CNTS 34
343 #define LD_ARRAY cfn_ld
344 
345 #endif
Population genotypes.
Definition: population.h:13
Definition: constants.h:7
std::vector< Genotype > likelihoods
Genotypic likelihood.
Definition: population.h:23
float_t m
minor allele frequency
Definition: population.h:27
float_t f
departure from HWE
Definition: population.h:28