mapgd  0.4
A program for the Maximum-likelihood analysis of population genomic data.
 All Data Structures Functions Variables Friends Groups Pages
newton-method-theta.h
1 /*This code was automatically generated by extras/automatic_src/make_newton_theta.py*/
2 
3 #include "allele.h"
4 #include "quartet.h"
5 #include "typedef.h"
6 
7 inline float_t H0 (const quartet_t &q, const Allele &a) {
8  const float_t M=q.base[a.major];
9  const float_t m=q.base[a.minor];
10  const float_t E=q.base[a.e1]+q.base[a.e2];
11  return -9.375e-8*pow(a.freq, 5) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
12 ;
13 }
14 
15 inline float_t H1 (const quartet_t &q, const Allele &a) {
16  const float_t M=q.base[a.major];
17  const float_t m=q.base[a.minor];
18  const float_t E=q.base[a.e1]+q.base[a.e2];
19  return -9.375e-8*pow(a.error, 5) + (-0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
20 ;
21 }
22 
23 inline float_t H2 (const quartet_t &q, const Allele &a) {
24  const float_t M=q.base[a.major];
25  const float_t m=q.base[a.minor];
26  const float_t E=q.base[a.e1]+q.base[a.e2];
27  return -9.375e-8*pow(a.f, 5) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
28 ;
29 }
30 
31 inline float_t J00 (const quartet_t &q, const Allele &a) {
32  const float_t M=q.base[a.major];
33  const float_t m=q.base[a.minor];
34  const float_t E=q.base[a.e1]+q.base[a.e2];
35  return -4.6875e-7*pow(a.freq, 4) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 2.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(2*a.freq)/pow(exp(a.freq) + 1.0, 3) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 2.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(2*a.freq)/pow(exp(a.freq) + 1.0, 3) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)) - 2.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(2*a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 3)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0)) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))*(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2)
36 ;
37 }
38 
39 inline float_t J01 (const quartet_t &q, const Allele &a) {
40  const float_t M=q.base[a.major];
41  const float_t m=q.base[a.minor];
42  const float_t E=q.base[a.e1]+q.base[a.e2];
43  return (0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.error)*exp(a.freq)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.error)*exp(a.freq)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)*exp(a.freq)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0)) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))*(0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) - 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) + 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2)
44 ;
45 }
46 
47 inline float_t J02 (const quartet_t &q, const Allele &a) {
48  const float_t M=q.base[a.major];
49  const float_t m=q.base[a.minor];
50  const float_t E=q.base[a.e1]+q.base[a.e2];
51  return (-1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
52 ;
53 }
54 
55 inline float_t J10 (const quartet_t &q, const Allele &a) {
56  const float_t M=q.base[a.major];
57  const float_t m=q.base[a.minor];
58  const float_t E=q.base[a.e1]+q.base[a.e2];
59  return (0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.error)*exp(a.freq)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.error)*exp(a.freq)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)*exp(a.freq)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*pow(exp(a.freq) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.freq)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0)) + (1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))*(-0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2)
60 ;
61 }
62 
63 inline float_t J11 (const quartet_t &q, const Allele &a) {
64  const float_t M=q.base[a.major];
65  const float_t m=q.base[a.minor];
66  const float_t E=q.base[a.e1]+q.base[a.e2];
67  return -4.6875e-7*pow(a.error, 4) + (0.0625*pow(E, 2)*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*pow(4.0*exp(a.error) + 4.0, 2)*exp(2*a.error)/(pow(exp(a.error) + 1, 4)*(exp(a.freq) + 1.0)) - 0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 1.0*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(2*a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.5*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/(pow(exp(a.error) + 1, 3)*(exp(a.freq) + 1.0)) - 0.125*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 4)*(exp(a.freq) + 1.0)) + 0.5625*pow(M, 2)*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/(pow(1.0 - 0.75/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)*(exp(a.f) + 1.0)) - 0.375*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 4)*(exp(a.f) + 1.0)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) - 1.5*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 3)*(exp(a.f) + 1.0)) - 0.5625*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/(pow(1.0 - 0.75/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)*(exp(a.f) + 1.0)) + 0.5625*pow(m, 2)*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/(pow(1.0 - 0.75/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)) - 0.375*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 4)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) - 1.5*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 3)) - 0.5625*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.error)/(pow(1.0 - 0.75/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)) + 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.5*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(2*a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 3)*(exp(a.freq) + 1.0)) + 0.0625*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*pow(M + m, 2)*exp(2*a.error)/(pow(0.5 - 0.25/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)*(exp(a.freq) + 1.0)) - 0.0625*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(2*a.error)/(pow(0.5 - 0.25/(exp(a.error) + 1), 2)*pow(exp(a.error) + 1, 4)*(exp(a.freq) + 1.0)) + 0.0625*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*pow(E + M, 2)*pow(4.0*exp(a.error) + 4.0, 2)*exp(2*a.error)/pow(exp(a.error) + 1, 4) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*exp(2*a.error)/pow(exp(a.error) + 1, 2) + 0.5*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/pow(exp(a.error) + 1, 3) + 0.0625*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*pow(E + m, 2)*pow(4.0*exp(a.error) + 4.0, 2)*exp(2*a.error)/(pow(exp(a.error) + 1, 4)*(exp(a.f) + 1.0)) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*exp(2*a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) + 0.5*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(2*a.error)/(pow(exp(a.error) + 1, 3)*(exp(a.f) + 1.0)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0)) + (-0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))*(0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) - 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) + 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2)
68 ;
69 }
70 
71 inline float_t J12 (const quartet_t &q, const Allele &a) {
72  const float_t M=q.base[a.major];
73  const float_t m=q.base[a.minor];
74  const float_t E=q.base[a.e1]+q.base[a.e2];
75  return (-1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(-0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) - 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (-0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)*exp(a.f)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)*exp(a.f)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.f)/(pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.f)/(pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
76 ;
77 }
78 
79 inline float_t J20 (const quartet_t &q, const Allele &a) {
80  const float_t M=q.base[a.major];
81  const float_t m=q.base[a.minor];
82  const float_t E=q.base[a.e1]+q.base[a.e2];
83  return (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
84 ;
85 }
86 
87 inline float_t J21 (const quartet_t &q, const Allele &a) {
88  const float_t M=q.base[a.major];
89  const float_t m=q.base[a.minor];
90  const float_t E=q.base[a.e1]+q.base[a.e2];
91  return (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(0.25*E*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) - 0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)) - 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*(M + m)*exp(a.error)/((0.5 - 0.25/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*(exp(a.freq) + 1.0)) + 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)/pow(exp(a.error) + 1, 2) + 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)/(pow(exp(a.error) + 1, 2)*(exp(a.f) + 1.0)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (-0.75*M*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)*exp(a.f)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) + 0.75*m*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.error)*exp(a.f)/((1.0 - 0.75/(exp(a.error) + 1))*pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) - 0.25*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + M)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.f)/(pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)) + 0.25*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*(E + m)*(4.0*exp(a.error) + 4.0)*exp(a.error)*exp(a.f)/(pow(exp(a.error) + 1, 2)*pow(exp(a.f) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
92 ;
93 }
94 
95 inline float_t J22 (const quartet_t &q, const Allele &a) {
96  const float_t M=q.base[a.major];
97  const float_t m=q.base[a.minor];
98  const float_t E=q.base[a.e1]+q.base[a.e2];
99  return -4.6875e-7*pow(a.f, 4) + (-1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 2.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.f)/pow(exp(a.f) + 1.0, 3) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 2.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.f)/pow(exp(a.f) + 1.0, 3))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
100 ;
101 }
102 
103 inline float_t lnL_NR (const quartet_t &q, const Allele &a) {
104  const float_t M=q.base[a.major];
105  const float_t m=q.base[a.minor];
106  const float_t E=q.base[a.e1]+q.base[a.e2];
107  return -1.5625e-8*pow(a.error, 6) - 1.5625e-8*pow(a.f, 6) - 1.5625e-8*pow(a.freq, 6) + log((1.0*pow(0.25/(exp(a.error) + 1), E)*pow((0.5*exp(a.error) + 0.25)/(exp(a.error) + 1), M + m)*(exp(a.f) + 1.0) + (pow(0.25/(exp(a.error) + 1), E + M)*pow((1.0*exp(a.error) + 0.25)/(exp(a.error) + 1), m)*exp(a.f) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow((1.0*exp(a.error) + 0.25)/(exp(a.error) + 1), M))*exp(a.freq))/((exp(a.f) + 1.0)*(exp(a.freq) + 1.0)))
108 ;
109 }
110 
111 inline float_t H0R (const quartet_t &q, const Allele &a) {
112  const float_t M=q.base[a.major];
113  const float_t m=q.base[a.minor];
114  const float_t E=q.base[a.e1]+q.base[a.e2];
115  return -9.375e-8*pow(a.freq, 5) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
116 ;
117 }
118 
119 inline float_t H1R (const quartet_t &q, const Allele &a) {
120  const float_t M=q.base[a.major];
121  const float_t m=q.base[a.minor];
122  const float_t E=q.base[a.e1]+q.base[a.e2];
123  return -9.375e-8*pow(a.f, 5) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
124 ;
125 }
126 
127 inline float_t J00R (const quartet_t &q, const Allele &a) {
128  const float_t M=q.base[a.major];
129  const float_t m=q.base[a.minor];
130  const float_t E=q.base[a.e1]+q.base[a.e2];
131  return -4.6875e-7*pow(a.freq, 4) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 2.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(2*a.freq)/pow(exp(a.freq) + 1.0, 3) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 2.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(2*a.freq)/pow(exp(a.freq) + 1.0, 3) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)) - 2.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(2*a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 3)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0)) + (-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))*(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2)
132 ;
133 }
134 
135 inline float_t J01R (const quartet_t &q, const Allele &a) {
136  const float_t M=q.base[a.major];
137  const float_t m=q.base[a.minor];
138  const float_t E=q.base[a.e1]+q.base[a.e2];
139  return (-1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(-1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
140 ;
141 }
142 
143 inline float_t J10R (const quartet_t &q, const Allele &a) {
144  const float_t M=q.base[a.major];
145  const float_t m=q.base[a.minor];
146  const float_t E=q.base[a.e1]+q.base[a.e2];
147  return (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*exp(a.freq)/pow(exp(a.freq) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.freq)/((exp(a.f) + 1.0)*pow(exp(a.freq) + 1.0, 2)))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*exp(a.f)*exp(a.freq)/(pow(exp(a.f) + 1.0, 2)*pow(exp(a.freq) + 1.0, 2)))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
148 ;
149 }
150 
151 inline float_t J11R (const quartet_t &q, const Allele &a) {
152  const float_t M=q.base[a.major];
153  const float_t m=q.base[a.minor];
154  const float_t E=q.base[a.e1]+q.base[a.e2];
155  return -4.6875e-7*pow(a.f, 4) + (-1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))*(1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2))/pow(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0), 2) + (1.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) - 2.0*pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.f)/pow(exp(a.f) + 1.0, 3) - 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(a.f)/pow(exp(a.f) + 1.0, 2) + 2.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))*exp(2*a.f)/pow(exp(a.f) + 1.0, 3))/(1.0*pow(0.25/(exp(a.error) + 1), E)*pow(0.5 - 0.25/(exp(a.error) + 1), M + m)/(exp(a.freq) + 1.0) + pow(0.25/(exp(a.error) + 1), E + M)*pow(1.0 - 0.75/(exp(a.error) + 1), m)*(1 - 1.0/(exp(a.f) + 1.0))*(1 - 1.0/(exp(a.freq) + 1.0)) + 1.0*pow(0.25/(exp(a.error) + 1), E + m)*pow(1.0 - 0.75/(exp(a.error) + 1), M)*(1 - 1.0/(exp(a.freq) + 1.0))/(exp(a.f) + 1.0))
156 ;
157 }
158 
Allele::e2
gt_t e2
identity of error2.
Definition: allele.h:45
quartet
A class that stores quartet information.
Definition: quartet.h:16
quartet::base
count_t base[5]
The count of the number of occurnaces of bases. What nucleotieds this cout represents is stored at th...
Definition: quartet.h:17
Allele::f
float_t f
HW statistic.
Definition: allele.h:61
Allele::minor
gt_t minor
identity of minor allele.
Definition: allele.h:42
Allele::error
float_t error
ml error rate.
Definition: allele.h:47
Allele::e1
gt_t e1
identity of error1
Definition: allele.h:44
Allele::freq
float_t freq
frequency of major allele.
Definition: allele.h:39
Allele::major
gt_t major
identity major allele.
Definition: allele.h:43
Allele
Summary statistics from the allele command.
Definition: allele.h:19
Generated on Mon Oct 1 2018 13:55:20 for mapgd by   doxygen 1.8.6