1 // This file is part of PLINK 1.90, copyright (C) 2005-2020 Shaun Purcell,
2 // Christopher Chang.
3 //
4 // This program is free software: you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation, either version 3 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 // GNU General Public License for more details.
13 //
14 // You should have received a copy of the GNU General Public License
15 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
16 
17 
18 #include "plink_common.h"
19 
20 #include "plink_stats.h"
21 #include "ipmpar.h"
22 #include "dcdflib.h"
23 
24 // 2^{-40} for now, since 2^{-44} was too small on real data
25 #define FISHER_EPSILON 0.0000000000009094947017729282379150390625
26 
chiprob_p(double xx,double df)27 double chiprob_p(double xx, double df) {
28   int st = 0;
29   int ww = 1;
30   double bnd = 1;
31   double pp;
32   double qq;
33   cdfchi(&ww, &pp, &qq, &xx, &df, &st, &bnd);
34   if (st) {
35     return -9;
36   }
37   return qq;
38 }
39 
inverse_chiprob(double qq,double df)40 double inverse_chiprob(double qq, double df) {
41   double pp = 1 - qq;
42   int32_t st = 0;
43   int32_t ww = 2;
44   double bnd = 1;
45   double xx;
46 
47   if (qq >= 1.0) {
48     return 0;
49   }
50   cdfchi(&ww, &pp, &qq, &xx, &df, &st, &bnd);
51   if (st != 0) {
52     return -9;
53   }
54   return xx;
55 }
56 
calc_tprob(double tt,double df)57 double calc_tprob(double tt, double df) {
58   int32_t st = 0;
59   int32_t ww = 1;
60   double bnd = 1;
61   double pp;
62   double qq;
63   if (!realnum(tt)) {
64     return -9;
65   }
66   tt = fabs(tt);
67   cdft(&ww, &pp, &qq, &tt, &df, &st, &bnd);
68   if (st != 0) {
69     return -9;
70   }
71   return 2 * qq;
72 }
73 
inverse_tprob(double dbl_qq,double df)74 double inverse_tprob(double dbl_qq, double df) {
75   double qq = dbl_qq * 0.5;
76   double pp = 1 - qq;
77   int32_t st = 0;
78   int32_t ww = 2;
79   double bnd = 1;
80   double tt;
81   cdft(&ww, &pp, &qq, &tt, &df, &st, &bnd);
82   if (st != 0) {
83     return -9;
84   }
85   return tt;
86 }
87 
88 // Inverse normal distribution
89 
90 //
91 // Lower tail quantile for standard normal distribution function.
92 //
93 // This function returns an approximation of the inverse cumulative
94 // standard normal distribution function.  I.e., given P, it returns
95 // an approximation to the X satisfying P = Pr{Z <= X} where Z is a
96 // random variable from the standard normal distribution.
97 //
98 // The algorithm uses a minimax approximation by rational functions
99 // and the result has a relative error whose absolute value is less
100 // than 1.15e-9.
101 //
102 // Author:      Peter J. Acklam
103 // Time-stamp:  2002-06-09 18:45:44 +0200
104 // E-mail:      jacklam@math.uio.no
105 // WWW URL:     http://www.math.uio.no/~jacklam
106 //
107 // C implementation adapted from Peter's Perl version
108 
109 // Coefficients in rational approximations.
110 
111 static const double ivn_a[] =
112   {
113     -3.969683028665376e+01,
114     2.209460984245205e+02,
115     -2.759285104469687e+02,
116     1.383577518672690e+02,
117     -3.066479806614716e+01,
118      2.506628277459239e+00
119   };
120 
121 static const double ivn_b[] =
122   {
123     -5.447609879822406e+01,
124     1.615858368580409e+02,
125     -1.556989798598866e+02,
126     6.680131188771972e+01,
127     -1.328068155288572e+01
128   };
129 
130 static const double ivn_c[] =
131   {
132     -7.784894002430293e-03,
133     -3.223964580411365e-01,
134     -2.400758277161838e+00,
135     -2.549732539343734e+00,
136     4.374664141464968e+00,
137      2.938163982698783e+00
138   };
139 
140 static const double ivn_d[] =
141   {
142     7.784695709041462e-03,
143     3.224671290700398e-01,
144     2.445134137142996e+00,
145     3.754408661907416e+00
146   };
147 
148 #define IVN_LOW 0.02425
149 #define IVN_HIGH 0.97575
150 
ltqnorm(double p)151 double ltqnorm(double p) {
152   // assumes 0 < p < 1
153   double q, r;
154 
155   if (p < IVN_LOW) {
156     // Rational approximation for lower region
157     q = sqrt(-2*log(p));
158     return (((((ivn_c[0]*q+ivn_c[1])*q+ivn_c[2])*q+ivn_c[3])*q+ivn_c[4])*q+ivn_c[5]) /
159       ((((ivn_d[0]*q+ivn_d[1])*q+ivn_d[2])*q+ivn_d[3])*q+1);
160   } else if (p > IVN_HIGH) {
161     // Rational approximation for upper region
162     q  = sqrt(-2*log(1-p));
163     return -(((((ivn_c[0]*q+ivn_c[1])*q+ivn_c[2])*q+ivn_c[3])*q+ivn_c[4])*q+ivn_c[5]) /
164       ((((ivn_d[0]*q+ivn_d[1])*q+ivn_d[2])*q+ivn_d[3])*q+1);
165   } else {
166     // Rational approximation for central region
167     q = p - 0.5;
168     r = q*q;
169     return (((((ivn_a[0]*r+ivn_a[1])*r+ivn_a[2])*r+ivn_a[3])*r+ivn_a[4])*r+ivn_a[5])*q /
170       (((((ivn_b[0]*r+ivn_b[1])*r+ivn_b[2])*r+ivn_b[3])*r+ivn_b[4])*r+1);
171   }
172 }
173 
SNPHWE2(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,uint32_t midp)174 double SNPHWE2(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, uint32_t midp) {
175   // This function implements an exact SNP test of Hardy-Weinberg
176   // Equilibrium as described in Wigginton, JE, Cutler, DJ, and
177   // Abecasis, GR (2005) A Note on Exact Tests of Hardy-Weinberg
178   // Equilibrium. American Journal of Human Genetics. 76: 000 - 000.
179   //
180   // The original version was written by Jan Wigginton.
181   //
182   // This version was written by Christopher Chang.  It contains the following
183   // improvements over the original SNPHWE():
184   // - Proper handling of >64k genotypes.  Previously, there was a potential
185   //   integer overflow.
186   // - Detection and efficient handling of floating point overflow and
187   //   underflow.  E.g. instead of summing a tail all the way down, the loop
188   //   stops once the latest increment underflows the partial sum's 53-bit
189   //   precision; this results in a large speedup when max heterozygote count
190   //   >1k.
191   // - No malloc() call.  It's only necessary to keep track of a few partial
192   //   sums.
193   // - Support for the mid-p variant of this test.  See Graffelman J, Moreno V
194   //   (2013) The mid p-value in exact tests for Hardy-Weinberg equilibrium.
195   //
196   // Note that the SNPHWE_t() function below is a lot more efficient for
197   // testing against a p-value inclusion threshold.  SNPHWE2() should only be
198   // used if you need the actual p-value.
199   intptr_t obs_homc;
200   intptr_t obs_homr;
201   if (obs_hom1 < obs_hom2) {
202     obs_homc = obs_hom2;
203     obs_homr = obs_hom1;
204   } else {
205     obs_homc = obs_hom1;
206     obs_homr = obs_hom2;
207   }
208   int64_t rare_copies = 2LL * obs_homr + obs_hets;
209   int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
210   int32_t tie_ct = 1;
211   double curr_hets_t2 = obs_hets;
212   double curr_homr_t2 = obs_homr;
213   double curr_homc_t2 = obs_homc;
214   double tailp = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
215   double centerp = 0;
216   double lastp2 = tailp;
217   double lastp1 = tailp;
218   double curr_hets_t1;
219   double curr_homr_t1;
220   double curr_homc_t1;
221   double preaddp;
222   if (!genotypes2) {
223     if (midp) {
224       return 0.5;
225     } else {
226       return 1;
227     }
228   }
229 
230   if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
231     // tail 1 = upper
232     while (curr_hets_t2 > 1.5) {
233       // het_probs[curr_hets] = 1
234       // het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0)
235       curr_homr_t2 += 1;
236       curr_homc_t2 += 1;
237       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
238       curr_hets_t2 -= 2;
239       if (lastp2 < EXACT_TEST_BIAS) {
240 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
241 	  tie_ct++;
242 	}
243 	tailp += lastp2;
244 	break;
245       }
246       centerp += lastp2;
247       if (centerp == INFINITY) {
248 	return 0;
249       }
250     }
251     if ((centerp == 0) && (!midp)) {
252       return 1;
253     }
254     while (curr_hets_t2 > 1.5) {
255       curr_homr_t2 += 1;
256       curr_homc_t2 += 1;
257       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
258       curr_hets_t2 -= 2;
259       preaddp = tailp;
260       tailp += lastp2;
261       if (tailp <= preaddp) {
262 	break;
263       }
264     }
265     curr_hets_t1 = obs_hets + 2;
266     curr_homr_t1 = obs_homr;
267     curr_homc_t1 = obs_homc;
268     while (curr_homr_t1 > 0.5) {
269       // het_probs[curr_hets + 2] = het_probs[curr_hets] * 4 * curr_homr * curr_homc / ((curr_hets + 2) * (curr_hets + 1))
270       lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
271       preaddp = tailp;
272       tailp += lastp1;
273       if (tailp <= preaddp) {
274 	break;
275       }
276       curr_hets_t1 += 2;
277       curr_homr_t1 -= 1;
278       curr_homc_t1 -= 1;
279     }
280   } else {
281     // tail 1 = lower
282     while (curr_homr_t2 > 0.5) {
283       curr_hets_t2 += 2;
284       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
285       curr_homr_t2 -= 1;
286       curr_homc_t2 -= 1;
287       if (lastp2 < EXACT_TEST_BIAS) {
288 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
289           tie_ct++;
290 	}
291 	tailp += lastp2;
292 	break;
293       }
294       centerp += lastp2;
295       if (centerp == INFINITY) {
296 	return 0;
297       }
298     }
299     if ((centerp == 0) && (!midp)) {
300       return 1;
301     }
302     while (curr_homr_t2 > 0.5) {
303       curr_hets_t2 += 2;
304       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
305       curr_homr_t2 -= 1;
306       curr_homc_t2 -= 1;
307       preaddp = tailp;
308       tailp += lastp2;
309       if (tailp <= preaddp) {
310 	break;
311       }
312     }
313     curr_hets_t1 = obs_hets;
314     curr_homr_t1 = obs_homr;
315     curr_homc_t1 = obs_homc;
316     while (curr_hets_t1 > 1.5) {
317       curr_homr_t1 += 1;
318       curr_homc_t1 += 1;
319       lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
320       preaddp = tailp;
321       tailp += lastp1;
322       if (tailp <= preaddp) {
323 	break;
324       }
325       curr_hets_t1 -= 2;
326     }
327   }
328   if (!midp) {
329     return tailp / (tailp + centerp);
330   } else {
331     return (tailp - ((1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (tailp + centerp);
332   }
333 }
334 
SNPHWE_t(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,double thresh)335 int32_t SNPHWE_t(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double thresh) {
336   // Threshold-test-only version of SNPHWE2() which is usually able to exit
337   // from the calculation earlier.  Returns 0 if these counts are close enough
338   // to Hardy-Weinberg equilibrium, 1 otherwise.
339   //
340   // Suppose, for definiteness, that the number of observed hets is no less
341   // than expectation.  (Same ideas apply for the other case.)  We proceed as
342   // follows:
343   // - Sum the *relative* likelihoods of more likely smaller het counts.
344   // - Determine the minimum tail mass to pass the threshold.
345   // - The majority of the time, the tail boundary elements are enough to pass
346   //   the threshold; we never need to sum the remainder of the tails.
347   // - And in the case of disequilibrium, we will often be able to immediately
348   //   determine that the tail sum cannot possibly pass the threshold, just by
349   //   looking at the tail boundary elements and using a geometric series to
350   //   upper-bound the tail sums.
351   // - Only when neither of these conditions hold do we start traveling down
352   //   the tails.
353   intptr_t obs_homc;
354   intptr_t obs_homr;
355   if (obs_hom1 < obs_hom2) {
356     obs_homc = obs_hom2;
357     obs_homr = obs_hom1;
358   } else {
359     obs_homc = obs_hom1;
360     obs_homr = obs_hom2;
361   }
362   int64_t rare_copies = 2LL * obs_homr + obs_hets;
363   int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
364   double curr_hets_t2 = obs_hets; // tail 2
365   double curr_homr_t2 = obs_homr;
366   double curr_homc_t2 = obs_homc;
367 
368   // Subtract epsilon from initial probability mass, so that we can compare to
369   // 1 when determining tail vs. center membership without floating point error
370   // biting us in the ass
371   double tailp1 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
372   double centerp = 0;
373   double lastp2 = tailp1;
374   double tailp2 = 0;
375   double tail1_ceil;
376   double tail2_ceil;
377   double lastp1;
378   double curr_hets_t1;
379   double curr_homr_t1;
380   double curr_homc_t1;
381 
382   // Initially, if center sum reaches this, the test can immediately fail.
383   // Once center is summed, this is recalculated, and when tail sum has reached
384   // this, we've passed.
385   double exit_thresh;
386   double exit_threshx;
387   double ratio;
388   double preaddp;
389   if (!genotypes2) {
390     return 0;
391   }
392 
393   // Convert thresh into reverse odds ratio.
394   thresh = (1 - thresh) / thresh;
395 
396   // Expected het count:
397   //   2 * rarefreq * (1 - rarefreq) * genotypes
398   // = 2 * (rare_copies / (2 * genotypes)) * (1 - rarefreq) * genotypes
399   // = rare_copies * (1 - (rare_copies / (2 * genotypes)))
400   // = (rare_copies * (2 * genotypes - rare_copies)) / (2 * genotypes)
401   //
402   // The computational identity is
403   //   P(nhets == n) := P(nhets == n+2) * (n+2) * (n+1) /
404   //                    (4 * homr(n) * homc(n))
405   // where homr() and homc() are the number of homozygous rares/commons needed
406   // to maintain the same allele frequencies.
407   // This probability is always decreasing when proceeding away from the
408   // expected het count.
409 
410   if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
411     // tail 1 = upper
412     if (obs_hets < 2) {
413       return 0;
414     }
415 
416     // An initial upper bound on the tail sum is useful, since it lets us
417     // report test failure before summing the entire center.  We use the
418     // trivial bound of 1 + floor(rare_copies / 2): that's the total number
419     // of possible het counts, and the relative probability for each count must
420     // be <= 1 if it's in the tail.
421     exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
422 
423     // het_probs[curr_hets] = 1
424     // het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1) / (4 * (curr_homr + 1) * (curr_homc + 1))
425     do {
426       curr_homr_t2 += 1;
427       curr_homc_t2 += 1;
428       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
429       curr_hets_t2 -= 2;
430       if (lastp2 < EXACT_TEST_BIAS) {
431 	tailp2 = lastp2;
432 	break;
433       }
434       centerp += lastp2;
435       if (centerp > exit_thresh) {
436 	return 1;
437       }
438     } while (curr_hets_t2 > 1.5);
439     exit_thresh = centerp / thresh;
440     if (tailp1 + tailp2 >= exit_thresh) {
441       return 0;
442     }
443     // c + cr + cr^2 + ... = c/(1-r), which is an upper bound for the tail sum
444     ratio = (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * (curr_homr_t2 + 1) * (curr_homc_t2 + 1));
445     tail2_ceil = tailp2 / (1 - ratio);
446     curr_hets_t1 = obs_hets + 2;
447     curr_homr_t1 = obs_homr;
448     curr_homc_t1 = obs_homc;
449     // ratio for the other tail
450     lastp1 = (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
451     tail1_ceil = tailp1 / (1 - lastp1);
452     if (tail1_ceil + tail2_ceil < exit_thresh) {
453       return 1;
454     }
455     lastp1 *= tailp1;
456     tailp1 += lastp1;
457 
458     if (obs_homr > 1) {
459       // het_probs[curr_hets + 2] = het_probs[curr_hets] * 4 * curr_homr * curr_homc / ((curr_hets + 2) * (curr_hets + 1))
460       exit_threshx = exit_thresh - tailp2;
461       do {
462 	curr_hets_t1 += 2;
463 	curr_homr_t1 -= 1;
464 	curr_homc_t1 -= 1;
465 	lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
466 	preaddp = tailp1;
467 	tailp1 += lastp1;
468 	if (tailp1 > exit_threshx) {
469 	  return 0;
470 	}
471 	if (tailp1 <= preaddp) {
472 	  break;
473 	}
474       } while (curr_homr_t1 > 1.5);
475     }
476     if (tailp1 + tail2_ceil < exit_thresh) {
477       return 1;
478     }
479     exit_threshx = exit_thresh - tailp1;
480     while (curr_hets_t2 > 1) {
481       curr_homr_t2 += 1;
482       curr_homc_t2 += 1;
483       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
484       preaddp = tailp2;
485       tailp2 += lastp2;
486       if (tailp2 >= exit_threshx) {
487 	return 0;
488       }
489       if (tailp2 <= preaddp) {
490 	return 1;
491       }
492       curr_hets_t2 -= 2;
493     }
494     return 1;
495   } else {
496     // tail 1 = lower
497     if (!obs_homr) {
498       return 0;
499     }
500     exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
501     do {
502       curr_hets_t2 += 2;
503       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
504       curr_homr_t2 -= 1;
505       curr_homc_t2 -= 1;
506       if (lastp2 < EXACT_TEST_BIAS) {
507 	tailp2 = lastp2;
508 	break;
509       }
510       centerp += lastp2;
511       if (centerp > exit_thresh) {
512 	return 1;
513       }
514     } while (curr_homr_t2 > 0.5);
515     exit_thresh = centerp / thresh;
516     if (tailp1 + tailp2 >= exit_thresh) {
517       return 0;
518     }
519     ratio = (4 * curr_homr_t2 * curr_homc_t2) / ((curr_hets_t2 + 2) * (curr_hets_t2 + 1));
520     tail2_ceil = tailp2 / (1 - ratio);
521     curr_hets_t1 = obs_hets;
522     curr_homr_t1 = obs_homr + 1;
523     curr_homc_t1 = obs_homc + 1;
524     lastp1 = (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
525     tail1_ceil = tailp1 / (1 - lastp1);
526     lastp1 *= tailp1;
527     tailp1 += lastp1;
528 
529     if (tail1_ceil + tail2_ceil < exit_thresh) {
530       return 1;
531     }
532     if (obs_hets >= 4) {
533       exit_threshx = exit_thresh - tailp2;
534       do {
535 	curr_hets_t1 -= 2;
536 	curr_homr_t1 += 1;
537 	curr_homc_t1 += 1;
538 	lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
539 	preaddp = tailp1;
540         tailp1 += lastp1;
541 	if (tailp1 > exit_threshx) {
542 	  return 0;
543 	}
544 	if (tailp1 <= preaddp) {
545 	  break;
546 	}
547       } while (curr_hets_t1 > 3.5);
548     }
549     if (tailp1 + tail2_ceil < exit_thresh) {
550       return 1;
551     }
552     exit_threshx = exit_thresh - tailp1;
553     while (curr_homr_t2 > 0.5) {
554       curr_hets_t2 += 2;
555       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
556       curr_homr_t2 -= 1;
557       curr_homc_t2 -= 1;
558       preaddp = tailp2;
559       tailp2 += lastp2;
560       if (tailp2 >= exit_threshx) {
561 	return 0;
562       }
563       if (tailp2 <= preaddp) {
564 	return 1;
565       }
566     }
567     return 1;
568   }
569 }
570 
SNPHWE_midp_t(int32_t obs_hets,int32_t obs_hom1,int32_t obs_hom2,double thresh)571 int32_t SNPHWE_midp_t(int32_t obs_hets, int32_t obs_hom1, int32_t obs_hom2, double thresh) {
572   // Mid-p version of SNPHWE_t().  (There are enough fiddly differences that I
573   // think it's better for this to be a separate function.)  Assumes threshold
574   // is smaller than 0.5.
575   intptr_t obs_homc;
576   intptr_t obs_homr;
577   if (obs_hom1 < obs_hom2) {
578     obs_homc = obs_hom2;
579     obs_homr = obs_hom1;
580   } else {
581     obs_homc = obs_hom1;
582     obs_homr = obs_hom2;
583   }
584   int64_t rare_copies = 2LL * obs_homr + obs_hets;
585   int64_t genotypes2 = (obs_hets + obs_homc + obs_homr) * 2LL;
586   double curr_hets_t2 = obs_hets; // tail 2
587   double curr_homr_t2 = obs_homr;
588   double curr_homc_t2 = obs_homc;
589   double tailp1 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5;
590   double centerp = tailp1;
591   double lastp2 = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
592   double tailp2 = 0;
593   double tail1_ceil;
594   double tail2_ceil;
595   double lastp1;
596   double curr_hets_t1;
597   double curr_homr_t1;
598   double curr_homc_t1;
599   double exit_thresh;
600   double exit_threshx;
601   double ratio;
602   double preaddp;
603   if (!genotypes2) {
604     return 0;
605   }
606   thresh = (1 - thresh) / thresh;
607   if (obs_hets * genotypes2 > rare_copies * (genotypes2 - rare_copies)) {
608     if (obs_hets < 2) {
609       return 0;
610     }
611     exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
612     do {
613       curr_homr_t2 += 1;
614       curr_homc_t2 += 1;
615       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
616       curr_hets_t2 -= 2;
617       if (lastp2 < EXACT_TEST_BIAS) {
618 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
619 	  // tie with original contingency table, apply mid-p correction here
620 	  // too
621           tailp2 = tailp1;
622           centerp += tailp1;
623 	} else {
624 	  tailp2 = lastp2;
625 	}
626 	break;
627       }
628       centerp += lastp2;
629       if (centerp > exit_thresh) {
630 	return 1;
631       }
632     } while (curr_hets_t2 > 1.5);
633     exit_thresh = centerp / thresh;
634     if (tailp1 + tailp2 >= exit_thresh) {
635       return 0;
636     }
637     ratio = (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * (curr_homr_t2 + 1) * (curr_homc_t2 + 1));
638     // this needs to work in both the tie and no-tie cases
639     tail2_ceil = tailp2 + lastp2 * ratio / (1 - ratio);
640     curr_hets_t1 = obs_hets + 2;
641     curr_homr_t1 = obs_homr;
642     curr_homc_t1 = obs_homc;
643     lastp1 = (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
644     // always a tie here
645     tail1_ceil = tailp1 * 2 / (1 - lastp1) - tailp1;
646     if (tail1_ceil + tail2_ceil < exit_thresh) {
647       return 1;
648     }
649     lastp1 *= tailp1 * 2;
650     tailp1 += lastp1;
651 
652     if (obs_homr > 1) {
653       exit_threshx = exit_thresh - tailp2;
654       do {
655 	curr_hets_t1 += 2;
656 	curr_homr_t1 -= 1;
657 	curr_homc_t1 -= 1;
658 	lastp1 *= (4 * curr_homr_t1 * curr_homc_t1) / (curr_hets_t1 * (curr_hets_t1 - 1));
659 	preaddp = tailp1;
660 	tailp1 += lastp1;
661 	if (tailp1 > exit_threshx) {
662 	  return 0;
663 	}
664 	if (tailp1 <= preaddp) {
665 	  break;
666 	}
667       } while (curr_homr_t1 > 1.5);
668     }
669     if (tailp1 + tail2_ceil < exit_thresh) {
670       return 1;
671     }
672     exit_threshx = exit_thresh - tailp1;
673     while (curr_hets_t2 > 1) {
674       curr_homr_t2 += 1;
675       curr_homc_t2 += 1;
676       lastp2 *= (curr_hets_t2 * (curr_hets_t2 - 1)) / (4 * curr_homr_t2 * curr_homc_t2);
677       preaddp = tailp2;
678       tailp2 += lastp2;
679       if (tailp2 >= exit_threshx) {
680 	return 0;
681       }
682       if (tailp2 <= preaddp) {
683 	return 1;
684       }
685       curr_hets_t2 -= 2;
686     }
687     return 1;
688   } else {
689     if (!obs_homr) {
690       return 0;
691     }
692     exit_thresh = (1 + (rare_copies / 2)) * thresh * EXACT_TEST_BIAS;
693     do {
694       curr_hets_t2 += 2;
695       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
696       curr_homr_t2 -= 1;
697       curr_homc_t2 -= 1;
698       if (lastp2 < EXACT_TEST_BIAS) {
699 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
700           tailp2 = tailp1;
701           centerp += tailp1;
702 	} else {
703 	  tailp2 = lastp2;
704 	}
705 	break;
706       }
707       centerp += lastp2;
708       if (centerp > exit_thresh) {
709 	return 1;
710       }
711     } while (curr_homr_t2 > 0.5);
712     exit_thresh = centerp / thresh;
713     if (tailp1 + tailp2 >= exit_thresh) {
714       return 0;
715     }
716     ratio = (4 * curr_homr_t2 * curr_homc_t2) / ((curr_hets_t2 + 2) * (curr_hets_t2 + 1));
717     tail2_ceil = tailp2 + lastp2 * ratio / (1 - ratio);
718     curr_hets_t1 = obs_hets;
719     curr_homr_t1 = obs_homr + 1;
720     curr_homc_t1 = obs_homc + 1;
721     lastp1 = (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
722     tail1_ceil = 2 * tailp1 / (1 - lastp1) - tailp1;
723     lastp1 *= 2 * tailp1;
724     tailp1 += lastp1;
725 
726     if (tail1_ceil + tail2_ceil < exit_thresh) {
727       return 1;
728     }
729     if (obs_hets >= 4) {
730       exit_threshx = exit_thresh - tailp2;
731       do {
732 	curr_hets_t1 -= 2;
733 	curr_homr_t1 += 1;
734 	curr_homc_t1 += 1;
735 	lastp1 *= (curr_hets_t1 * (curr_hets_t1 - 1)) / (4 * curr_homr_t1 * curr_homc_t1);
736 	preaddp = tailp1;
737         tailp1 += lastp1;
738 	if (tailp1 > exit_threshx) {
739 	  return 0;
740 	}
741 	if (tailp1 <= preaddp) {
742 	  break;
743 	}
744       } while (curr_hets_t1 > 3.5);
745     }
746     if (tailp1 + tail2_ceil < exit_thresh) {
747       return 1;
748     }
749     exit_threshx = exit_thresh - tailp1;
750     while (curr_homr_t2 > 0.5) {
751       curr_hets_t2 += 2;
752       lastp2 *= (4 * curr_homr_t2 * curr_homc_t2) / (curr_hets_t2 * (curr_hets_t2 - 1));
753       curr_homr_t2 -= 1;
754       curr_homc_t2 -= 1;
755       preaddp = tailp2;
756       tailp2 += lastp2;
757       if (tailp2 >= exit_threshx) {
758 	return 0;
759       }
760       if (tailp2 <= preaddp) {
761 	return 1;
762       }
763     }
764     return 1;
765   }
766 }
767 
fisher22(uint32_t m11,uint32_t m12,uint32_t m21,uint32_t m22,uint32_t midp)768 double fisher22(uint32_t m11, uint32_t m12, uint32_t m21, uint32_t m22, uint32_t midp) {
769   // Basic 2x2 Fisher exact test p-value calculation.
770   double tprob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
771   double cur_prob = tprob;
772   double cprob = 0;
773   int32_t tie_ct = 1;
774   uint32_t uii;
775   double cur11;
776   double cur12;
777   double cur21;
778   double cur22;
779   double preaddp;
780   // Ensure we are left of the distribution center, m11 <= m22, and m12 <= m21.
781   if (m12 > m21) {
782     uii = m12;
783     m12 = m21;
784     m21 = uii;
785   }
786   if (m11 > m22) {
787     uii = m11;
788     m11 = m22;
789     m22 = uii;
790   }
791   if ((((uint64_t)m11) * m22) > (((uint64_t)m12) * m21)) {
792     uii = m11;
793     m11 = m12;
794     m12 = uii;
795     uii = m21;
796     m21 = m22;
797     m22 = uii;
798   }
799   cur11 = m11;
800   cur12 = m12;
801   cur21 = m21;
802   cur22 = m22;
803   while (cur12 > 0.5) {
804     cur11 += 1;
805     cur22 += 1;
806     cur_prob *= (cur12 * cur21) / (cur11 * cur22);
807     cur12 -= 1;
808     cur21 -= 1;
809     if (cur_prob == INFINITY) {
810       return 0;
811     }
812     if (cur_prob < EXACT_TEST_BIAS) {
813       if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
814         tie_ct++;
815       }
816       tprob += cur_prob;
817       break;
818     }
819     cprob += cur_prob;
820   }
821   if ((cprob == 0) && (!midp)) {
822     return 1;
823   }
824   while (cur12 > 0.5) {
825     cur11 += 1;
826     cur22 += 1;
827     cur_prob *= (cur12 * cur21) / (cur11 * cur22);
828     cur12 -= 1;
829     cur21 -= 1;
830     preaddp = tprob;
831     tprob += cur_prob;
832     if (tprob <= preaddp) {
833       break;
834     }
835   }
836   if (m11) {
837     cur11 = m11;
838     cur12 = m12;
839     cur21 = m21;
840     cur22 = m22;
841     cur_prob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
842     do {
843       cur12 += 1;
844       cur21 += 1;
845       cur_prob *= (cur11 * cur22) / (cur12 * cur21);
846       cur11 -= 1;
847       cur22 -= 1;
848       preaddp = tprob;
849       tprob += cur_prob;
850       if (tprob <= preaddp) {
851         if (!midp) {
852 	  return preaddp / (cprob + preaddp);
853 	} else {
854           return (preaddp - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (cprob + preaddp);
855 	}
856       }
857     } while (cur11 > 0.5);
858   }
859   if (!midp) {
860     return tprob / (cprob + tprob);
861   } else {
862     return (tprob - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (cprob + tprob);
863   }
864 }
865 
fisher22_tail_pval(uint32_t m11,uint32_t m12,uint32_t m21,uint32_t m22,int32_t right_offset,double tot_prob_recip,double right_prob,uint32_t midp,uint32_t new_m11)866 double fisher22_tail_pval(uint32_t m11, uint32_t m12, uint32_t m21, uint32_t m22, int32_t right_offset, double tot_prob_recip, double right_prob, uint32_t midp, uint32_t new_m11) {
867   // Given that the left (w.r.t. m11) reference contingency table has
868   // likelihood 1/tot_prob, the contingency table with m11 increased by
869   // right_offset has likelihood right_prob/tot_prob, and the tails (up to but
870   // not including the two references) sum to tail_sum/tot_prob, this
871   // calculates the p-value of the given m11 (which must be on one tail).
872   double left_prob = 1.0;
873   double dxx = ((intptr_t)new_m11);
874   double cur11;
875   double cur12;
876   double cur21;
877   double cur22;
878   double psum;
879   double thresh;
880   if (new_m11 < m11) {
881     cur11 = ((intptr_t)m11);
882     cur12 = ((intptr_t)m12);
883     cur21 = ((intptr_t)m21);
884     cur22 = ((intptr_t)m22);
885     dxx += 0.5; // unnecessary (53 vs. 32 bits precision), but whatever
886     do {
887       cur12 += 1;
888       cur21 += 1;
889       left_prob *= cur11 * cur22 / (cur12 * cur21);
890       cur11 -= 1;
891       cur22 -= 1;
892     } while (cur11 > dxx);
893     if (left_prob == 0) {
894       return 0;
895     }
896     if (!midp) {
897       psum = left_prob;
898     } else {
899       psum = left_prob * 0.5;
900     }
901     thresh = left_prob * (1 + FISHER_EPSILON);
902     do {
903       if (cur11 < 0.5) {
904 	break;
905       }
906       cur12 += 1;
907       cur21 += 1;
908       left_prob *= cur11 * cur22 / (cur12 * cur21);
909       cur11 -= 1;
910       cur22 -= 1;
911       dxx = psum;
912       psum += left_prob;
913     } while (psum > dxx);
914     cur11 = ((intptr_t)(m11 + right_offset));
915     cur12 = ((intptr_t)(m12 - right_offset));
916     cur21 = ((intptr_t)(m21 - right_offset));
917     cur22 = ((intptr_t)(m22 + right_offset));
918     while (right_prob > thresh) {
919       cur11 += 1;
920       cur22 += 1;
921       right_prob *= cur12 * cur21 / (cur11 * cur22);
922       cur12 -= 1;
923       cur21 -= 1;
924     }
925     if (right_prob > 0) {
926       if (midp && (right_prob < thresh * (1 - 2 * FISHER_EPSILON))) {
927 	psum += right_prob * 0.5;
928       } else {
929         psum += right_prob;
930       }
931       do {
932 	cur11 += 1;
933 	cur22 += 1;
934 	right_prob *= cur12 * cur21 / (cur11 * cur22);
935 	cur12 -= 1;
936 	cur21 -= 1;
937 	dxx = psum;
938 	psum += right_prob;
939       } while (psum > dxx);
940     }
941   } else {
942     dxx -= 0.5;
943     cur11 = ((intptr_t)(m11 + right_offset));
944     cur12 = ((intptr_t)(m12 - right_offset));
945     cur21 = ((intptr_t)(m21 - right_offset));
946     cur22 = ((intptr_t)(m22 + right_offset));
947     do {
948       cur11 += 1;
949       cur22 += 1;
950       right_prob *= cur12 * cur21 / (cur11 * cur22);
951       cur12 -= 1;
952       cur21 -= 1;
953     } while (cur11 < dxx);
954     if (right_prob == 0) {
955       return 0;
956     }
957     if (!midp) {
958       psum = right_prob;
959     } else {
960       psum = right_prob * 0.5;
961     }
962     thresh = right_prob * (1 + FISHER_EPSILON);
963     do {
964       if (cur12 < 0.5) {
965 	break;
966       }
967       cur11 += 1;
968       cur22 += 1;
969       right_prob *= cur12 * cur21 / (cur11 * cur22);
970       cur12 -= 1;
971       cur21 -= 1;
972       dxx = psum;
973       psum += right_prob;
974     } while (psum > dxx);
975     cur11 = ((intptr_t)m11);
976     cur12 = ((intptr_t)m12);
977     cur21 = ((intptr_t)m21);
978     cur22 = ((intptr_t)m22);
979     while (left_prob > thresh) {
980       cur12 += 1;
981       cur21 += 1;
982       left_prob *= cur11 * cur22 / (cur12 * cur21);
983       cur11 -= 1;
984       cur22 -= 1;
985     }
986     if (left_prob > 0) {
987       if (midp && (left_prob < thresh * (1 - 2 * FISHER_EPSILON))) {
988 	psum += left_prob * 0.5;
989       } else {
990         psum += left_prob;
991       }
992       do {
993 	cur12 += 1;
994 	cur21 += 1;
995 	left_prob *= cur11 * cur22 / (cur12 * cur21);
996 	cur11 -= 1;
997 	cur22 -= 1;
998 	dxx = psum;
999 	psum += left_prob;
1000       } while (psum > dxx);
1001     }
1002   }
1003   return psum * tot_prob_recip;
1004 }
1005 
fisher22_precomp_pval_bounds(double pval,uint32_t midp,uint32_t row1_sum,uint32_t col1_sum,uint32_t total,uint32_t * bounds,double * tprobs)1006 void fisher22_precomp_pval_bounds(double pval, uint32_t midp, uint32_t row1_sum, uint32_t col1_sum, uint32_t total, uint32_t* bounds, double* tprobs) {
1007   // bounds[0] = m11 min
1008   // bounds[1] = m11 (max + 1)
1009   // bounds[2] = m11 min after including ties
1010   // bounds[3] = m11 (max + 1) after including ties
1011   // Treating m11 as the only variable, this returns the minimum and (maximum +
1012   // 1) values of m11 which are less extreme than the observed result, and
1013   // notes ties (2^{-40} tolerance).  Also, returns the values necessary for
1014   // invoking fisher22_tail_pval() with
1015   //   m11 := bounds[2] and
1016   //   right_offset := bounds[3] - bounds[2] - 1
1017   // in tprobs[0], [1], and [2] (if tprobs is not NULL).
1018   //
1019   // Algorithm:
1020   // 1. Determine center.
1021   // 2. Sum unscaled probabilities in both directions to precision limit.
1022   // 3. Proceed further outwards to (pval * unscaled_psum) precision limit,
1023   //    fill in the remaining return values.
1024   double tot_prob = 1.0 / EXACT_TEST_BIAS;
1025   double left_prob = tot_prob;
1026   double right_prob = tot_prob;
1027   intptr_t m11_offset = 0;
1028   double tail_prob = 0;
1029   double cmult = midp? 0.5 : 1.0;
1030   double dxx;
1031   double left11;
1032   double left12;
1033   double left21;
1034   double left22;
1035   double right11;
1036   double right12;
1037   double right21;
1038   double right22;
1039   double cur_prob;
1040   double cur11;
1041   double cur12;
1042   double cur21;
1043   double cur22;
1044   double threshold;
1045   intptr_t lii;
1046   uint32_t uii;
1047   if (!total) {
1048     // hardcode this to avoid divide-by-zero
1049     bounds[0] = 0;
1050     bounds[1] = 0;
1051     bounds[2] = 0;
1052     bounds[3] = 1;
1053     // no need to initialize the other return values, they should never be used
1054     return;
1055   } else {
1056     if (pval == 0) {
1057       if (total >= row1_sum + col1_sum) {
1058 	bounds[0] = 0;
1059 	bounds[1] = MINV(row1_sum, col1_sum) + 1;
1060       } else {
1061 	bounds[0] = row1_sum + col1_sum - total;
1062 	bounds[1] = total - MAXV(row1_sum, col1_sum) + 1;
1063       }
1064       bounds[2] = bounds[0];
1065       bounds[3] = bounds[1];
1066       return;
1067     }
1068   }
1069   // Center must be adjacent to the x which satisfies
1070   //   (m11 + x)(m22 + x) = (m12 - x)(m21 - x), so
1071   //   x = (m12 * m21 - m11 * m22) / (m11 + m12 + m21 + m22)
1072   if (total >= row1_sum + col1_sum) {
1073     // m11 = 0;
1074     // m12 = row1_sum;
1075     // m21 = col1_sum;
1076     // m22 = total - row1_sum - col1_sum;
1077     lii = (((uint64_t)row1_sum) * col1_sum) / total;
1078     left11 = lii;
1079     left12 = row1_sum - lii;
1080     left21 = col1_sum - lii;
1081     left22 = (total - row1_sum - col1_sum) + lii;
1082   } else {
1083     // m11 = row1_sum + col1_sum - total;
1084     // m12 = row1_sum - m11;
1085     // m21 = col1_sum - m11;
1086     // m22 = 0;
1087     lii = (((uint64_t)(total - row1_sum)) * (total - col1_sum)) / total;
1088     // Force m11 <= m22 for internal calculation, then adjust at end.
1089     m11_offset = row1_sum + col1_sum - total;
1090     left11 = lii;
1091     left12 = total - col1_sum - lii;
1092     left21 = total - row1_sum - lii;
1093     left22 = m11_offset + lii;
1094   }
1095   // We rounded x down.  Should we have rounded up instead?
1096   if ((left11 + 1) * (left22 + 1) < left12 * left21) {
1097     left11 += 1;
1098     left12 -= 1;
1099     left21 -= 1;
1100     left22 += 1;
1101   }
1102   // Safe to force m12 <= m21.
1103   if (left12 > left21) {
1104     dxx = left12;
1105     left12 = left21;
1106     left21 = dxx;
1107   }
1108   // Sum right side to limit, then left.
1109   right11 = left11;
1110   right12 = left12;
1111   right21 = left21;
1112   right22 = left22;
1113   do {
1114     if (right12 < 0.5) {
1115       break;
1116     }
1117     right11 += 1;
1118     right22 += 1;
1119     right_prob *= (right12 * right21) / (right11 * right22);
1120     right12 -= 1;
1121     right21 -= 1;
1122     dxx = tot_prob;
1123     tot_prob += right_prob;
1124   } while (tot_prob > dxx);
1125   do {
1126     if (left11 < 0.5) {
1127       break;
1128     }
1129     left12 += 1;
1130     left21 += 1;
1131     left_prob *= (left11 * left22) / (left12 * left21);
1132     left11 -= 1;
1133     left22 -= 1;
1134     dxx = tot_prob;
1135     tot_prob += left_prob;
1136   } while (tot_prob > dxx);
1137   // Now traverse the tails to p-value precision limit.
1138   // Upper bound for tail sum, if current element c is included:
1139   //   (c + cr + cr^2 + ...) + (c + cs + cs^2 + ...)
1140   // = c(1/(1 - r) + 1/(1 - s))
1141   // Compare this to pval * tot_prob.
1142   // I.e. compare c to pval * tot_prob * (1-r)(1-s) / (2-r-s)
1143   dxx = 1 - (left11 * left22) / ((left12 + 1) * (left21 + 1));
1144   threshold = 1 - (right12 * right21) / ((right11 + 1) * (right22 + 1));
1145   threshold = pval * tot_prob * dxx * threshold / (dxx + threshold);
1146   while (left11 > 0.5) {
1147     if (left_prob < threshold) {
1148       tail_prob = left_prob;
1149       cur11 = left11;
1150       cur12 = left12;
1151       cur21 = left21;
1152       cur22 = left22;
1153       cur_prob = left_prob;
1154       do {
1155 	cur12 += 1;
1156 	cur21 += 1;
1157 	cur_prob *= (cur11 * cur22) / (cur12 * cur21);
1158 	cur11 -= 1;
1159 	cur22 -= 1;
1160 	dxx = tail_prob;
1161 	tail_prob += cur_prob;
1162       } while (dxx < tail_prob);
1163       left11 += 1;
1164       left22 += 1;
1165       left_prob *= (left12 * left21) / (left11 * left22);
1166       left12 -= 1;
1167       left21 -= 1;
1168       break;
1169     }
1170     left12 += 1;
1171     left21 += 1;
1172     left_prob *= (left11 * left22) / (left12 * left21);
1173     left11 -= 1;
1174     left22 -= 1;
1175   }
1176   while (right12 > 0.5) {
1177     if (right_prob < threshold) {
1178       tail_prob += right_prob;
1179       cur11 = right11;
1180       cur12 = right12;
1181       cur21 = right21;
1182       cur22 = right22;
1183       cur_prob = right_prob;
1184       do {
1185 	cur11 += 1;
1186 	cur22 += 1;
1187 	cur_prob *= (cur12 * cur21) / (cur11 * cur22);
1188 	cur12 -= 1;
1189 	cur21 -= 1;
1190 	dxx = tail_prob;
1191 	tail_prob += cur_prob;
1192       } while (dxx < tail_prob);
1193       right12 += 1;
1194       right21 += 1;
1195       right_prob *= (right11 * right22) / (right12 * right21);
1196       right11 -= 1;
1197       right22 -= 1;
1198       break;
1199     }
1200     right11 += 1;
1201     right22 += 1;
1202     right_prob *= (right12 * right21) / (right11 * right22);
1203     right12 -= 1;
1204     right21 -= 1;
1205   }
1206   dxx = pval * tot_prob * (1 - FISHER_EPSILON / 2);
1207   threshold = pval * tot_prob * (1 + FISHER_EPSILON / 2);
1208   lii = 0;
1209   while (1) {
1210     if (left_prob < right_prob * (1 - FISHER_EPSILON / 2)) {
1211       cur_prob = tail_prob + left_prob * cmult;
1212       if (cur_prob > threshold) {
1213 	break;
1214       }
1215       tail_prob += left_prob;
1216       uii = 1;
1217     } else if (right_prob < left_prob * (1 - FISHER_EPSILON / 2)) {
1218       cur_prob = tail_prob + right_prob * cmult;
1219       if (cur_prob > threshold) {
1220 	break;
1221       }
1222       tail_prob += right_prob;
1223       uii = 2;
1224     } else {
1225       cur_prob = tail_prob + (left_prob + right_prob) * cmult;
1226       if (cur_prob > threshold) {
1227 	if (left11 == right11) {
1228 	  cur_prob = tail_prob + left_prob * cmult;
1229 	  // center: left and right refer to same table.  subcases:
1230 	  // 1. cur_prob > threshold: center table has less extreme pval.
1231 	  //    lii = 0, both intervals size 1
1232 	  // 2. dxx < cur_prob < threshold: center table has equal pval.
1233 	  //    lii = 1, less-than interval size 0 but leq interval size 1
1234 	  // 3. cur_prob < dxx: even centermost table has more extreme pval
1235 	  //    (only possible with mid-p adj).
1236 	  //    lii = 0, we increment left11 so both intervals size 0
1237 	  if (cur_prob < threshold) {
1238 	    if (cur_prob > dxx) {
1239 	      lii = 1;
1240 	    } else {
1241 	      left11++;
1242 	      left22++;
1243 	      left_prob *= (left12 * left21) / (left11 * left22);
1244 	    }
1245 	  }
1246 	}
1247 	break;
1248       }
1249       tail_prob += left_prob + right_prob;
1250       uii = 3;
1251     }
1252     if (cur_prob > dxx) {
1253       lii = uii;
1254       break;
1255     }
1256     // if more speed is necessary, we could use a buffer to save all unscaled
1257     // probabilities during the initial outward traversal.
1258     if (uii & 1) {
1259       left11 += 1;
1260       left22 += 1;
1261       left_prob *= (left12 * left21) / (left11 * left22);
1262       left12 -= 1;
1263       left21 -= 1;
1264     }
1265     if (uii & 2) {
1266       right12 += 1;
1267       right21 += 1;
1268       right_prob *= (right11 * right22) / (right12 * right21);
1269       right11 -= 1;
1270       right22 -= 1;
1271     }
1272   }
1273   bounds[2] = m11_offset + ((intptr_t)left11);
1274   bounds[3] = m11_offset + ((intptr_t)right11) + 1;
1275   bounds[0] = bounds[2] + (lii & 1);
1276   bounds[1] = bounds[3] - (lii >> 1);
1277   if (!tprobs) {
1278     return;
1279   }
1280   dxx = 1.0 / left_prob;
1281   tprobs[0] = left_prob / tot_prob;
1282   tprobs[1] = right_prob * dxx;
1283   /*
1284   if (lii & 1) {
1285     tail_prob -= left_prob;
1286   }
1287   if (lii >> 1) {
1288     tail_prob -= right_prob;
1289   }
1290   tprobs[2] = tail_prob * dxx;
1291   */
1292 }
1293 
fisher23_tailsum(double * base_probp,double * saved12p,double * saved13p,double * saved22p,double * saved23p,double * totalp,uint32_t * tie_ctp,uint32_t right_side)1294 int32_t fisher23_tailsum(double* base_probp, double* saved12p, double* saved13p, double* saved22p, double* saved23p, double *totalp, uint32_t* tie_ctp, uint32_t right_side) {
1295   double total = 0;
1296   double cur_prob = *base_probp;
1297   double tmp12 = *saved12p;
1298   double tmp13 = *saved13p;
1299   double tmp22 = *saved22p;
1300   double tmp23 = *saved23p;
1301   double tmps12;
1302   double tmps13;
1303   double tmps22;
1304   double tmps23;
1305   double prev_prob;
1306   // identify beginning of tail
1307   if (right_side) {
1308     if (cur_prob > EXACT_TEST_BIAS) {
1309       prev_prob = tmp13 * tmp22;
1310       while (prev_prob > 0.5) {
1311 	tmp12 += 1;
1312 	tmp23 += 1;
1313 	cur_prob *= prev_prob / (tmp12 * tmp23);
1314 	tmp13 -= 1;
1315 	tmp22 -= 1;
1316 	if (cur_prob <= EXACT_TEST_BIAS) {
1317 	  break;
1318 	}
1319 	prev_prob = tmp13 * tmp22;
1320       }
1321       *base_probp = cur_prob;
1322       tmps12 = tmp12;
1323       tmps13 = tmp13;
1324       tmps22 = tmp22;
1325       tmps23 = tmp23;
1326     } else {
1327       tmps12 = tmp12;
1328       tmps13 = tmp13;
1329       tmps22 = tmp22;
1330       tmps23 = tmp23;
1331       while (1) {
1332 	prev_prob = cur_prob;
1333 	tmp13 += 1;
1334 	tmp22 += 1;
1335 	cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1336 	if (cur_prob < prev_prob) {
1337 	  return 1;
1338 	}
1339 	tmp12 -= 1;
1340 	tmp23 -= 1;
1341 	if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1342 	  // throw in extra (1 - SMALL_EPSILON) multiplier to prevent rounding
1343 	  // errors from causing this to keep going when the left-side test
1344 	  // stopped
1345 	  if (cur_prob > (1 - SMALL_EPSILON) * EXACT_TEST_BIAS) {
1346 	    break;
1347 	  }
1348           *tie_ctp += 1;
1349 	}
1350 	total += cur_prob;
1351       }
1352       prev_prob = cur_prob;
1353       cur_prob = *base_probp;
1354       *base_probp = prev_prob;
1355     }
1356   } else {
1357     if (cur_prob > EXACT_TEST_BIAS) {
1358       prev_prob = tmp12 * tmp23;
1359       while (prev_prob > 0.5) {
1360 	tmp13 += 1;
1361 	tmp22 += 1;
1362 	cur_prob *= prev_prob / (tmp13 * tmp22);
1363 	tmp12 -= 1;
1364 	tmp23 -= 1;
1365 	if (cur_prob <= EXACT_TEST_BIAS) {
1366 	  break;
1367 	}
1368 	prev_prob = tmp12 * tmp23;
1369       }
1370       *base_probp = cur_prob;
1371       tmps12 = tmp12;
1372       tmps13 = tmp13;
1373       tmps22 = tmp22;
1374       tmps23 = tmp23;
1375     } else {
1376       tmps12 = tmp12;
1377       tmps13 = tmp13;
1378       tmps22 = tmp22;
1379       tmps23 = tmp23;
1380       while (1) {
1381 	prev_prob = cur_prob;
1382 	tmp12 += 1;
1383 	tmp23 += 1;
1384 	cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1385 	if (cur_prob < prev_prob) {
1386 	  return 1;
1387 	}
1388 	tmp13 -= 1;
1389 	tmp22 -= 1;
1390 	if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1391 	  if (cur_prob > EXACT_TEST_BIAS) {
1392 	    break;
1393 	  }
1394           *tie_ctp += 1;
1395 	}
1396 	total += cur_prob;
1397       }
1398       prev_prob = cur_prob;
1399       cur_prob = *base_probp;
1400       *base_probp = prev_prob;
1401     }
1402   }
1403   *saved12p = tmp12;
1404   *saved13p = tmp13;
1405   *saved22p = tmp22;
1406   *saved23p = tmp23;
1407   if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1408     if (cur_prob > EXACT_TEST_BIAS) {
1409       // even most extreme table on this side is too probable
1410       *totalp = 0;
1411       return 0;
1412     }
1413     *tie_ctp += 1;
1414   }
1415   // sum tail to floating point precision limit
1416   if (right_side) {
1417     prev_prob = total;
1418     total += cur_prob;
1419     while (total > prev_prob) {
1420       tmps12 += 1;
1421       tmps23 += 1;
1422       cur_prob *= (tmps13 * tmps22) / (tmps12 * tmps23);
1423       tmps13 -= 1;
1424       tmps22 -= 1;
1425       prev_prob = total;
1426       total += cur_prob;
1427     }
1428   } else {
1429     prev_prob = total;
1430     total += cur_prob;
1431     while (total > prev_prob) {
1432       tmps13 += 1;
1433       tmps22 += 1;
1434       cur_prob *= (tmps12 * tmps23) / (tmps13 * tmps22);
1435       tmps12 -= 1;
1436       tmps23 -= 1;
1437       prev_prob = total;
1438       total += cur_prob;
1439     }
1440   }
1441   *totalp = total;
1442   return 0;
1443 }
1444 
fisher23(uint32_t m11,uint32_t m12,uint32_t m13,uint32_t m21,uint32_t m22,uint32_t m23,uint32_t midp)1445 double fisher23(uint32_t m11, uint32_t m12, uint32_t m13, uint32_t m21, uint32_t m22, uint32_t m23, uint32_t midp) {
1446   // 2x3 Fisher-Freeman-Halton exact test p-value calculation.
1447   // The number of tables involved here is still small enough that the network
1448   // algorithm (and the improved variants thereof that I've seen) are
1449   // suboptimal; a 2-dimensional version of the SNPHWE2 strategy has higher
1450   // performance.
1451   // 2x4, 2x5, and 3x3 should also be practical with this method, but beyond
1452   // that I doubt it's worth the trouble.
1453   // Complexity of approach is O(n^{df/2}), where n is number of observations.
1454   double cur_prob = (1 - FISHER_EPSILON) * EXACT_TEST_BIAS;
1455   double tprob = cur_prob;
1456   double cprob = 0;
1457   double dyy = 0;
1458   uint32_t tie_ct = 1;
1459   uint32_t dir = 0; // 0 = forwards, 1 = backwards
1460   double base_probl;
1461   double base_probr;
1462   double orig_base_probl;
1463   double orig_base_probr;
1464   double orig_row_prob;
1465   double row_prob;
1466   uint32_t uii;
1467   uint32_t ujj;
1468   uint32_t ukk;
1469   double cur11;
1470   double cur21;
1471   double savedl12;
1472   double savedl13;
1473   double savedl22;
1474   double savedl23;
1475   double savedr12;
1476   double savedr13;
1477   double savedr22;
1478   double savedr23;
1479   double orig_savedl12;
1480   double orig_savedl13;
1481   double orig_savedl22;
1482   double orig_savedl23;
1483   double orig_savedr12;
1484   double orig_savedr13;
1485   double orig_savedr22;
1486   double orig_savedr23;
1487   double tmp12;
1488   double tmp13;
1489   double tmp22;
1490   double tmp23;
1491   double dxx;
1492   double preaddp;
1493   // Ensure m11 + m21 <= m12 + m22 <= m13 + m23.
1494   uii = m11 + m21;
1495   ujj = m12 + m22;
1496   if (uii > ujj) {
1497     ukk = m11;
1498     m11 = m12;
1499     m12 = ukk;
1500     ukk = m21;
1501     m21 = m22;
1502     m22 = ukk;
1503     ukk = uii;
1504     uii = ujj;
1505     ujj = ukk;
1506   }
1507   ukk = m13 + m23;
1508   if (ujj > ukk) {
1509     ujj = ukk;
1510     ukk = m12;
1511     m12 = m13;
1512     m13 = ukk;
1513     ukk = m22;
1514     m22 = m23;
1515     m23 = ukk;
1516   }
1517   if (uii > ujj) {
1518     ukk = m11;
1519     m11 = m12;
1520     m12 = ukk;
1521     ukk = m21;
1522     m21 = m22;
1523     m22 = ukk;
1524   }
1525   // Ensure majority of probability mass is in front of m11.
1526   if ((((uint64_t)m11) * (m22 + m23)) > (((uint64_t)m21) * (m12 + m13))) {
1527     ukk = m11;
1528     m11 = m21;
1529     m21 = ukk;
1530     ukk = m12;
1531     m12 = m22;
1532     m22 = ukk;
1533     ukk = m13;
1534     m13 = m23;
1535     m23 = ukk;
1536   }
1537   if ((((uint64_t)m12) * m23) > (((uint64_t)m13) * m22)) {
1538     base_probr = cur_prob;
1539     savedr12 = m12;
1540     savedr13 = m13;
1541     savedr22 = m22;
1542     savedr23 = m23;
1543     tmp12 = savedr12;
1544     tmp13 = savedr13;
1545     tmp22 = savedr22;
1546     tmp23 = savedr23;
1547     // m12 and m23 must be nonzero
1548     dxx = tmp12 * tmp23;
1549     do {
1550       tmp13 += 1;
1551       tmp22 += 1;
1552       cur_prob *= dxx / (tmp13 * tmp22);
1553       tmp12 -= 1;
1554       tmp23 -= 1;
1555       if (cur_prob <= EXACT_TEST_BIAS) {
1556 	if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1557           tie_ct++;
1558 	}
1559 	tprob += cur_prob;
1560 	break;
1561       }
1562       cprob += cur_prob;
1563       if (cprob == INFINITY) {
1564 	return 0;
1565       }
1566       dxx = tmp12 * tmp23;
1567       // must enforce tmp12 >= 0 and tmp23 >= 0 since we're saving these
1568     } while (dxx > 0.5);
1569     savedl12 = tmp12;
1570     savedl13 = tmp13;
1571     savedl22 = tmp22;
1572     savedl23 = tmp23;
1573     base_probl = cur_prob;
1574     do {
1575       tmp13 += 1;
1576       tmp22 += 1;
1577       cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1578       tmp12 -= 1;
1579       tmp23 -= 1;
1580       preaddp = tprob;
1581       tprob += cur_prob;
1582     } while (tprob > preaddp);
1583     tmp12 = savedr12;
1584     tmp13 = savedr13;
1585     tmp22 = savedr22;
1586     tmp23 = savedr23;
1587     cur_prob = base_probr;
1588     do {
1589       tmp12 += 1;
1590       tmp23 += 1;
1591       cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1592       tmp13 -= 1;
1593       tmp22 -= 1;
1594       preaddp = tprob;
1595       tprob += cur_prob;
1596     } while (tprob > preaddp);
1597   } else {
1598     base_probl = cur_prob;
1599     savedl12 = m12;
1600     savedl13 = m13;
1601     savedl22 = m22;
1602     savedl23 = m23;
1603     if (!((((uint64_t)m12) * m23) + (((uint64_t)m13) * m22))) {
1604       base_probr = cur_prob;
1605       savedr12 = savedl12;
1606       savedr13 = savedl13;
1607       savedr22 = savedl22;
1608       savedr23 = savedl23;
1609     } else {
1610       tmp12 = savedl12;
1611       tmp13 = savedl13;
1612       tmp22 = savedl22;
1613       tmp23 = savedl23;
1614       dxx = tmp13 * tmp22;
1615       do {
1616 	tmp12 += 1;
1617 	tmp23 += 1;
1618 	cur_prob *= dxx / (tmp12 * tmp23);
1619 	tmp13 -= 1;
1620 	tmp22 -= 1;
1621 	if (cur_prob <= EXACT_TEST_BIAS) {
1622           if (cur_prob > (1 - 2 * FISHER_EPSILON) * EXACT_TEST_BIAS) {
1623             tie_ct++;
1624 	  }
1625 	  tprob += cur_prob;
1626 	  break;
1627 	}
1628 	cprob += cur_prob;
1629 	if (cprob == INFINITY) {
1630 	  return 0;
1631 	}
1632 	dxx = tmp13 * tmp22;
1633       } while (dxx > 0.5);
1634       savedr12 = tmp12;
1635       savedr13 = tmp13;
1636       savedr22 = tmp22;
1637       savedr23 = tmp23;
1638       base_probr = cur_prob;
1639       do {
1640 	tmp12 += 1;
1641 	tmp23 += 1;
1642 	cur_prob *= (tmp13 * tmp22) / (tmp12 * tmp23);
1643 	tmp13 -= 1;
1644 	tmp22 -= 1;
1645 	preaddp = tprob;
1646 	tprob += cur_prob;
1647       } while (tprob > preaddp);
1648       tmp12 = savedl12;
1649       tmp13 = savedl13;
1650       tmp22 = savedl22;
1651       tmp23 = savedl23;
1652       cur_prob = base_probl;
1653       do {
1654 	tmp13 += 1;
1655 	tmp22 += 1;
1656 	cur_prob *= (tmp12 * tmp23) / (tmp13 * tmp22);
1657 	tmp12 -= 1;
1658 	tmp23 -= 1;
1659 	preaddp = tprob;
1660 	tprob += cur_prob;
1661       } while (tprob > preaddp);
1662     }
1663   }
1664   row_prob = tprob + cprob;
1665   orig_base_probl = base_probl;
1666   orig_base_probr = base_probr;
1667   orig_row_prob = row_prob;
1668   orig_savedl12 = savedl12;
1669   orig_savedl13 = savedl13;
1670   orig_savedl22 = savedl22;
1671   orig_savedl23 = savedl23;
1672   orig_savedr12 = savedr12;
1673   orig_savedr13 = savedr13;
1674   orig_savedr22 = savedr22;
1675   orig_savedr23 = savedr23;
1676   for (; dir < 2; dir++) {
1677     cur11 = m11;
1678     cur21 = m21;
1679     if (dir) {
1680       base_probl = orig_base_probl;
1681       base_probr = orig_base_probr;
1682       row_prob = orig_row_prob;
1683       savedl12 = orig_savedl12;
1684       savedl13 = orig_savedl13;
1685       savedl22 = orig_savedl22;
1686       savedl23 = orig_savedl23;
1687       savedr12 = orig_savedr12;
1688       savedr13 = orig_savedr13;
1689       savedr22 = orig_savedr22;
1690       savedr23 = orig_savedr23;
1691       ukk = m11;
1692       if (ukk > m22 + m23) {
1693 	ukk = m22 + m23;
1694       }
1695     } else {
1696       ukk = m21;
1697       if (ukk > m12 + m13) {
1698 	ukk = m12 + m13;
1699       }
1700     }
1701     ukk++;
1702     while (--ukk) {
1703       if (dir) {
1704 	cur21 += 1;
1705 	if (savedl23) {
1706 	  savedl13 += 1;
1707 	  row_prob *= (cur11 * (savedl22 + savedl23)) / (cur21 * (savedl12 + savedl13));
1708 	  base_probl *= (cur11 * savedl23) / (cur21 * savedl13);
1709 	  savedl23 -= 1;
1710 	} else {
1711 	  savedl12 += 1;
1712 	  row_prob *= (cur11 * (savedl22 + savedl23)) / (cur21 * (savedl12 + savedl13));
1713 	  base_probl *= (cur11 * savedl22) / (cur21 * savedl12);
1714 	  savedl22 -= 1;
1715 	}
1716 	cur11 -= 1;
1717       } else {
1718 	cur11 += 1;
1719 	if (savedl12) {
1720 	  savedl22 += 1;
1721 	  row_prob *= (cur21 * (savedl12 + savedl13)) / (cur11 * (savedl22 + savedl23));
1722 	  base_probl *= (cur21 * savedl12) / (cur11 * savedl22);
1723 	  savedl12 -= 1;
1724 	} else {
1725 	  savedl23 += 1;
1726 	  row_prob *= (cur21 * (savedl12 + savedl13)) / (cur11 * (savedl22 + savedl23));
1727 	  base_probl *= (cur21 * savedl13) / (cur11 * savedl23);
1728 	  savedl13 -= 1;
1729 	}
1730 	cur21 -= 1;
1731       }
1732       if (fisher23_tailsum(&base_probl, &savedl12, &savedl13, &savedl22, &savedl23, &dxx, &tie_ct, 0)) {
1733 	break;
1734       }
1735       tprob += dxx;
1736       if (dir) {
1737 	if (savedr22) {
1738 	  savedr12 += 1;
1739 	  base_probr *= ((cur11 + 1) * savedr22) / (cur21 * savedr12);
1740 	  savedr22 -= 1;
1741 	} else {
1742 	  savedr13 += 1;
1743 	  base_probr *= ((cur11 + 1) * savedr23) / (cur21 * savedr13);
1744 	  savedr23 -= 1;
1745 	}
1746       } else {
1747 	if (savedr13) {
1748 	  savedr23 += 1;
1749 	  base_probr *= ((cur21 + 1) * savedr13) / (cur11 * savedr23);
1750 	  savedr13 -= 1;
1751 	} else {
1752 	  savedr22 += 1;
1753 	  base_probr *= ((cur21 + 1) * savedr12) / (cur11 * savedr22);
1754 	  savedr12 -= 1;
1755 	}
1756       }
1757       fisher23_tailsum(&base_probr, &savedr12, &savedr13, &savedr22, &savedr23, &dyy, &tie_ct, 1);
1758       tprob += dyy;
1759       cprob += row_prob - dxx - dyy;
1760       if (cprob == INFINITY) {
1761 	return 0;
1762       }
1763     }
1764     if (!ukk) {
1765       continue;
1766     }
1767     savedl12 += savedl13;
1768     savedl22 += savedl23;
1769     if (dir) {
1770       while (1) {
1771 	preaddp = tprob;
1772 	tprob += row_prob;
1773 	if (tprob <= preaddp) {
1774 	  break;
1775 	}
1776 	cur21 += 1;
1777 	savedl12 += 1;
1778 	row_prob *= (cur11 * savedl22) / (cur21 * savedl12);
1779 	cur11 -= 1;
1780 	savedl22 -= 1;
1781       }
1782     } else {
1783       while (1) {
1784 	preaddp = tprob;
1785 	tprob += row_prob;
1786 	if (tprob <= preaddp) {
1787 	  break;
1788 	}
1789 	cur11 += 1;
1790 	savedl22 += 1;
1791 	row_prob *= (cur21 * savedl12) / (cur11 * savedl22);
1792 	cur21 -= 1;
1793 	savedl12 -= 1;
1794       }
1795     }
1796   }
1797   if (!midp) {
1798     return tprob / (tprob + cprob);
1799   } else {
1800     return (tprob - ((1 - FISHER_EPSILON) * EXACT_TEST_BIAS * 0.5) * ((int32_t)tie_ct)) / (tprob + cprob);
1801   }
1802 }
1803 
chi22_get_coeffs(intptr_t row1_sum,intptr_t col1_sum,intptr_t total,double * expm11p,double * recip_sump)1804 void chi22_get_coeffs(intptr_t row1_sum, intptr_t col1_sum, intptr_t total, double* expm11p, double* recip_sump) {
1805   // chisq = (m11 - expm11)^2 * recip_sum
1806   // (see discussion for chi22_precomp_val_bounds() below.)
1807   //
1808   // expm11 = row1_sum * col1_sum / total
1809   // expm12 = row1_sum * col2_sum / total, etc.
1810   // recip_sum = 1 / expm11 + 1 / expm12 + 1 / expm21 + 1 / expm22
1811   // = total * (1 / (row1_sum * col1_sum) + 1 / (row1_sum * col2_sum) +
1812   //            1 / (row2_sum * col1_sum) + 1 / (row2_sum * col2_sum))
1813   // = total^3 / (row1_sum * col1_sum * row2_sum * col2_sum)
1814   double m11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1815   double denom = m11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1816   double dxx;
1817   if (denom != 0) {
1818     dxx = total;
1819     *expm11p = m11_numer / dxx;
1820     *recip_sump = dxx * dxx * dxx / denom;
1821   } else {
1822     // since an entire row or column is zero, either m11 or m22 is zero
1823     // row1_sum + col1_sum - total = m11 - m22
1824     if (row1_sum + col1_sum < total) {
1825       *expm11p = 0;
1826     } else {
1827       *expm11p = row1_sum + col1_sum - total;
1828     }
1829     *recip_sump = 0;
1830   }
1831 }
1832 
chi22_eval(intptr_t m11,intptr_t row1_sum,intptr_t col1_sum,intptr_t total)1833 double chi22_eval(intptr_t m11, intptr_t row1_sum, intptr_t col1_sum, intptr_t total) {
1834   double expm11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1835   double denom = expm11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1836   double dxx;
1837   double dyy;
1838   if (denom != 0) {
1839     dxx = total;
1840     dyy = m11 * dxx - expm11_numer; // total * (m11 - expm11)
1841     return (dyy * dyy * dxx) / denom;
1842   } else {
1843     return 0;
1844   }
1845 }
1846 
chi22_evalx(intptr_t m11,intptr_t row1_sum,intptr_t col1_sum,intptr_t total)1847 double chi22_evalx(intptr_t m11, intptr_t row1_sum, intptr_t col1_sum, intptr_t total) {
1848   // PLINK emulation.  returns -9 instead of 0 if row1_sum, row2_sum, col1_sum,
1849   // or col2_sum is zero, for identical "NA" reporting.
1850   double expm11_numer = ((uint64_t)row1_sum) * ((uint64_t)col1_sum);
1851   double denom = expm11_numer * (((uint64_t)(total - row1_sum)) * ((uint64_t)(total - col1_sum)));
1852   double dxx;
1853   double dyy;
1854   if (denom != 0) {
1855     dxx = total;
1856     dyy = m11 * dxx - expm11_numer; // total * (m11 - expm11)
1857     return (dyy * dyy * dxx) / denom;
1858   } else {
1859     return -9;
1860   }
1861 }
1862 
chi22_precomp_val_bounds(double chisq,intptr_t row1_sum,intptr_t col1_sum,intptr_t total,uint32_t * bounds,double * coeffs)1863 void chi22_precomp_val_bounds(double chisq, intptr_t row1_sum, intptr_t col1_sum, intptr_t total, uint32_t* bounds, double* coeffs) {
1864   // Treating m11 as the only variable, this returns the minimum and (maximum +
1865   // 1) values of m11 which produce smaller chisq statistics than given in
1866   // bounds[0] and bounds[1] respectively, and smaller-or-equal interval
1867   // bounds in bounds[2] and bounds[3].
1868   double expm11;
1869   double recip_sum;
1870   double cur11;
1871   double dxx;
1872   intptr_t ceil11;
1873   intptr_t lii;
1874   chi22_get_coeffs(row1_sum, col1_sum, total, &expm11, &recip_sum);
1875   if (coeffs) {
1876     coeffs[0] = expm11;
1877     coeffs[1] = recip_sum;
1878   }
1879   if (recip_sum == 0) {
1880     // sum-0 row or column, no freedom at all
1881     bounds[0] = (intptr_t)expm11;
1882     bounds[1] = bounds[0];
1883     bounds[2] = bounds[0];
1884     if (chisq == 0) {
1885       bounds[3] = bounds[0] + 1;
1886     } else {
1887       bounds[3] = bounds[0];
1888     }
1889     return;
1890   }
1891 
1892   // double cur_stat = (cur11 - exp11) * (cur11 - exp11) * recipx11 + (cur12 - exp12) * (cur12 - exp12) * recipx12 + (cur21 - exp21) * (cur21 - exp21) * recipx21 + (cur22 - exp22) * (cur22 - exp22) * recipx22;
1893   // However, we have
1894   //   cur11 - exp11 = -(cur12 - exp12) = -(cur21 - exp21) = cur22 - exp22
1895   // So the chisq statistic reduces to
1896   //   (cur11 - exp11)^2 * (recipx11 + recipx12 + recipx21 + recipx22).
1897 
1898   ceil11 = MINV(row1_sum, col1_sum);
1899   // chisq = (cur11 - expm11)^2 * recip_sum
1900   // -> expm11 +/- sqrt(chisq / recip_sum) = cur11
1901   recip_sum = sqrt(chisq / recip_sum);
1902   cur11 = expm11 - recip_sum;
1903   dxx = cur11 + 1 - BIG_EPSILON;
1904   if (dxx < 0) {
1905     bounds[0] = 0;
1906     bounds[2] = 0;
1907   } else {
1908     lii = (intptr_t)dxx;
1909     bounds[2] = lii;
1910     if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
1911       bounds[0] = lii + 1;
1912     } else {
1913       bounds[0] = lii;
1914     }
1915   }
1916   cur11 = expm11 + recip_sum;
1917   if (cur11 > ceil11 + BIG_EPSILON) {
1918     bounds[1] = ceil11 + 1;
1919     bounds[3] = bounds[1];
1920   } else {
1921     dxx = cur11 + 1 - BIG_EPSILON;
1922     lii = (intptr_t)dxx;
1923     bounds[1] = lii;
1924     if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
1925       bounds[3] = lii + 1;
1926     } else {
1927       bounds[3] = lii;
1928     }
1929   }
1930 }
1931 
chi23_eval(intptr_t m11,intptr_t m12,intptr_t row1_sum,intptr_t col1_sum,intptr_t col2_sum,intptr_t total)1932 double chi23_eval(intptr_t m11, intptr_t m12, intptr_t row1_sum, intptr_t col1_sum, intptr_t col2_sum, intptr_t total) {
1933   // assumes no sum-zero row
1934   intptr_t m13 = row1_sum - m11 - m12;
1935   intptr_t col3_sum = total - col1_sum - col2_sum;
1936   double col1_sumd;
1937   double col2_sumd;
1938   double col3_sumd;
1939   double tot_recip;
1940   double dxx;
1941   double expect;
1942   double delta;
1943   double chisq;
1944   col1_sumd = col1_sum;
1945   col2_sumd = col2_sum;
1946   col3_sumd = col3_sum;
1947   tot_recip = 1.0 / ((double)total);
1948   dxx = row1_sum * tot_recip;
1949   expect = dxx * col1_sumd;
1950   delta = m11 - expect;
1951   chisq = delta * delta / expect;
1952   expect = dxx * col2_sumd;
1953   delta = m12 - expect;
1954   chisq += delta * delta / expect;
1955   expect = dxx * col3_sumd;
1956   delta = m13 - expect;
1957   chisq += delta * delta / expect;
1958   dxx = (total - row1_sum) * tot_recip;
1959   expect = dxx * col1_sumd;
1960   delta = (col1_sum - m11) - expect;
1961   chisq += delta * delta / expect;
1962   expect = dxx * col2_sumd;
1963   delta = (col2_sum - m12) - expect;
1964   chisq += delta * delta / expect;
1965   expect = dxx * col3_sumd;
1966   delta = (col3_sum - m13) - expect;
1967   chisq += delta * delta / expect;
1968   if (chisq < (SMALL_EPSILON * SMALL_EPSILON)) {
1969     return 0;
1970   }
1971   return chisq;
1972 }
1973 
chi23_evalx(intptr_t m11,intptr_t m12,intptr_t m13,intptr_t m21,intptr_t m22,intptr_t m23,double * chip,uint32_t * dfp)1974 void chi23_evalx(intptr_t m11, intptr_t m12, intptr_t m13, intptr_t m21, intptr_t m22, intptr_t m23, double* chip, uint32_t* dfp) {
1975   // Slightly different from PLINK calculation, since it detects lone nonzero
1976   // columns.
1977   intptr_t row1_sum = m11 + m12 + m13;
1978   intptr_t row2_sum = m21 + m22 + m23;
1979   intptr_t col1_sum = m11 + m21;
1980   intptr_t col2_sum = m12 + m22;
1981   intptr_t col3_sum = m13 + m23;
1982   intptr_t total;
1983   double col1_sumd;
1984   double col2_sumd;
1985   double col3_sumd;
1986   double tot_recip;
1987   double dxx;
1988   double expect;
1989   double delta;
1990   double chisq;
1991   if ((!row1_sum) || (!row2_sum)) {
1992     *chip = -9;
1993     *dfp = 0;
1994     return;
1995   }
1996   total = row1_sum + row2_sum;
1997   if (!col1_sum) {
1998     *chip = chi22_evalx(m12, row1_sum, col2_sum, total);
1999     if (*chip != -9) {
2000       *dfp = 1;
2001     } else {
2002       *dfp = 0;
2003     }
2004     return;
2005   } else if ((!col2_sum) || (!col3_sum)) {
2006     *chip = chi22_evalx(m11, row1_sum, col1_sum, total);
2007     if (*chip != -9) {
2008       *dfp = 1;
2009     } else {
2010       *dfp = 0;
2011     }
2012     return;
2013   }
2014   col1_sumd = col1_sum;
2015   col2_sumd = col2_sum;
2016   col3_sumd = col3_sum;
2017   tot_recip = 1.0 / ((double)total);
2018   dxx = row1_sum * tot_recip;
2019   expect = dxx * col1_sumd;
2020   delta = m11 - expect;
2021   chisq = delta * delta / expect;
2022   expect = dxx * col2_sumd;
2023   delta = m12 - expect;
2024   chisq += delta * delta / expect;
2025   expect = dxx * col3_sumd;
2026   delta = m13 - expect;
2027   chisq += delta * delta / expect;
2028   dxx = row2_sum * tot_recip;
2029   expect = dxx * col1_sumd;
2030   delta = m21 - expect;
2031   chisq += delta * delta / expect;
2032   expect = dxx * col2_sumd;
2033   delta = m22 - expect;
2034   chisq += delta * delta / expect;
2035   expect = dxx * col3_sumd;
2036   delta = m23 - expect;
2037   chisq += delta * delta / expect;
2038   if (chisq < (SMALL_EPSILON * SMALL_EPSILON)) {
2039     chisq = 0;
2040   }
2041   *chip = chisq;
2042   *dfp = 2;
2043 }
2044 
ca_trend_eval(intptr_t case_dom_ct,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total)2045 double ca_trend_eval(intptr_t case_dom_ct, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total) {
2046   // case_dom_ct is an allele count (2 * homa2 + het), while other inputs are
2047   // observation counts.
2048   //
2049   // If case_missing_ct is fixed,
2050   //   row1_sum = case ct
2051   //   col1_sum = A2 ct
2052   //   case_ct * ctrl_ct * REC_ct * DOM_ct
2053   //   REC_ct = 2 * obs_11 + obs_12
2054   //   DOM_ct = 2 * obs_22 + obs_12
2055   //   CA = (obs_U / obs_T) * (case REC ct) - (obs_A / obs_T) * (ctrl DOM ct)
2056   //      = (case A2) * (obs_U / obs_T) - (obs_A / obs_T) * (DOM ct - case DOM)
2057   //      = (case A2) * (obs_U / obs_T) + (case DOM) * (obs_A / obs_T) - DOM*(A/T)
2058   //      = (case A2 ct) - total A2 ct * (A/T)
2059   //   CAT = CA * obs_T
2060   //   varCA_recip = obs_T * obs_T * obs_T /
2061   //     (obs_A * obs_U * (obs_T * (obs_12 + 4 * obs_22) - DOMct * DOMct))
2062   //   trend statistic = CAT * CAT * [varCA_recip / obs_T^2]
2063   double dom_ct = het_ct + 2 * homdom_ct;
2064   double totald = total;
2065   double case_ctd = case_ct;
2066   double cat = case_dom_ct * totald - dom_ct * case_ctd;
2067   double dxx = totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ct * dom_ct;
2068 
2069   // This should never be called with dxx == 0 (which happens when two columns
2070   // are all-zero).  Use ca_trend_evalx() to check for that.
2071   dxx *= case_ctd * (totald - case_ctd);
2072   return cat * cat * totald / dxx;
2073 }
2074 
ca_trend_evalx(intptr_t case_dom_ct,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total)2075 double ca_trend_evalx(intptr_t case_dom_ct, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total) {
2076   double dom_ct = het_ct + 2 * homdom_ct;
2077   double totald = total;
2078   double case_ctd = case_ct;
2079   double cat = case_dom_ct * totald - dom_ct * case_ctd;
2080   double dxx = totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ct * dom_ct;
2081   if (dxx != 0) {
2082     dxx *= case_ctd * (totald - case_ctd);
2083     return cat * cat * totald / dxx;
2084   } else {
2085     return -9;
2086   }
2087 }
2088 
ca_trend_precomp_val_bounds(double chisq,intptr_t case_ct,intptr_t het_ct,intptr_t homdom_ct,intptr_t total,uint32_t * bounds,double * coeffs)2089 void ca_trend_precomp_val_bounds(double chisq, intptr_t case_ct, intptr_t het_ct, intptr_t homdom_ct, intptr_t total, uint32_t* bounds, double* coeffs) {
2090   // If case_missing_ct is fixed,
2091   //   row1_sum = case ct
2092   //   col1_sum = DOM ct
2093   //   case_ct * ctrl_ct * REC_ct * DOM_ct
2094   //   REC_ct = 2 * obs_11 + obs_12
2095   //   DOM_ct = 2 * obs_22 + obs_12
2096   //   CA = (obs_U / obs_T) * (case DOM ct) - (obs_A / obs_T) * (ctrl DOM ct)
2097   //      = (case DOM) * (obs_U / obs_T) - (obs_A / obs_T) * (DOM ct - case DOM)
2098   //      = (case DOM) * (obs_U / obs_T) + (case DOM) * (obs_A / obs_T) - DOM*(A/T)
2099   //      = (case DOM ct) - total DOM ct * (A/T)
2100   //   varCA_recip = obs_T * obs_T * obs_T /
2101   //     (obs_A * obs_U * (obs_T * (obs_12 + 4 * obs_22) - DOMct * DOMct))
2102   //   trend statistic = CA * CA * varCA_recip
2103   intptr_t dom_ct = het_ct + 2 * homdom_ct;
2104   double dom_ctd = dom_ct;
2105   double totald = total;
2106   double case_ctd = case_ct;
2107   double tot_recip = 1.0 / totald;
2108   double expm11 = dom_ctd * case_ctd * tot_recip;
2109   double dxx = case_ctd * (totald - case_ctd) * (totald * (het_ct + 4 * ((int64_t)homdom_ct)) - dom_ctd * dom_ctd);
2110   double varca_recip;
2111   double cur11;
2112   intptr_t ceil11;
2113   intptr_t lii;
2114   if (dxx == 0) {
2115     // bounds/coefficients should never be referenced in this case
2116     return;
2117   }
2118   varca_recip = totald * totald * totald / dxx;
2119   if (coeffs) {
2120     coeffs[0] = expm11;
2121     coeffs[1] = varca_recip;
2122   }
2123 
2124   // statistic: (cur11 - expm11)^2 * varca_recip
2125   ceil11 = case_ct * 2;
2126   if (dom_ct < ceil11) {
2127     ceil11 = dom_ct;
2128   }
2129   // chisq = (cur11 - expm11)^2 * varca_recip
2130   // -> expm11 +/- sqrt(chisq / varca_recip) = cur11
2131   varca_recip = sqrt(chisq / varca_recip);
2132   cur11 = expm11 - varca_recip;
2133   dxx = cur11 + 1 - BIG_EPSILON;
2134   if (dxx < 0) {
2135     bounds[0] = 0;
2136     bounds[2] = 0;
2137   } else {
2138     lii = (intptr_t)dxx;
2139     bounds[2] = lii;
2140     if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
2141       bounds[0] = lii + 1;
2142     } else {
2143       bounds[0] = lii;
2144     }
2145   }
2146   cur11 = expm11 + varca_recip;
2147   if (cur11 > ceil11 + BIG_EPSILON) {
2148     bounds[1] = ceil11 + 1;
2149     bounds[3] = bounds[1];
2150   } else {
2151     dxx = cur11 + 1 - BIG_EPSILON;
2152     lii = (intptr_t)dxx;
2153     bounds[1] = lii;
2154     if (lii == (intptr_t)(cur11 + BIG_EPSILON)) {
2155       bounds[3] = lii + 1;
2156     } else {
2157       bounds[3] = lii;
2158     }
2159   }
2160 }
2161 
linear_hypothesis_chisq(uintptr_t constraint_ct,uintptr_t param_ct,double * constraints_con_major,double * coef,double * cov_matrix,double * param_df_buf,double * param_df_buf2,double * df_df_buf,MATRIX_INVERT_BUF1_TYPE * mi_buf,double * df_buf,double * chisq_ptr)2162 uint32_t linear_hypothesis_chisq(uintptr_t constraint_ct, uintptr_t param_ct, double* constraints_con_major, double* coef, double* cov_matrix, double* param_df_buf, double* param_df_buf2, double* df_df_buf, MATRIX_INVERT_BUF1_TYPE* mi_buf, double* df_buf, double* chisq_ptr) {
2163   // See PLINK model.cpp Model::linearHypothesis().
2164   //
2165   // outer = df_buf
2166   // inner = df_df_buf
2167   // tmp = param_df_buf
2168   // mi_buf only needs to be of length constraint_ct
2169   //
2170   // Since no PLINK function ever calls this with nonzero h[] values, this just
2171   // takes a df parameter for now; it's trivial to switch to the more general
2172   // interface later.
2173   double* dptr = constraints_con_major;
2174   double* dptr2;
2175   uintptr_t constraint_idx;
2176   uintptr_t constraint_idx2;
2177   uintptr_t param_idx;
2178   double dxx;
2179   double dyy;
2180   for (constraint_idx = 0; constraint_idx < constraint_ct; constraint_idx++) {
2181     dxx = 0;
2182     dptr2 = coef;
2183     for (param_idx = 0; param_idx < param_ct; param_idx++) {
2184       dxx += (*dptr++) * (*dptr2++);
2185     }
2186     df_buf[constraint_idx] = dxx;
2187   }
2188   // temporarily set param_df_buf2[][] to H-transpose
2189   transpose_copy(constraint_ct, param_ct, constraints_con_major, param_df_buf2);
2190   col_major_matrix_multiply(constraint_ct, param_ct, param_ct, param_df_buf2, cov_matrix, param_df_buf);
2191   // tmp[][] is now param-major
2192   col_major_matrix_multiply(constraint_ct, constraint_ct, param_ct, param_df_buf, constraints_con_major, df_df_buf);
2193 
2194   if (invert_matrix((uint32_t)constraint_ct, df_df_buf, mi_buf, param_df_buf2)) {
2195     return 1;
2196   }
2197   dxx = 0; // result
2198   dptr = df_df_buf;
2199   for (constraint_idx = 0; constraint_idx < constraint_ct; constraint_idx++) {
2200     dyy = 0; // tmp2[c]
2201     dptr2 = df_buf;
2202     for (constraint_idx2 = 0; constraint_idx2 < constraint_ct; constraint_idx2++) {
2203       dyy += (*dptr++) * (*dptr2++);
2204     }
2205     dxx += dyy * df_buf[constraint_idx];
2206   }
2207   *chisq_ptr = dxx;
2208   return 0;
2209 }
2210 
binom_2sided(uint32_t succ,uint32_t obs,uint32_t midp)2211 double binom_2sided(uint32_t succ, uint32_t obs, uint32_t midp) {
2212   // straightforward to generalize this to any success probability
2213   double cur_succ_t2 = (int32_t)succ;
2214   double cur_fail_t2 = (int32_t)(obs - succ);
2215   double tailp = (1 - SMALL_EPSILON) * EXACT_TEST_BIAS;
2216   double centerp = 0;
2217   double lastp2 = tailp;
2218   double lastp1 = tailp;
2219   int32_t tie_ct = 1;
2220   // double rate_mult_incr = rate / (1 - rate);
2221   // double rate_mult_decr = (1 - rate) / rate;
2222   double cur_succ_t1;
2223   double cur_fail_t1;
2224   double preaddp;
2225   if (!obs) {
2226     return midp? 0.5 : 1;
2227   }
2228   if (obs < succ * 2) {
2229   // if (obs * rate < succ) {
2230     while (cur_succ_t2 > 0.5) {
2231       cur_fail_t2 += 1;
2232       lastp2 *= cur_succ_t2 / cur_fail_t2;
2233       // lastp2 *= rate_mult_decr * cur_succ_t2 / cur_fail_t2;
2234       cur_succ_t2 -= 1;
2235       if (lastp2 < EXACT_TEST_BIAS) {
2236 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
2237 	  tie_ct++;
2238 	}
2239 	tailp += lastp2;
2240 	break;
2241       }
2242       centerp += lastp2;
2243       if (centerp == INFINITY) {
2244 	return 0;
2245       }
2246     }
2247     if ((centerp == 0) && (!midp)) {
2248       return 1;
2249     }
2250     while (cur_succ_t2 > 0.5) {
2251       cur_fail_t2 += 1;
2252       lastp2 *= cur_succ_t2 / cur_fail_t2;
2253       // lastp2 *= rate_mult_decr * cur_succ_t2 / cur_fail_t2;
2254       cur_succ_t2 -= 1;
2255       preaddp = tailp;
2256       tailp += lastp2;
2257       if (tailp <= preaddp) {
2258 	break;
2259       }
2260     }
2261     cur_succ_t1 = (int32_t)(succ + 1);
2262     cur_fail_t1 = (int32_t)(obs - succ);
2263     while (cur_fail_t1 > 0.5) {
2264       lastp1 *= cur_fail_t1 / cur_succ_t1;
2265       // lastp1 *= rate_mult_incr * cur_fail_t1 / cur_succ_t1;
2266       preaddp = tailp;
2267       tailp += lastp1;
2268       if (tailp <= preaddp) {
2269 	break;
2270       }
2271       cur_succ_t1 += 1;
2272       cur_fail_t1 -= 1;
2273     }
2274   } else {
2275     while (cur_fail_t2 > 0.5) {
2276       cur_succ_t2++;
2277       lastp2 *= cur_fail_t2 / cur_succ_t2;
2278       // lastp2 *= rate_mult_incr * cur_fail_t2 / cur_succ_t2;
2279       cur_fail_t2--;
2280       if (lastp2 < EXACT_TEST_BIAS) {
2281 	if (lastp2 > (1 - 2 * SMALL_EPSILON) * EXACT_TEST_BIAS) {
2282 	  tie_ct++;
2283 	}
2284 	tailp += lastp2;
2285 	break;
2286       }
2287       centerp += lastp2;
2288       if (centerp == INFINITY) {
2289 	return 0;
2290       }
2291     }
2292     if ((centerp == 0) && (!midp)) {
2293       return 1;
2294     }
2295     while (cur_fail_t2 > 0.5) {
2296       cur_succ_t2 += 1;
2297       lastp2 *= cur_fail_t2 / cur_succ_t2;
2298       // lastp2 *= rate_mult_incr * cur_fail_t2 / cur_succ_t2;
2299       cur_fail_t2 -= 1;
2300       preaddp = tailp;
2301       tailp += lastp2;
2302       if (tailp <= preaddp) {
2303 	break;
2304       }
2305     }
2306     cur_succ_t1 = (int32_t)succ;
2307     cur_fail_t1 = (int32_t)(obs - succ);
2308     while (cur_succ_t1 > 0.5) {
2309       cur_fail_t1 += 1;
2310       lastp1 *= cur_succ_t1 / cur_fail_t1;
2311       // lastp1 *= rate_mult_decr * cur_succ_t1 / cur_fail_t1;
2312       preaddp = tailp;
2313       tailp += lastp1;
2314       if (tailp <= preaddp) {
2315 	break;
2316       }
2317       cur_succ_t1 -= 1;
2318     }
2319   }
2320   if (!midp) {
2321     return tailp / (tailp + centerp);
2322   } else {
2323     return (tailp - ((1 - SMALL_EPSILON) * EXACT_TEST_BIAS * 0.5) * tie_ct) / (tailp + centerp);
2324   }
2325 }
2326