-#define PREPARE_BF_F84\
- A = (BF[0]*BF[2])/(BF[0] + BF[2]) + (BF[1]*BF[3])/(BF[1] + BF[3]);\
- B = BF[0]*BF[2] + BF[1]*BF[3];\
- C = (BF[0] + BF[2])*(BF[1] + BF[3]);
-
-#define COMPUTE_DIST_F84\
- P = ((double) Ns/L);\
- Q = ((double) (Nd - Ns)/L);\
- d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\
- if (*variance) {\
- t1 = A*C;\
- t2 = C*P/2;\
- t3 = (A - B)*Q/2;\
- a = t1/(t1 - t2 - t3);\
- b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\
- var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/2;\
- }
-
-void distDNA_F84(unsigned char *x, int *n, int *s, double *d,
- double *BF, int *variance, double *var)
-{
- int i1, i2, Nd, Ns, L, target, s1, s2;
- double P, Q, A, B, C, a, b, t1, t2, t3;
-
- PREPARE_BF_F84
- L = *s;
-
- target = 0;
- for (i1 = 1; i1 < *n; i1++) {
- for (i2 = i1 + 1; i2 <= *n; i2++) {
- Nd = Ns = 0;
- for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
- COUNT_TS_TV
- }
- COMPUTE_DIST_F84
- target++;
- }
- }
-}
-
-void distDNA_F84_pairdel(unsigned char *x, int *n, int *s, double *d,
- double *BF, int *variance, double *var)
-{
- int i1, i2, Nd, Ns, L, target, s1, s2;
- double P, Q, A, B, C, a, b, t1, t2, t3;
-
- PREPARE_BF_F84
-
- target = 0;
- for (i1 = 1; i1 < *n; i1++) {
- for (i2 = i1 + 1; i2 <= *n; i2++) {
- Nd = Ns = L = 0;
- for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
- CHECK_PAIRWISE_DELETION
- COUNT_TS_TV
- }
- COMPUTE_DIST_F84
- target++;
- }
- }
-}