diff --git a/validation/gwas/method/hwe/Makefile b/validation/gwas/method/hwe/Makefile new file mode 100644 index 000000000..eedc04d4d --- /dev/null +++ b/validation/gwas/method/hwe/Makefile @@ -0,0 +1,12 @@ +all: clean chwe.so + +clean: + rm -f *.o *.so + +chwe.so: chwe.o + gcc -shared -o libchwe.so chwe.o + +chwe.o: chwe.c + gcc -c -Wall -Werror -fpic chwe.c + + diff --git a/validation/gwas/method/hwe/README.md b/validation/gwas/method/hwe/README.md new file mode 100644 index 000000000..491eb3c45 --- /dev/null +++ b/validation/gwas/method/hwe/README.md @@ -0,0 +1,24 @@ +## HWE Exact Test Validation + +This validation produces simulated genotype counts and corresponding HWE statistics from the (C) implementation described in [Wigginton et al. 2005](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1199378). + +The `invoke` [tasks](tasks.py) will compile the C code, simulate genotype counts (inputs for unit tests), and attach p values (outputs for unit tests) from the C code to the genotype counts, as a dataframe. + +The [hwe_unit_test.ipynb](hwe_unit_test.ipynb) is only instructive and shows how to debug and possibly extend test cases, perhaps validating inputs/outputs on a scale that wouldn't be included in unit testing. + +To export the unit test data, all steps can be run as follows: + +```bash +> invoke compile simulate export +Building reference C library +rm -f *.o *.so +gcc -c -Wall -Werror -fpic chwe.c +gcc -shared -o libchwe.so chwe.o +Build complete +Generating unit test data +Unit test data written to data/sim_01.csv +Exporting test data to /home/jovyan/work/repos/sgkit/sgkit/tests/test_hwe +Clearing test datadir at /home/jovyan/work/repos/sgkit/sgkit/tests/test_hwe +Copying data/sim_01.csv to /home/jovyan/work/repos/sgkit/sgkit/tests/test_hwe/sim_01.csv +Export complete +``` \ No newline at end of file diff --git a/validation/gwas/method/hwe/__init__.py b/validation/gwas/method/hwe/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/validation/gwas/method/hwe/chwe.c b/validation/gwas/method/hwe/chwe.c new file mode 100644 index 000000000..50c327e7d --- /dev/null +++ b/validation/gwas/method/hwe/chwe.c @@ -0,0 +1,98 @@ +// Lift from http://csg.sph.umich.edu/abecasis/Exact/snp_hwe.c +#include +#include + +double hwep(int obs_hets, int obs_hom1, int obs_hom2){ + if (obs_hom1 < 0 || obs_hom2 < 0 || obs_hets < 0) + { + printf("FATAL ERROR - SNP-HWE: Current genotype configuration (%d %d %d ) includes a" + " negative count", obs_hets, obs_hom1, obs_hom2); + exit(EXIT_FAILURE); + } + + int obs_homc = obs_hom1 < obs_hom2 ? obs_hom2 : obs_hom1; + int obs_homr = obs_hom1 < obs_hom2 ? obs_hom1 : obs_hom2; + + int rare_copies = 2 * obs_homr + obs_hets; + int genotypes = obs_hets + obs_homc + obs_homr; + + double * het_probs = (double *) malloc((size_t) (rare_copies + 1) * sizeof(double)); + if (het_probs == NULL) + { + printf("FATAL ERROR - SNP-HWE: Unable to allocate array for heterozygote probabilities" ); + exit(EXIT_FAILURE); + } + + int i; + for (i = 0; i <= rare_copies; i++) + het_probs[i] = 0.0; + + /* start at midpoint */ + int mid = rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes); + + /* check to ensure that midpoint and rare alleles have same parity */ + if ((rare_copies & 1) ^ (mid & 1)) + mid++; + + int curr_hets = mid; + int curr_homr = (rare_copies - mid) / 2; + int curr_homc = genotypes - curr_hets - curr_homr; + + het_probs[mid] = 1.0; + double sum = het_probs[mid]; + for (curr_hets = mid; curr_hets > 1; curr_hets -= 2) + { + het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0) + / (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0)); + sum += het_probs[curr_hets - 2]; + + /* 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote */ + curr_homr++; + curr_homc++; + } + + curr_hets = mid; + curr_homr = (rare_copies - mid) / 2; + curr_homc = genotypes - curr_hets - curr_homr; + for (curr_hets = mid; curr_hets <= rare_copies - 2; curr_hets += 2) + { + het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc + /((curr_hets + 2.0) * (curr_hets + 1.0)); + sum += het_probs[curr_hets + 2]; + + /* add 2 heterozygotes for next iteration -> subtract one rare, one common homozygote */ + curr_homr--; + curr_homc--; + } + + for (i = 0; i <= rare_copies; i++) + het_probs[i] /= (sum > 0 ? sum : 1e-128); + + /* alternate p-value calculation for p_hi/p_lo + double p_hi = het_probs[obs_hets]; + for (i = obs_hets + 1; i <= rare_copies; i++) + p_hi += het_probs[i]; + + double p_lo = het_probs[obs_hets]; + for (i = obs_hets - 1; i >= 0; i--) + p_lo += het_probs[i]; + + + double p_hi_lo = p_hi < p_lo ? 2.0 * p_hi : 2.0 * p_lo; + */ + + double p_hwe = 0.0; + /* p-value calculation for p_hwe */ + for (i = 0; i <= rare_copies; i++) + { + if (het_probs[i] > het_probs[obs_hets]) + continue; + p_hwe += het_probs[i]; + } + + p_hwe = p_hwe > 1.0 ? 1.0 : p_hwe; + + free(het_probs); + + return p_hwe; +} diff --git a/validation/gwas/method/hwe/chwe.o b/validation/gwas/method/hwe/chwe.o new file mode 100644 index 000000000..a88353de1 Binary files /dev/null and b/validation/gwas/method/hwe/chwe.o differ diff --git a/validation/gwas/method/hwe/data/sim_01.csv b/validation/gwas/method/hwe/data/sim_01.csv new file mode 100644 index 000000000..97bbb05d7 --- /dev/null +++ b/validation/gwas/method/hwe/data/sim_01.csv @@ -0,0 +1,201 @@ +n_het,n_hom_1,n_hom_2,p +1,0,0,1.0 +51,27,26,0.845926828898329 +101,47,56,0.88930424473698 +151,71,99,0.3684874920023832 +201,137,91,0.28344213097688165 +251,154,128,0.1937805064539818 +301,158,201,0.03464156123462855 +351,115,117,9.778051416840305e-07 +401,123,253,0.09001690989721362 +451,275,292,0.00027739277561481717 +501,346,310,5.785950777451914e-06 +551,267,337,0.15629004579821143 +601,208,334,0.031198370505333042 +651,232,441,0.7782347792381512 +701,356,326,0.6282550560763566 +751,304,457,0.9172428161320041 +801,386,422,0.8810870885296063 +851,261,465,0.00012095475067049487 +901,490,492,0.06519606837210713 +951,644,544,4.879275092149919e-07 +1001,444,475,0.0614320395959403 +1051,608,340,0.0016721676932256468 +1101,623,625,0.002583239824235713 +1151,442,404,9.22666204121466e-12 +1201,511,535,0.0011559180266171758 +1251,660,594,0.999999999999999 +1301,904,443,0.5213032029206021 +1351,518,492,2.4059583160333022e-12 +1401,786,562,0.17901712606288267 +1451,705,577,0.0008509642369958977 +1501,545,516,4.0878258413432674e-18 +1551,872,551,0.0026562564901761523 +1601,606,716,1.7304716281165935e-07 +1651,1037,559,0.026125733725720014 +1701,1080,575,0.03402821665047714 +1751,1209,853,5.572653418323456e-06 +1801,1243,976,2.2818771131803835e-10 +1851,1102,584,3.710099688663844e-05 +1901,785,661,2.0279602738807187e-15 +1951,816,677,2.8569078183446102e-15 +2001,854,931,0.00044371590402349367 +2051,667,1183,1.2488023350848009e-05 +2101,1106,853,0.01396809036282682 +2151,1095,726,6.691337471001894e-09 +2201,1167,1478,9.432197243064005e-10 +2251,962,1276,0.610161450028642 +2301,811,1349,0.0019349995547350155 +2351,977,877,1.444515271924519e-14 +2401,1283,739,1.1565179920072941e-11 +2451,1547,739,6.164163620718031e-06 +2501,1428,1020,0.21851919021140662 +2551,1515,1747,2.6790289570097594e-20 +2601,1039,1379,0.003812650444672208 +2651,1423,1402,0.018717042880622592 +2701,1051,1839,0.2872761239911692 +2751,1317,1756,0.00014693980061811332 +2801,1623,1173,0.5905238038361554 +2851,1783,1307,0.009007126080770289 +2901,1892,1544,7.821435236313817e-11 +2951,1926,1702,1.5502927534087137e-16 +3001,1770,1502,0.0010079660986331209 +3051,2082,1701,6.517752718327402e-18 +3101,1456,1682,0.7226272913884884 +3151,969,1325,3.506219413030305e-33 +3201,1805,1331,0.21569403737115994 +3251,1778,1532,0.5365544398754494 +3301,1169,1384,4.982094499906052e-23 +3351,1769,1797,0.010072541488350021 +3401,1801,1908,0.00028325611912255233 +3451,1935,1630,0.23164761482387491 +3501,2305,1565,0.0005408313105776889 +3551,1684,2332,2.2869325074804336e-06 +3601,2241,2094,2.3394325551477485e-16 +3651,1241,2438,0.04636664490239822 +3701,2167,2588,1.9557166448260108e-29 +3751,1349,2427,0.13201350023523056 +3801,1387,2076,1.9301027903365754e-06 +3851,1346,2461,0.017224681698312864 +3901,2430,2058,4.440082083118905e-10 +3951,1828,1294,2.647713527040051e-25 +4001,2316,1926,0.014256792586706637 +4051,2385,2619,2.3465822237767657e-23 +4101,2830,2634,5.149856509652088e-44 +4151,1264,1843,4.48467880717666e-38 +4201,2486,1548,0.002389012903737995 +4251,2161,1367,3.000197368786019e-20 +4301,1634,1322,9.919986792554593e-58 +4351,2686,1695,0.38577132434045974 +4401,1928,2954,0.0001199016650457927 +4451,2589,1391,1.0958652901993421e-12 +4501,1646,2469,4.5269662611402833e-07 +4551,2416,1798,4.521174352826995e-05 +4601,3099,2510,3.742694664287275e-22 +4651,2391,2492,0.01847696637280268 +4701,2783,1996,0.9012359250831754 +4751,2182,1824,4.711674213296269e-16 +4801,1797,3253,0.7262132846779064 +4851,2890,2406,2.680676415905425e-05 +4901,1916,1968,1.9576351291926123e-27 +4951,1600,2345,3.4468103297849734e-30 +5001,2124,2893,0.6729732567000968 +5051,2278,1878,1.9099808504594497e-21 +5101,1580,1667,2.7112967783602924e-92 +5151,2945,2480,0.0133474541803742 +5201,2676,3425,9.423282617882861e-16 +5251,3655,2030,0.06045595125773339 +5301,2996,2148,0.027356813588423446 +5351,1649,3228,3.454974198424328e-13 +5401,2311,2448,1.8596430888102382e-10 +5451,2918,3447,2.2579116525279354e-16 +5501,3034,3570,8.880020019491782e-23 +5551,2272,3437,0.7316458577549785 +5601,2096,3814,0.6201032727285236 +5651,3249,2182,0.002020139572176394 +5701,3870,3376,3.265438252095903e-41 +5751,2309,2216,9.550965561735892e-34 +5801,2942,1799,1.584359459942145e-31 +5851,2240,2749,1.6772397283662627e-17 +5901,2653,2864,0.0003015718221877463 +5951,2446,3182,0.0005946478430517058 +6001,3873,2082,0.0035099937981355705 +6051,3067,2134,1.4269228740338587e-18 +6101,3579,2796,0.044070159340087074 +6151,3236,2296,1.0207282058200817e-10 +6201,2219,3070,6.128444533922714e-20 +6251,2764,4226,4.0114938634006145e-07 +6301,3819,3777,4.523296332293788e-28 +6351,4201,2117,0.0006366210109203425 +6401,3334,3416,0.0024069172689109197 +6451,4417,2689,0.0001510795635284454 +6501,2576,2211,1.464394314932261e-59 +6551,2008,4401,1.1960349880178966e-07 +6601,3749,4053,2.3049912620301826e-23 +6651,2744,3555,0.00039048361950001163 +6701,2181,3311,3.0049171175443936e-33 +6751,4664,4392,4.673341342704634e-75 +6801,2960,4656,2.1444652549020045e-07 +6851,2690,4656,0.05748086029092399 +6901,4668,4276,7.549019391941723e-59 +6951,3838,4516,1.0082038188372868e-28 +7001,2920,4477,0.05455036043254638 +7051,3857,2152,2.004089851444536e-29 +7101,3116,2551,2.600550198892281e-38 +7151,4953,3513,1.3795033294064868e-21 +7201,3592,4002,0.0016733670587886574 +7251,3244,2572,5.114970910833606e-38 +7301,4591,2744,0.09630014560631561 +7351,3708,2864,1.850482626161518e-12 +7401,2509,4772,7.605121039219745e-05 +7451,5134,5098,1.9026070284002625e-97 +7501,4970,4572,8.417077321620368e-55 +7551,3271,2510,5.027720092225525e-56 +7601,3518,2986,4.747697469921551e-21 +7651,2700,2458,1.5342703779996893e-108 +7701,4545,2345,3.6554706456314347e-22 +7751,4714,2780,3.593455829002233e-05 +7801,2588,2619,4.290089679393479e-115 +7851,4465,3125,0.002296908375565855 +7901,3699,4131,0.5129622712973023 +7951,5122,4697,2.8201026369032667e-44 +8001,3265,2821,1.2548891947495432e-59 +8051,2593,3386,4.265167768158848e-72 +8101,3279,3908,2.925977657425333e-14 +8151,4673,4713,1.1393802442050531e-20 +8201,3390,3706,2.6510169444939377e-19 +8251,3073,5077,0.006392608069466881 +8301,2679,4804,3.445753907043592e-19 +8351,5106,5102,2.4696466905259013e-42 +8401,3392,3776,2.188930822752493e-23 +8451,4521,3457,2.325359213936893e-05 +8501,3811,3220,3.675484972066898e-33 +8551,4138,2717,6.2817846465740635e-50 +8601,5331,2845,4.01855250200448e-10 +8651,4390,3656,8.142102766041389e-07 +8701,4620,5949,7.199420468330689e-38 +8751,4885,2749,1.232938976024364e-28 +8801,4155,4435,0.10455510731861506 +8851,4553,5067,2.8097763112411485e-08 +8901,3658,3129,3.388296324761072e-65 +8951,4091,6109,3.676828438505e-14 +9001,3374,5955,0.7977941668959887 +9051,4684,4369,0.9881371290166111 +9101,5941,4399,7.830542740608026e-16 +9151,5396,4205,0.006102055167403756 +9201,6087,5299,3.5473799808013304e-51 +9251,5364,3987,1.0 +9301,5605,5156,1.2301526743633959e-24 +9351,3703,3405,7.065360404513176e-69 +9401,5815,6427,2.27403687143121e-82 +9451,4567,5069,0.2128530705427019 +9501,6109,4587,2.098685170824773e-14 +9551,6501,5064,4.208356112420144e-40 +9601,6032,6370,1.7234955440937876e-79 +9651,6043,3510,0.001588282559753373 +9701,5350,4456,0.6461087776569374 +9751,3169,4579,1.8055731034783026e-58 +9801,3954,6268,0.270043231915957 +9851,3086,6734,1.8497805574159475e-07 +9901,4377,4383,6.979118182818968e-17 +9951,3050,3722,6.099767234370392e-137 diff --git a/validation/gwas/method/hwe/hwe_unit_test.ipynb b/validation/gwas/method/hwe/hwe_unit_test.ipynb new file mode 100644 index 000000000..704fe5385 --- /dev/null +++ b/validation/gwas/method/hwe/hwe_unit_test.ipynb @@ -0,0 +1,1577 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### HWE Unit Test Development" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import ctypes\n", + "import numpy as np\n", + "import pandas as pd\n", + "from sgkit.stats.hwe import hardy_weinberg_p_value_jit as hwep" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import glob\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['data/sim_01.csv']" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "glob.glob(str(Path('data')/ '*'))" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "rm -f *.o *.so \n", + "gcc -c -Wall -Werror -fpic chwe.c\n", + "gcc -shared -o libchwe.so chwe.o\n" + ] + } + ], + "source": [ + "!make clean\n", + "!make" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "libc = ctypes.CDLL(\"./libchwe.so\")\n", + "chwep = libc.hwep\n", + "chwep.restype = ctypes.c_double" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.8422797565707926" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "chwep(57, 14, 50)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.8422797565707925" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hwep(57, 14, 50)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Scalar Tests" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1, 0, 0],\n", + " [ 51, 27, 26],\n", + " [101, 47, 56],\n", + " [151, 71, 99],\n", + " [201, 137, 91]])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rs = np.random.RandomState(0)\n", + "n, s = 10_000, 50\n", + "n_het = np.expand_dims(np.arange(n, step=s) + 1, -1)\n", + "frac = rs.uniform(.3, .7, size=(n // s, 2))\n", + "n_hom = frac * n_het\n", + "n_hom = n_hom.astype(int)\n", + "args = np.concatenate((n_het, n_hom), axis=1)\n", + "args[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import time\n", + "df = []\n", + "for i in range(args.shape[0]):\n", + " a = [int(x) for x in args[i]]\n", + " p1 = chwep(a[0], a[1], a[2])\n", + " p2 = hwep(a[0], a[1], a[2])\n", + " df.append(dict(\n", + " p_true=p1,\n", + " p_pred=p2\n", + " ))\n", + "df = pd.DataFrame(df)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
p_truep_pred
01.0000001.000000
10.8459270.845927
20.8893040.889304
30.3684870.368487
40.2834420.283442
\n", + "
" + ], + "text/plain": [ + " p_true p_pred\n", + "0 1.000000 1.000000\n", + "1 0.845927 0.845927\n", + "2 0.889304 0.889304\n", + "3 0.368487 0.368487\n", + "4 0.283442 0.283442" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df.head()" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD4CAYAAADrRI2NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAATX0lEQVR4nO3df5BlZX3n8ffHgQgaXaBoyESYjFrEBN04aDNxy2RLQVzARKCMSayNobImY7KypZtYCWJqg1u1VWSjIUm5ZRwjlfFHkoVFhDUxcWCDllUK9uiATEYXNaMCU0xrKoUYAwt89497Olx6+sftnjn33uZ5v6pu3XOee3585s7Mt0+f55znpKqQJLXjKZMOIEkaLwu/JDXGwi9JjbHwS1JjLPyS1JhjJh1gFCeffHJt3bp10jEkaUPZs2fPt6pqZnH7hij8W7duZW5ubtIxJGlDSfL1pdp7P9WTZFOSLyT5WDd/UpLdSe7u3k/sO4Mk6XHjOMf/ZmD/0PzlwC1VdQZwSzcvSRqTXgt/ktOAVwF/MtR8EbCrm94FXNxnBknSE/V9xP8HwG8Cjw21nVpVBwG691OWWjHJjiRzSebm5+d7jilJ7eit8Cf5KeBQVe1Zz/pVtbOqZqtqdmbmsE5pSdI69XlVz0uBVye5EDgOeGaSDwH3J9lcVQeTbAYO9ZhBkrRIb0f8VfW2qjqtqrYCPw/8n6r6BeAm4NJusUuBG/vKIEk63CTu3L0KOC/J3cB53bwkaUzGcgNXVd0K3NpNfxs4dxz7lSQdbkPcuStJLdl6+V/+y/SBq1511LfvIG2S1BgLvyQ1xsIvSY2x8EtSYyz8ktQYC78kNcbCL0mNsfBLUmMs/JLUGAu/JDXGwi9JjbHwS1JjLPyS1BgLvyQ1xsIvSY2x8EtSYyz8ktSY3gp/kuOS3J7kjiT7kryja78yyb1J9navC/vKIEk6XJ+PXnwIOKeqHkxyLPDpJB/vPru6qt7Z474lScvorfBXVQEPdrPHdq/qa3+SpNH0eo4/yaYke4FDwO6quq376LIkdya5JsmJy6y7I8lckrn5+fk+Y0pSU3ot/FX1aFVtA04Dtid5AfAe4LnANuAg8K5l1t1ZVbNVNTszM9NnTElqyliu6qmqfwRuBc6vqvu7HwiPAe8Dto8jgyRpoM+remaSnNBNHw+8AvhSks1Di10C3NVXBknS4fq8qmczsCvJJgY/YK6tqo8l+WCSbQw6eg8Ab+wxgyRpkT6v6rkTOGuJ9tf3tU9J0uq8c1eSGmPhl6TGWPglqTEWfklqjIVfkhpj4Zekxlj4JakxFn5JaoyFX5IaY+GXpMZY+CWpMRZ+SWqMhV+SGmPhl6TGWPglqTEWfklqjIVfkhrT5zN3j0tye5I7kuxL8o6u/aQku5Pc3b2f2FcGSdLh+jzifwg4p6peCGwDzk/yEuBy4JaqOgO4pZuXJI1Jb4W/Bh7sZo/tXgVcBOzq2ncBF/eVQZJ0uF7P8SfZlGQvcAjYXVW3AadW1UGA7v2UPjNIkp6o18JfVY9W1TbgNGB7kheMum6SHUnmkszNz8/3F1KSGjOWq3qq6h+BW4HzgfuTbAbo3g8ts87OqpqtqtmZmZlxxJSkJvR5Vc9MkhO66eOBVwBfAm4CLu0WuxS4sa8MkqTDHdPjtjcDu5JsYvAD5tqq+liSzwDXJnkD8A3gtT1mkCQt0lvhr6o7gbOWaP82cG5f+5Ukrcw7dyWpMRZ+SWqMhV+SGmPhl6TGWPglqTEWfklqjIVfkhpj4Zekxlj4JakxFn5JaoyFX5IaY+GXpMZY+CWpMRZ+SWqMhV+SGmPhl6TGWPglqTEWfklqTJ8PWz89yd8m2Z9kX5I3d+1XJrk3yd7udWFfGSRJh+vzYeuPAL9RVZ9P8gxgT5Ld3WdXV9U7e9y3JGkZfT5s/SBwsJv+TpL9wLP62p8kaTRjOcefZCtwFnBb13RZkjuTXJPkxGXW2ZFkLsnc/Pz8OGJKUhNGKvxJXrDeHST5fuB64C1V9QDwHuC5wDYGvxG8a6n1qmpnVc1W1ezMzMx6dy9JWmTUI/4/TnJ7kv+Y5IRRN57kWAZF/8NV9RGAqrq/qh6tqseA9wHb15xakrRuIxX+qvoJ4N8DpwNzSf4syXkrrZMkwPuB/VX1+0Ptm4cWuwS4a82pJUnrNnLnblXdneS3gTngj4CzuuJ+xcLR/CIvBV4PfDHJ3q7tCuB1SbYBBRwA3ngE+SVJazRS4U/yY8AvAa8CdgM/3V2m+YPAZ4DDCn9VfRrIEpv7q/XHlSQdqVGP+N/N4Hz8FVX1vYXGqrqv+y1AkrRBjFr4LwS+V1WPAiR5CnBcVf1TVX2wt3SSpKNu1Kt6bgaOH5p/WtcmSdpgRi38x1XVgwsz3fTT+okkSerTqIX/u0letDCT5MXA91ZYXpI0pUY9x/8W4Lok93Xzm4Gf6yeSJKlPIxX+qvpckh8BnsfgEs0vVdX/6zWZJKkXaxmd82xga7fOWUmoqg/0kkqS1JtRb+D6IIOB1fYCj3bNBVj4JWmDGfWIfxY4s6qqzzCSpP6NelXPXcAP9BlEkjQeox7xnwz8XZLbgYcWGqvq1b2kkiT1ZtTCf2WfISRJ4zPq5ZyfTPJDwBlVdXOSpwGb+o0mSerDqI9e/BXgfwHv7ZqeBXy0r1CSpP6M2rn7JgYPVnkABg9lAU7pK5QkqT+jFv6HqurhhZkkxzC4jl+StMGMWvg/meQK4PjuWbvXAf97pRWSnJ7kb5PsT7IvyZu79pOS7E5yd/d+4pH9ESRJazFq4b8cmAe+yOAZuX8FrPbkrUeA36iqHwVeArwpyZndtm6pqjOAW7p5SdKYjHpVz2MMHr34vlE3XFUHgYPd9HeS7GfQKXwR8LJusV3ArcBvjZxYknRERh2r5+9Z4px+VT1nxPW3AmcBtwGndj8UqKqDSewklqQxWstYPQuOA14LnDTKikm+H7geeEtVPZBkpB0m2QHsANiyZcuIMSVp49p6+V+OZT8jneOvqm8Pve6tqj8AzlltvSTHMij6H66qj3TN9yfZ3H2+GTi0zD53VtVsVc3OzMyM9IeRJK1u1FM9LxqafQqD3wCesco6Ad4P7K+q3x/66CbgUuCq7v3GtQSWJB2ZUU/1vGto+hHgAPCzq6zzUuD1wBeT7O3armBQ8K9N8gbgGwxOG0mSxmTUq3pevtYNV9WnGTymcSnnrnV7kqSjY9RTPb++0ueLTuVIkqbYWq7qOZvB+XmAnwY+BXyzj1CSpP6s5UEsL6qq7wAkuRK4rqp+ua9gkqR+jDpkwxbg4aH5h4GtRz2NJKl3ox7xfxC4PckNDO7gvQT4QG+pJEm9GfWqnv+W5OPAT3ZNv1RVX+gvliSpL6Oe6gF4GvBAVf0hcE+SZ/eUSZLUo1Efvfg7DEbQfFvXdCzwob5CSZL6M+oR/yXAq4HvAlTVfawyZIMkaTqNWvgfrqqiG5o5ydP7iyRJ6tOohf/aJO8FTkjyK8DNrOGhLJKk6bHqVT3dKJv/E/gR4AHgecB/qardPWeTpCe9cY3BP2zVwl9VleSjVfViwGIvSRvcqKd6Ppvk7F6TSJLGYtQ7d18O/GqSAwyu7AmDXwZ+rK9gkqR+rFj4k2ypqm8AF4wpjySpZ6sd8X+UwaicX09yfVW9ZhyhJEn9We0c//ATtJ7TZxBJ0nisVvhrmelVJbkmyaEkdw21XZnk3iR7u9eFa9mmJOnIrXaq54VJHmBw5H98Nw2Pd+4+c4V1/xR4N4cP33x1Vb1zPWElSUduxcJfVZvWu+Gq+lSSretdX5LUj7UMy3y0XJbkzu5U0InLLZRkR5K5JHPz8/PjzCdJT2rjLvzvAZ4LbAMOAu9absGq2llVs1U1OzMzM658kvSkN9bCX1X3V9WjVfUYg0Heto9z/5KkMRf+JJuHZi8B7lpuWUlSP0YdsmHNkvw58DLg5CT3AL8DvCzJNgaXhh4A3tjX/iVJS+ut8FfV65Zofn9f+5MkjWYSV/VIkibIwi9JjbHwS1JjLPyS1BgLvyQ1prereiRJy5vEQ9YXeMQvSY2x8EtSYyz8ktQYC78kNcbCL0mNsfBLUmMs/JLUGAu/JDXGwi9JjbHwS1JjLPyS1BgLvyQ1prfCn+SaJIeS3DXUdlKS3Unu7t5P7Gv/kqSl9XnE/6fA+YvaLgduqaozgFu6eUnSGPVW+KvqU8A/LGq+CNjVTe8CLu5r/5KkpY37HP+pVXUQoHs/ZbkFk+xIMpdkbn5+fmwBJenJbmo7d6tqZ1XNVtXszMzMpONI0pPGuAv//Uk2A3Tvh8a8f0lq3rgL/03Apd30pcCNY96/JDWvz8s5/xz4DPC8JPckeQNwFXBekruB87p5SdIY9faw9ap63TIfndvXPiVp2k3yIesLprZzV5LUDwu/JDXGwi9JjbHwS1JjLPyS1BgLvyQ1xsIvSY2x8EtSYyz8ktQYC78kNcbCL0mNsfBLUmMs/JLUGAu/JDXGwi9JjeltPH5J2ogWxss/cNWrjvo2p4VH/JLUmIkc8Sc5AHwHeBR4pKpmJ5FDklo0yVM9L6+qb01w/5LUJE/1SFJjJnXEX8AnkhTw3qrauXiBJDuAHQBbtmwZczxJOnLT1qm7YFJH/C+tqhcBFwBvSvJvFy9QVTuraraqZmdmZsafUJKepCZS+Kvqvu79EHADsH0SOSSpRWMv/EmenuQZC9PAK4G7xp1Dklo1iXP8pwI3JFnY/59V1V9PIIckNWnshb+qvga8cNz7lSQNOGSDpKk1fFXMWodQWDz0wpFsay372Qi8jl+SGmPhl6TGWPglqTEWfklqjJ27ko6KtXZyjmvc+5WGTVgpw1qHW+ir87gPHvFLUmMs/JLUGAu/JDXGwi9JjbFzV5oSS3UmTrKTcNruSD2aeRZ/10d73PxpHYd/gUf8ktQYC78kNcbCL0mNsfBLUmOe9J27G+luuqWst0NrXB1zo36/feZZKsNa/96XyjfKXZ2j/JlHzTDqNkbpBD6aHcWrdVSutO9Rcq3l8/UuqyfyiF+SGmPhl6TGTKTwJzk/yZeTfCXJ5ZPIIEmtGnvhT7IJ+B/ABcCZwOuSnDnuHJLUqkkc8W8HvlJVX6uqh4G/AC6aQA5JalKqarw7TH4GOL+qfrmbfz3w41V12aLldgA7utnnAV9eZpMnA9/qKW4fNlpeMPM4bLS8YOZxOZLMP1RVM4sbJ3E5Z5ZoO+ynT1XtBHauurFkrqpmj0awcdhoecHM47DR8oKZx6WPzJM41XMPcPrQ/GnAfRPIIUlNmkTh/xxwRpJnJ/k+4OeBmyaQQ5KaNPZTPVX1SJLLgL8BNgHXVNW+I9jkqqeDpsxGywtmHoeNlhfMPC5HPfPYO3clSZPlnbuS1BgLvyQ1ZsMU/iSvTbIvyWNJDru0KcmWJA8meetQ24uTfLEbGuKPkix1KenYMyc5L8meLtueJOdMQ+aVvuMkb+syfTnJv5uGvIsl2Zbks0n2JplLsn3osyXzT4Mk/6nLtS/Jfx9qn9rMAEnemqSSnDzUNnWZk/xeki8luTPJDUlOGPps6vIu6HVom6raEC/gRxncyHUrMLvE59cD1wFvHWq7Hfg3DO4d+DhwwTRkBs4CfrCbfgFw7zRkXiHvmcAdwFOBZwNfBTZNOu8S+T+xsH/gQuDW1fJP+gW8HLgZeGo3f8q0Z+7ync7gAo2vAydPc2bglcAx3fTvAr87zXm7bJu6PM8Bvq/LeebR2v6GOeKvqv1VteTdu0kuBr4G7Btq2ww8s6o+U4Nv8gPAxWMJ21kuc1V9oaoW7l3YBxyX5KmTzrzCd3wR8BdV9VBV/T3wFWD7pPMuoYBndtP/isfvD1ky/wTyLeXXgKuq6iGAqjrUtU9zZoCrgd/kiTdfTmXmqvpEVT3SzX6Wwb1DMKV5O70ObbNhCv9ykjwd+C3gHYs+ehaDm8UW3NO1TZvXAF/o/uNPa+ZnAd8cml/INW153wL8XpJvAu8E3ta1L5d/Gvww8JNJbkvyySRnd+1TmznJqxn8lnrHoo+mNvOQ/8DgN1OY7ry9ZpuqJ3AluRn4gSU+entV3bjMau8Arq6qBxedXh5paIgjtc7MC+s+n8Gvnq9caFpisaOaeZ15l8s1lu/4CUFWyA+cC/znqro+yc8C7wdewQRyDlsl8zHAicBLgLOBa5M8h+nOfAWP/5t9wmpLtI0l8yj/rpO8HXgE+PDCakssPy3Xt/eabaoKf1W9Yh2r/TjwM12n2AnAY0n+mcE5/9OGlutlaIh1ZibJacANwC9W1Ve75nvoOfM68y43zEbveRdbKX+SDwBv7mavA/6km57oMCGrZP414CPdqbLbkzzGYFCuqcyc5F8zOB9+R3egdRrw+a4jfWKZV/t3neRS4KeAc7vvGqZ7+Jh+s026E2MdnR63skTnbvfZlTyxc/dzDI6kFjoeL5yGzAx+QN0BvGaJZSeeeYm8z+eJnWBf4/HO3YnnHcq5H3hZN30usGe1/JN+Ab8K/Ndu+ocZ/Hqfac68KP8BHu/cncrMwPnA3wEzi9qnMm+X7Zguz7N5vHP3+Udt+5P+A67hi7iEwU/Bh4D7gb9ZYpnFhX8WuItB7/i76e5UnnRm4LeB7wJ7h16nTDrzSt8xg1/xv8pgeOwLhton+h0vyv8TwJ7uP8ltwItXyz/pV/ef+kPdd/h54Jxpz7wo/78U/mnNzKDT9ptD/9f+eJrzDmW7EPi/Xb63H81tO2SDJDVmw1/VI0laGwu/JDXGwi9JjbHwS1JjLPyS1BgLvyQ1xsIvSY35/+aPfCzCrGUxAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "df['p_true'].apply(np.log10).plot(kind='hist', bins=128)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEHCAYAAAC5u6FsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAe0klEQVR4nO3df5iUdb3/8ed7f7CQgBAoHlhWVKACRcxR4JD9UDTyW6JiRtnBc6rDycy072UoUX3TjpV4rjz9OnY41tXXssNBKfEcM4zSTBN1yQVBsxZRWSiM/YLuFgzs7vv7x9wDwzKzs/cyM/c9M6/Hdc3l3J/7np13G/ri8+P+3ObuiIiI9FdN1AWIiEh5UXCIiEgoCg4REQlFwSEiIqEoOEREJBQFh4iIhFIXdQGlMHr0aJ8wYULUZYiIlI3169fvcvfjsp2riuCYMGECzc3NUZchIlI2zOzlXOc0VCUiIqEoOEREJBQFh4iIhKLgEBGRUMoyOMxsrpm9YGatZnZj1PWIiFSTsltVZWa1wLeB84E24Gkzu9/dn4u2MhGReGjd2cGl//YYryd7GDd8EI9/9vyC/vxy7HGcDbS6+4vuvh9YAcyLuCYRkVi4Yvk65tz+KK8newDY/vp+Jtz4QEG/oxyDYxywLeO4LWgTEalazVvbOenGB3j8xfas52d/+ecF+66yG6oCLEvbEU+jMrNFwCKApqamYtckIhKZ9/zrIzz/p7/0ec321/cX7PvKscfRBozPOG4EdvS+yN2Xu3vC3RPHHZf1rnkRkbLWvLWdk298IG9oAIwbPqhg31uOPY6ngUlmdhKwHVgAfCjakkRESuuCrz3C71/NHxhphZwgL7vgcPcuM/sksAaoBb7n7psjLktEpCR+8dyfWHTXerpDfOalr/6vgtZQdsEB4O4/BX4adR0iIqXSurODD/z7b2j/a1e/P/POSaP5/kdnFLyWsgwOEZFqsvieFlau397v6+uBdZ+bw6ihDUWpR8EhIhJjn/vJs6FC418vn8bFbx2f/8KjoOAQEYmpC772ML9/9a/9unboINh0c2HnMnJRcIiIxEx7Z5ILv/4oOzvy33txzCDjGwveynlTTihBZSkKDhGRGLl73cv8n/s30dXT93UGXHveRK47/00lqSuTgkNEJGKtOzt4rPXPPP/HDv6ruS3v9RdP/xs+/96pRZv8zkfBISISoWv/87es3vDHfl9/eaKRZZedXsSK8lNwiIhE5FP/+Vvu72dojBxSyz0fn83EMcOKXFV+Cg4RkQj84rk/9Ss0zmwawSfeeUpJJ7/zUXCIiJRIe2eSJ7bs4vu/eZnml3f3eW19LXzxfadyxcwTS1Rd/yk4RERK4O51L/O5+zYd+QyILK56x8l87JyTI5v8zkfBISJSRK07O7jjkVZWPXPE0x+y+ux73syid5xS5KqOjoJDRKRIFt+7gZX9WF4LUGtw88WncsWM+A1N9abgEBEpsPbOJLf97Hf9Do0LTxvDl+adFtuhqd4UHCIiBbT8V1v48oO/69e1s08exU3zpsZiiW0YCg4RkQJo70yy6K5m1r+yp8/raoBLzhjLVe+cWHaBkabgEBE5Sv1dMVVj8NB1by/bwEhTcIiIDFB7Z5Jv/uIPfP+Jl/NeW1djfO3y08s+NEDBISIyIHeve5nP37eJPJvYAvD3s07kmvMmlc3kdz4KDhGRkG5/6AW+/svWfl17yyXlscQ2jFgGh5ndBrwP2A9sAf7B3fcE55YAHwW6gU+5+5rIChWRqpHaLqSdOx97kZZtr+W9/s1jjuHuf5xVMb2MTLEMDuDnwBJ37zKzW4ElwA1mNgVYAEwFxgJrzWyyu3dHWKuIVLi7173M0vs29eva2aeM4tNzJpE4aVSRq4pOLIPD3R/KOFwHXBa8nwescPcksNXMWoGzgSdKXKKIVIn+DkvVGHxpXjw3JSy0WAZHLx8B/it4P45UkKS1BW0iIgXVurODG3/8bN5dbAEWJBr5zNw3V+SwVDaRBYeZrQWybTC/1N1XB9csBbqAu9Mfy3J91qXTZrYIWATQ1NR01PWKSHVo70yy+N4N/OJ3f+7X9ZcnGvlqxE/kK7XIgsPd5/R13syuBN4LnOfu6XBoA8ZnXNYIZN1y0t2XA8sBEolEf3YyFpEqt7plO9euaOnXtWdNGMFXLplWEfdlhBXLoSozmwvcALzD3f+acep+4Edm9jVSk+OTgKciKFFEKkx7Z7JfoVFjcOPc+G99XkyxDA7gW0AD8HMzA1jn7h93981mthJ4jtQQ1tVaUSUiR6u9M8l3H3sx73Vz3nIct84/vWrmMnKJZXC4+8Q+zt0C3FLCckSkQrV3Jrn7yVf49sOt9HjuEe2zTxzBly+tzmGpbGIZHCIixba6ZTuL791AsqvvKdAPz2jiny85rURVlYeaqAsQESm19s4kN6zamDc0Zp40UqGRhXocIlI12juTtO3ey2t7D1BfU8O+LFsUDqqtYeGsE5k7dUxF3/19NBQcIlLRWnd28Fjrn9m666+sePoVBtXWsr+7h+6eI0Ojoa6G2y6bxkXTdV9xXxQcIlKRWnd28IXVm/nNi+2HtSe7ugCorzUa6giCpJtPvmsSH5rRVPUrpvpDwSEiFecL9z3LXete6fOawXW1fPuKMzh2yCAaRw5RYISg4BCRitK8tT1vaAAc6Olh6thjFRgDoFVVIlIxVrdsZ8GdT+a9rqHOWDZ/mkJjgNTjEJGy196ZZPOO11h87wa6unMvsTXgf58/WXMZR0nBISJlbXXLdm5YtZEas5z3ZZx14giu/NuTmHXKKAVGASg4RKQsNW9tZ81zO/n+41s5cOTKWgDqao0VH5uh+zEKTMEhImWlvTPJR77/NBvacj/3+w31tfTgLJs/TaFRBAoOESkbq1u2c/3Klpw9DEjdxPedvzuTqWOHa1iqSBQcIlIW0vtL5QqN+lqjtia1Wurtk48rbXFVRsEhImWhbffenPtLAdx80VQumHqCehkloOAQkbLQOHIIB7LsLwVwzsRRfHDGiSWuqHrpBkARKQujhjawbP40BtfXMKyhjvoamDtlDPf+00x+8LGZUZdXVdTjEJFYaO9M8tDmP/Hirr/w7inZtzS/aPo4Zk8cTdvuvdpfKkIKDhGJ3OqW7Vy7ouXg8X/8eivnTByVtScxamiDAiNiGqoSkUi1dya5fmXLEe2/bm2neWt7lk9I1GIdHGZ2vZm5mY3OaFtiZq1m9oKZvTvK+kRk4No7k2zYtofNO14Hs6zXPPqHXSWuSvojtkNVZjYeOB94JaNtCrAAmAqMBdaa2WR3746mShEZiPT+UvU1Namn8eXYmPDtk0ZnbZdoxbnHcTuwGMj8EzUPWOHuSXffCrQCZ0dRnIgMTOvODj5z70b2HeihI9lFsquH2tojexznTByl7UJiKpY9DjO7CNju7hvs8C7sOGBdxnFb0CYiMdbemaRt9142bX+Nm/57M/t79TAG19Xy3SvPYPvuvX2uqpJ4iCw4zGwtcEKWU0uBzwIXZPtYlrasfVwzWwQsAmhqahpglSJytNLDUnU1Rmcy+6hy+ml8b598fImrk4GILDjcfU62djM7DTgJSPc2GoHfmtnZpHoY4zMubwR25Pj5y4HlAIlEIveTXUSkaNL7S+3rY1fCQXU1ehpfmYndUJW7Pwsc/GuHmb0EJNx9l5ndD/zIzL5GanJ8EvBUJIWKSE6tOzto2baHwfW1fe4vNajW+Ok1b2PimGElrlCORuyCoy/uvtnMVgLPAV3A1VpRJRIvX7jvWe5ad3AxJLU1R44wHzOolm5PPS9DoVF+Yh8c7j6h1/EtwC3RVCMifWnd2XFYaAB09zgNdTUMqq3hQE8Pn3/vFE4de6y2DCljsQ8OESkfLdv2ZG1f8p43c0bTSIVFhVBwiEjBTB8/Imv72yaO1pBUBYnzDYAiUmYmjhnGwlmHL39fOKtJoVFh1OMQkYK6ed5pLJw5gZZte5g+foRCowIpOESk4CaOGabAqGAaqhKRPqV3sW3vTEZdisSEehwiklPmLrYHenpYNn8aF03X9nDVTj0OETlCe2eSR3//Kovv3XBwF9t9B3pYvGqjeh6iHoeIHC7dy6gxI9l1+DZv9TU1tO3eq3sxqpyCQ0QOyrcp4YGeHhpHDilxVRI3Cg4ROaht996smxK+ob6WHly72Aqg4BCRDI0jh3Cg5/DQaKir4Tt/dyZTxw5XaAigyXGRqpRrie2ooQ0smz+NwfU1DGuoY3B9DbddNo23Tz5OoSEHqcchUmXyLbG9aPo4Zk8cTdvuvdqUULJScIhUkdadHXzmng3s7/aD8xiLV21k9sTRhwXEqKENCgzJSUNVIlVidct2LvzmY+zvzr7EVqS/FBwiVSC9zHZ/15HLbLXEVsLSUJVIBWvvTNK2ey+v7d2fdZntoFrTElsJTcEhUmHaO5Ns3vE6T2zZxfce38qg2lr2d3fTc/gIFYPqavjpNW/TLrYSmoJDpIIs/9UWvvKz3+EZIZHs6gKgrobDnv29bP40hYYMSGyDw8yuAT4JdAEPuPvioH0J8FGgG/iUu6+JrkqR+Fh87wZWNrflPD+kvo5vX/FWjh1Sr2W2clRiGRxm9i5gHjDN3ZNmdnzQPgVYAEwFxgJrzWyyu3dHV61I9Fp3dvQZGpCaBNfd31IIeYPDzDoAz3Xe3YcXtKKUq4Cvunsy+I5Xg/Z5wIqgfauZtQJnA08UoQaRstGybU+f5xvqNAkuhZM3ONx9GICZ3Qz8CfgBYMAVQLEGSCcD55jZLcA+4Hp3fxoYB6zLuK4taBOpatPHj8jaXldjXHveJD40o0mhIQUTZqjq3e4+I+P4DjN7Elg2kC82s7XACVlOLQ3qGgnMBM4CVprZyaQCq7esvSEzWwQsAmhqahpIiSJlY+KYYSyc1cRdT7xysG3OW47nVvUypAjCBEe3mV0BrCD1H+sPkpqgHhB3n5PrnJldBfzY3R14ysx6gNGkehjjMy5tBHbk+PnLgeUAiUQi51CbSKW4ed5pLJw5gZZte5g+foRWTEnRhLlz/EPA5cDO4PX+oK0Y7gPOBTCzycAgYBdwP7DAzBrM7CRgEvBUkWoQKTsTxwzjssR4hYYUVb97HO7+EqnJ6VL4HvA9M9sE7AeuDHofm81sJfAcqWW6V2tFlVSq9F3fWjorcdPv4Aj+5n8HMMbdTzWzacBF7v7PhS7K3fcDH85x7hbglkJ/p0ic5Nv6XCRKYYaq/gNYAhwAcPeNpO6pEJECynzud0eyi30Heli8auMRD10SiUqY4HiDu/eeT+gqZDEicui535m09bnESZjg2GVmpxAsfzWzy4A/FqUqkSqW7bnf2vpc4iRMcFwN/DvwZjPbDlwHfLwoVYlUsWzP/dZd3xIn/ZocN7Na4Cp3n2NmxwA17t5R3NJEKltfq6b03G+Js34Fh7t3m9mZwfu/FLckkcrXn1VTeu63xFWYO8efMbP7gXuAg+Hh7j8ueFUiFSxz1VT6iXyLV21k9sTRCgopC2GC441AO8Ed3QEHFBwiIaRXTWU+xjW9akrBIeUgzJ3j/1DMQkQqUbZ5DK2aknLX71VVZnaymf23mf3ZzF41s9XBflEiksXqlu3MvvWXfPjOJ5l96y+5v2U7oFVTUv7CDFX9CPg2cElwvIDUTrkzcn5CpErlm8fQqikpZ2Hu4zB3/4G7dwWvH9LHkwFFqll/7v4eNbSB08ePUGhI2QnT43jYzG7k0PM4PgA8YGZvBHD3/1eE+kTKkuYxpJKFCY4PBP/8p17tHyEVJCcXpCKRMpNtAjw9j7G4170a6l1IJQizqqrPiXAzO9/df370JYmUj75u5NM8hlSqMHMc+dxawJ8lEnv92f5c8xhSiQoZHFbAnyUSe9r+XKpVIYNDK6ykqmgCXKpVIYNDpKroRj6pVmGeOT4Y+ATwNlK9i8eAO9x9X3DJSwWvTiRCfW17nqYJcKlGYZbj3gV0AN8Mjj8I/AB4P4C7X1qoosxsOvAdYDCpx9N+Iv3YWjNbAnwU6AY+5e5rCvW9Imn92fY8TdufS7UJExxvcvfTM44fNrMNhS4osAy4yd0fNLMLg+N3mtkUUludTAXGAmvNbLK7dxepDqlC2vZcpG9h5jieMbOZ6QMzmwE8XviSgNRQ2PDg/bHAjuD9PGCFuyfdfSvQCpxdpBqkSmm1lEjfwvQ4ZgALzeyV4LgJeN7MngXc3acVsK7rgDVm9i+kwu1vg/ZxwLqM69qCNpEB6z2XodVSIn0LExxzC/nFZrYWOCHLqaXAecCn3X2VmV0OfBeYQ/Z7RbIuAzazRcAigKampoLULJUn11yGtgsRyc3c43f7hZm9BoxwdzczA15z9+HBxDju/pXgujXAF939ib5+XiKR8Obm5qLXLeWlvTPJ7Ft/yb4Dh3oXg+trePyGcxk1tKFfq6pEKpWZrXf3RLZzcb2PYwfwjuD9ucAfgvf3AwvMrCF4iNQk4KkI6pMKkG8uQ9uFiGQXZqiqlP4R+LqZ1QH7CIac3H2zma0EniO1TPdqraiSgdJchsjAxLLH4e6PufuZ7n66u89w9/UZ525x91Pc/U3u/mCUdUp5053fIgMT1x6HSEnozm+R8BQcUvV057dIOLEcqhIRkfhScIiISCgKDhERCUXBISIioSg4REQkFAWHiIiEouAQEZFQFBxSNto7k2zYtof2zmTUpYhUNd0AKGUhzKNcRaS41OOQ2Mt8lGtHsot9B3pYvGqjeh4iEVFwSOzpUa4i8aLgkNjT9uci8aLgkNjT9uci8aLJcSkL2v5cJD4UHBKZsM/01vbnIvGg4JBIaHmtSPnSHIeUnJbXipQ3BYeUnJbXipS3yILDzN5vZpvNrMfMEr3OLTGzVjN7wczendF+ppk9G5z7hplZ6SuXo6XltSLlLcoexybgUuDRzEYzmwIsAKYCc4F/M7Pa4PQdwCJgUvCaW7JqpWC0vFakvEU2Oe7uzwNk6TTMA1a4exLYamatwNlm9hIw3N2fCD53F3Ax8GDJipaC0fJakfIVx1VV44B1GcdtQduB4H3vdilTWl4rUp6KGhxmthY4Icuppe6+OtfHsrR5H+25vnsRqWEtmpqa8lQqIiL9VdTgcPc5A/hYGzA+47gR2BG0N2Zpz/Xdy4HlAIlEImfAiIhIOHFcjns/sMDMGszsJFKT4E+5+x+BDjObGaymWgjk6rWIiEiRRLkc9xIzawNmAQ+Y2RoAd98MrASeA34GXO3u3cHHrgLuBFqBLWhiXESk5My98kdxEomENzc3R12GiEjZMLP17p7Idi6OQ1UiIhJjCg4REQlFwSEiIqEoOCSr9s4kG7bt0Y61InKEON45LhHTszJEpC/qcchh9KwMEclHwSGH0bMyRCQfBYccRs/KEJF8FBxyGD0rQ0Ty0eS4HEHPyhCRvig4JCs9K0NEctFQlYiIhKLgEBGRUBQcFU53gItIoWmOo4LpDnARKQb1OCqU7gAXkWJRcFQo3QEuIsWi4KhQugNcRIpFwVGhdAe4iBSLJscrmO4AF5FiiKzHYWbvN7PNZtZjZomM9vPNbL2ZPRv889yMc2cG7a1m9g0zs2iqLx+jhjZw+vgRCg0RKZgoh6o2AZcCj/Zq3wW8z91PA64EfpBx7g5gETApeM0tQZ0iIpIhsqEqd38eoHenwd2fyTjcDAw2swbgjcBwd38i+NxdwMXAgyUpWEREgPhPjs8HnnH3JDAOaMs41xa0iYhICRW1x2Fma4ETspxa6u6r83x2KnArcEG6Kctl3sfnF5Ea1qKpqalf9YqISH5FDQ53nzOQz5lZI/ATYKG7bwma24DGjMsagR19fPdyYDlAIpHIGTAiIhJO7IaqzGwE8ACwxN0fT7e7+x+BDjObGaymWgj02WsREZHCi3I57iVm1gbMAh4wszXBqU8CE4HPm1lL8Do+OHcVcCfQCmxBE+MiIiVn7pU/ipNIJLy5uTnqMkREyoaZrXf3RLZzsRuqEhGReFNwiIhIKAoOEREJRcEhIiKhKDhERCQUBYeIiISi4BARkVAUHCIiEoqCQ0REQlFwiIhIKAoOEREJRcEhIiKhKDhERCQUBYeIiISi4BARkVAUHCIiEoqCow/tnUk2bNtDe2cy6lJERGKjLuoC4mp1y3ZuWLWR+poaDvT0sGz+NC6aPi7qskREIqceRxbtnUluWLWRfQd66Eh2se9AD4tXbVTPQ0QEBUdWbbv3Ul9z+K+mvqaGtt17I6pIRCQ+IgsOM3u/mW02sx4zO+KB6GbWZGadZnZ9RtuZZvasmbWa2TfMzIpRW+PIIRzo6Tms7UBPD40jhxTj60REykqUPY5NwKXAoznO3w482KvtDmARMCl4zS1GYaOGNrBs/jQG19cwrKGOwfU1LJs/jVFDG4rxdSIiZSWyyXF3fx4gW6fBzC4GXgT+ktH2N8Bwd38iOL4LuJgjw6UgLpo+jtkTR9O2ey+NI4coNEREArFbVWVmxwA3AOcD12ecGge0ZRy3BW1FM2pogwJDRKSXogaHma0FTshyaqm7r87xsZuA2929s1dvJNt8hvfx3YtIDWvR1NTUv4JFRCSvogaHu88ZwMdmAJeZ2TJgBNBjZvuAVUBjxnWNwI4+vns5sBwgkUjkDBgREQkndkNV7n5O+r2ZfRHodPdvBccdZjYTeBJYCHwzkiJFRKpYlMtxLzGzNmAW8ICZrenHx64C7gRagS0UaWJcRERyM/fKH8Uxsz8DL2c5NRrYVeJyjpZqLr5yqxdUcymUW71wdDWf6O7HZTtRFcGRi5k1u/sRNx/GmWouvnKrF1RzKZRbvVC8mrXliIiIhKLgEBGRUKo9OJZHXcAAqObiK7d6QTWXQrnVC0WquarnOEREJLxq73GIiEhICg4REQmlaoIjzs//CFOvmZ1vZuuDutab2blxqLevmoNzS4K6XjCzd8el5l41TjezdWbWYmbNZnZ2xrms9ceBmV0T1LU52Kon3R7nmq83Mzez0RltsazXzG4zs9+Z2UYz+4mZjcg4F8uaAcxsblBXq5ndWNAf7u5V8QLeArwJeARIZDm/CrgHuD6j7SlSd7YbqbvU3xN1vcAZwNjg/anA9jjUm6fmKcAGoAE4idRd/7VxqLlX/Q+lvx+4EHgkX/1Rv4B3AWuBhuD4+DKoeTywhtRNuaPLoN4LgLrg/a3ArWVQc21Qz8nAoKDOKYX6+VXT43D35939hWznMp7/sTmj7eDzPzz1/0T6+R8lkated3/G3dObO24GBptZQ9T1BrXl+h3PA1a4e9Ldt5LaMubsONTciwPDg/fHcmgTzaz1R1BfNlcBX3X3JIC7vxq0x7nm24HFHL67dWzrdfeH3L0rOFzHoc1WY1szqTpa3f1Fd98PrCBVb0FUTXDkkvH8j5t6nSr58z8GYD7wTPAfjTjXOw7YlnGcri1uNV8H3GZm24B/AZYE7bnqj4PJwDlm9qSZ/crMzgraY1mzmV1Eqpe8odepWNabxUc4tEdenGsuam2x2x33aET5/I+BGGC96c9OJdVtviDdlOWygq+1HmDNuWorSc2HFdJH/cB5wKfdfZWZXQ58F5hDBHVmylNzHTASmAmcBaw0s5OJsOY89X6WQ39mD/tYlrZY/I7Tf67NbCnQBdyd/liW6+Nyf0NRa6uo4PAIn/8xEAOsFzNrBH4CLHT3LUFzG0WuFwZccxupce20dG0lqTlTX/UHjyO+Nji8h9ROzJC7/pLIU/NVwI+Dob6nzKyH1MZ2kdWcq14zO43UXMCG4C9pjcBvg0UIsf0dA5jZlcB7gfOC3zVEXHMexa0t6kmcUr/IMTkenPsih0+OP03qb3LpidsLo66XVLhtAOZnuTbyenPUPJXDJxFf5NDkeCxqDmp5Hnhn8P48YH2++qN+AR8Hbg7eTyY1PGFxrjmj9pc4NDke23qBucBzwHG92uNcc11Qz0kcmhyfWrCfH/X/wBL+Ii8hlcJJYCewJss1vYMjAWwitTrhWwR32kdZL/A54C9AS8br+Kjrzfc7JjVMsQV4gYyVU1HX3Kv+twHrg3/JngTOzFd/1K/gPwo/DH6HvwXOjXvNGfUdDI4410tq0ntbxr9v34l7zUFtFwK/D+pbWsifrS1HREQklKpfVSUiIuEoOEREJBQFh4iIhKLgEBGRUBQcIiVkZteZ2RuirkPkaGhVlUgJmdlLpO5x2ZXlXK27d5e+KpFw1OMQOUpmNiHYdvv/Bltv35utV2FmnwLGAg+b2cNBW6eZ3WxmTwKzzOyl9FbjZpYws0eC98eY2ffM7Gkze8bMCrZhnUhYCg6RwngTsNzdpwGvA5/ofYG7f4PUtg/vcvd3Bc3HAJvcfYa7P9bHz18K/NLdzyK1lfptwQadIiWn4BApjG3u/njw/oek7kLvj25S+6LlcwFwo5m1kNrSZTDQFLZIkUKoqE0ORSLUe7Kwv5OH+3rNa3Rx6C90gzPajdT+ZFmfKSNSSupxiBRGk5nNCt5/EMg17NQBDOvj57wEnBm8n5/Rvga4Jv1oXTM7Y+ClihwdBYdIYTwPXGlmG4E3AnfkuG458GB6cjyLm4Cvm9mvSQ1jpX0JqAc2mtmm4FgkElqOK3KUzGwC8D/ufmrEpYiUhHocIiISinocIkVgZj8h9RCdTDe4+5oo6hEpJAWHiIiEoqEqEREJRcEhIiKhKDhERCQUBYeIiISi4BARkVAUHCIiEsr/B7cMUL4iUxY6AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "(\n", + " df\n", + " .assign(\n", + " p_true=lambda df: np.log10(df['p_true']),\n", + " p_pred=lambda df: np.log10(df['p_pred'])\n", + " )\n", + " .plot(kind='scatter', x='p_true', y='p_pred')\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "np.testing.assert_allclose(df['p_true'], df['p_pred'])" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1.00000000e+000, 8.45926829e-001, 8.89304245e-001, 3.68487492e-001,\n", + " 2.83442131e-001, 1.93780506e-001, 3.46415612e-002, 9.77805142e-007,\n", + " 9.00169099e-002, 2.77392776e-004, 5.78595078e-006, 1.56290046e-001,\n", + " 3.11983705e-002, 7.78234779e-001, 6.28255056e-001, 9.17242816e-001,\n", + " 8.81087089e-001, 1.20954751e-004, 6.51960684e-002, 4.87927509e-007,\n", + " 6.14320396e-002, 1.67216769e-003, 2.58323982e-003, 9.22666204e-012,\n", + " 1.15591803e-003, 1.00000000e+000, 5.21303203e-001, 2.40595832e-012,\n", + " 1.79017126e-001, 8.50964237e-004, 4.08782584e-018, 2.65625649e-003,\n", + " 1.73047163e-007, 2.61257337e-002, 3.40282167e-002, 5.57265342e-006,\n", + " 2.28187711e-010, 3.71009969e-005, 2.02796027e-015, 2.85690782e-015,\n", + " 4.43715904e-004, 1.24880234e-005, 1.39680904e-002, 6.69133747e-009,\n", + " 9.43219724e-010, 6.10161450e-001, 1.93499955e-003, 1.44451527e-014,\n", + " 1.15651799e-011, 6.16416362e-006, 2.18519190e-001, 2.67902896e-020,\n", + " 3.81265044e-003, 1.87170429e-002, 2.87276124e-001, 1.46939801e-004,\n", + " 5.90523804e-001, 9.00712608e-003, 7.82143524e-011, 1.55029275e-016,\n", + " 1.00796610e-003, 6.51775272e-018, 7.22627291e-001, 3.50621941e-033,\n", + " 2.15694037e-001, 5.36554440e-001, 4.98209450e-023, 1.00725415e-002,\n", + " 2.83256119e-004, 2.31647615e-001, 5.40831311e-004, 2.28693251e-006,\n", + " 2.33943256e-016, 4.63666449e-002, 1.95571664e-029, 1.32013500e-001,\n", + " 1.93010279e-006, 1.72246817e-002, 4.44008208e-010, 2.64771353e-025,\n", + " 1.42567926e-002, 2.34658222e-023, 5.14985651e-044, 4.48467881e-038,\n", + " 2.38901290e-003, 3.00019737e-020, 9.91998679e-058, 3.85771324e-001,\n", + " 1.19901665e-004, 1.09586529e-012, 4.52696626e-007, 4.52117435e-005,\n", + " 3.74269466e-022, 1.84769664e-002, 9.01235925e-001, 4.71167421e-016,\n", + " 7.26213285e-001, 2.68067642e-005, 1.95763513e-027, 3.44681033e-030,\n", + " 6.72973257e-001, 1.90998085e-021, 2.71129678e-092, 1.33474542e-002,\n", + " 9.42328262e-016, 6.04559513e-002, 2.73568136e-002, 3.45497420e-013,\n", + " 1.85964309e-010, 2.25791165e-016, 8.88002002e-023, 7.31645858e-001,\n", + " 6.20103273e-001, 2.02013957e-003, 3.26543825e-041, 9.55096556e-034,\n", + " 1.58435946e-031, 1.67723973e-017, 3.01571822e-004, 5.94647843e-004,\n", + " 3.50999380e-003, 1.42692287e-018, 4.40701593e-002, 1.02072821e-010,\n", + " 6.12844453e-020, 4.01149386e-007, 4.52329633e-028, 6.36621011e-004,\n", + " 2.40691727e-003, 1.51079564e-004, 1.46439431e-059, 1.19603499e-007,\n", + " 2.30499126e-023, 3.90483620e-004, 3.00491712e-033, 4.67334134e-075,\n", + " 2.14446525e-007, 5.74808603e-002, 7.54901939e-059, 1.00820382e-028,\n", + " 5.45503604e-002, 2.00408985e-029, 2.60055020e-038, 1.37950333e-021,\n", + " 1.67336706e-003, 5.11497091e-038, 9.63001456e-002, 1.85048263e-012,\n", + " 7.60512104e-005, 1.90260703e-097, 8.41707732e-055, 5.02772009e-056,\n", + " 4.74769747e-021, 1.53427038e-108, 3.65547065e-022, 3.59345583e-005,\n", + " 4.29008968e-115, 2.29690838e-003, 5.12962271e-001, 2.82010264e-044,\n", + " 1.25488919e-059, 4.26516777e-072, 2.92597766e-014, 1.13938024e-020,\n", + " 2.65101694e-019, 6.39260807e-003, 3.44575391e-019, 2.46964669e-042,\n", + " 2.18893082e-023, 2.32535921e-005, 3.67548497e-033, 6.28178465e-050,\n", + " 4.01855250e-010, 8.14210277e-007, 7.19942047e-038, 1.23293898e-028,\n", + " 1.04555107e-001, 2.80977631e-008, 3.38829632e-065, 3.67682844e-014,\n", + " 7.97794167e-001, 9.88137129e-001, 7.83054274e-016, 6.10205517e-003,\n", + " 3.54737998e-051, 1.00000000e+000, 1.23015267e-024, 7.06536040e-069,\n", + " 2.27403687e-082, 2.12853071e-001, 2.09868517e-014, 4.20835611e-040,\n", + " 1.72349554e-079, 1.58828256e-003, 6.46108778e-001, 1.80557310e-058,\n", + " 2.70043232e-001, 1.84978056e-007, 6.97911818e-017, 6.09976723e-137])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df['p_true'].values" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "nan" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hwep(0, 0, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1000 0.9999999999999996 2.423934816528403e-34\n", + "10000 1.0 0.0\n", + "100000 0.9999999999999997 0.0\n", + "1000000 1.0 0.0\n", + "10000000 1.0 0.0\n" + ] + } + ], + "source": [ + "for n_het in 10**np.arange(3, 8):\n", + " p1 = hwep(n_het, n_het//2, n_het//2)\n", + " p2 = hwep(n_het, n_het//10, n_het//2 + n_het//10)\n", + " print(n_het, p1, p2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Vectorized Tests" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "Show/Hide data repr\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Show/Hide attributes\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
xarray.Dataset
    • alleles: 2
    • ploidy: 2
    • samples: 1000
    • variants: 50
      • variant/contig
        (variants)
        int64
        0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0
        array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
+       "       0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,\n",
+       "       0, 0, 0, 0, 0, 0])
      • variant/position
        (variants)
        int64
        0 1 2 3 4 5 6 ... 44 45 46 47 48 49
        array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,\n",
+       "       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,\n",
+       "       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49])
      • variant/alleles
        (variants, alleles)
        |S1
        b'T' b'C' b'T' ... b'A' b'T' b'C'
        array([[b'T', b'C'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'C', b'T'],\n",
+       "       [b'G', b'G'],\n",
+       "       [b'A', b'C'],\n",
+       "       [b'G', b'A'],\n",
+       "       [b'T', b'A'],\n",
+       "       [b'T', b'C'],\n",
+       "       [b'G', b'A'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'C', b'G'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'C', b'C'],\n",
+       "       [b'T', b'A'],\n",
+       "       [b'C', b'G'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'A', b'C'],\n",
+       "       [b'C', b'T'],\n",
+       "       [b'T', b'A'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'T', b'A'],\n",
+       "       [b'C', b'C'],\n",
+       "       [b'C', b'T'],\n",
+       "       [b'G', b'G'],\n",
+       "       [b'G', b'T'],\n",
+       "       [b'T', b'C'],\n",
+       "       [b'A', b'C'],\n",
+       "       [b'G', b'T'],\n",
+       "       [b'G', b'C'],\n",
+       "       [b'T', b'A'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'G', b'G'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'A', b'A'],\n",
+       "       [b'T', b'T'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'G', b'A'],\n",
+       "       [b'C', b'G'],\n",
+       "       [b'T', b'T'],\n",
+       "       [b'G', b'A'],\n",
+       "       [b'A', b'T'],\n",
+       "       [b'G', b'A'],\n",
+       "       [b'A', b'G'],\n",
+       "       [b'C', b'C'],\n",
+       "       [b'T', b'G'],\n",
+       "       [b'C', b'A'],\n",
+       "       [b'T', b'C']], dtype='|S1')
      • sample/id
        (samples)
        <U4
        'S0' 'S1' 'S2' ... 'S998' 'S999'
        array(['S0', 'S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10',\n",
+       "       'S11', 'S12', 'S13', 'S14', 'S15', 'S16', 'S17', 'S18', 'S19',\n",
+       "       'S20', 'S21', 'S22', 'S23', 'S24', 'S25', 'S26', 'S27', 'S28',\n",
+       "       'S29', 'S30', 'S31', 'S32', 'S33', 'S34', 'S35', 'S36', 'S37',\n",
+       "       'S38', 'S39', 'S40', 'S41', 'S42', 'S43', 'S44', 'S45', 'S46',\n",
+       "       'S47', 'S48', 'S49', 'S50', 'S51', 'S52', 'S53', 'S54', 'S55',\n",
+       "       'S56', 'S57', 'S58', 'S59', 'S60', 'S61', 'S62', 'S63', 'S64',\n",
+       "       'S65', 'S66', 'S67', 'S68', 'S69', 'S70', 'S71', 'S72', 'S73',\n",
+       "       'S74', 'S75', 'S76', 'S77', 'S78', 'S79', 'S80', 'S81', 'S82',\n",
+       "       'S83', 'S84', 'S85', 'S86', 'S87', 'S88', 'S89', 'S90', 'S91',\n",
+       "       'S92', 'S93', 'S94', 'S95', 'S96', 'S97', 'S98', 'S99', 'S100',\n",
+       "       'S101', 'S102', 'S103', 'S104', 'S105', 'S106', 'S107', 'S108',\n",
+       "       'S109', 'S110', 'S111', 'S112', 'S113', 'S114', 'S115', 'S116',\n",
+       "       'S117', 'S118', 'S119', 'S120', 'S121', 'S122', 'S123', 'S124',\n",
+       "       'S125', 'S126', 'S127', 'S128', 'S129', 'S130', 'S131', 'S132',\n",
+       "       'S133', 'S134', 'S135', 'S136', 'S137', 'S138', 'S139', 'S140',\n",
+       "       'S141', 'S142', 'S143', 'S144', 'S145', 'S146', 'S147', 'S148',\n",
+       "       'S149', 'S150', 'S151', 'S152', 'S153', 'S154', 'S155', 'S156',\n",
+       "       'S157', 'S158', 'S159', 'S160', 'S161', 'S162', 'S163', 'S164',\n",
+       "       'S165', 'S166', 'S167', 'S168', 'S169', 'S170', 'S171', 'S172',\n",
+       "       'S173', 'S174', 'S175', 'S176', 'S177', 'S178', 'S179', 'S180',\n",
+       "       'S181', 'S182', 'S183', 'S184', 'S185', 'S186', 'S187', 'S188',\n",
+       "       'S189', 'S190', 'S191', 'S192', 'S193', 'S194', 'S195', 'S196',\n",
+       "       'S197', 'S198', 'S199', 'S200', 'S201', 'S202', 'S203', 'S204',\n",
+       "       'S205', 'S206', 'S207', 'S208', 'S209', 'S210', 'S211', 'S212',\n",
+       "       'S213', 'S214', 'S215', 'S216', 'S217', 'S218', 'S219', 'S220',\n",
+       "       'S221', 'S222', 'S223', 'S224', 'S225', 'S226', 'S227', 'S228',\n",
+       "       'S229', 'S230', 'S231', 'S232', 'S233', 'S234', 'S235', 'S236',\n",
+       "       'S237', 'S238', 'S239', 'S240', 'S241', 'S242', 'S243', 'S244',\n",
+       "       'S245', 'S246', 'S247', 'S248', 'S249', 'S250', 'S251', 'S252',\n",
+       "       'S253', 'S254', 'S255', 'S256', 'S257', 'S258', 'S259', 'S260',\n",
+       "       'S261', 'S262', 'S263', 'S264', 'S265', 'S266', 'S267', 'S268',\n",
+       "       'S269', 'S270', 'S271', 'S272', 'S273', 'S274', 'S275', 'S276',\n",
+       "       'S277', 'S278', 'S279', 'S280', 'S281', 'S282', 'S283', 'S284',\n",
+       "       'S285', 'S286', 'S287', 'S288', 'S289', 'S290', 'S291', 'S292',\n",
+       "       'S293', 'S294', 'S295', 'S296', 'S297', 'S298', 'S299', 'S300',\n",
+       "       'S301', 'S302', 'S303', 'S304', 'S305', 'S306', 'S307', 'S308',\n",
+       "       'S309', 'S310', 'S311', 'S312', 'S313', 'S314', 'S315', 'S316',\n",
+       "       'S317', 'S318', 'S319', 'S320', 'S321', 'S322', 'S323', 'S324',\n",
+       "       'S325', 'S326', 'S327', 'S328', 'S329', 'S330', 'S331', 'S332',\n",
+       "       'S333', 'S334', 'S335', 'S336', 'S337', 'S338', 'S339', 'S340',\n",
+       "       'S341', 'S342', 'S343', 'S344', 'S345', 'S346', 'S347', 'S348',\n",
+       "       'S349', 'S350', 'S351', 'S352', 'S353', 'S354', 'S355', 'S356',\n",
+       "       'S357', 'S358', 'S359', 'S360', 'S361', 'S362', 'S363', 'S364',\n",
+       "       'S365', 'S366', 'S367', 'S368', 'S369', 'S370', 'S371', 'S372',\n",
+       "       'S373', 'S374', 'S375', 'S376', 'S377', 'S378', 'S379', 'S380',\n",
+       "       'S381', 'S382', 'S383', 'S384', 'S385', 'S386', 'S387', 'S388',\n",
+       "       'S389', 'S390', 'S391', 'S392', 'S393', 'S394', 'S395', 'S396',\n",
+       "       'S397', 'S398', 'S399', 'S400', 'S401', 'S402', 'S403', 'S404',\n",
+       "       'S405', 'S406', 'S407', 'S408', 'S409', 'S410', 'S411', 'S412',\n",
+       "       'S413', 'S414', 'S415', 'S416', 'S417', 'S418', 'S419', 'S420',\n",
+       "       'S421', 'S422', 'S423', 'S424', 'S425', 'S426', 'S427', 'S428',\n",
+       "       'S429', 'S430', 'S431', 'S432', 'S433', 'S434', 'S435', 'S436',\n",
+       "       'S437', 'S438', 'S439', 'S440', 'S441', 'S442', 'S443', 'S444',\n",
+       "       'S445', 'S446', 'S447', 'S448', 'S449', 'S450', 'S451', 'S452',\n",
+       "       'S453', 'S454', 'S455', 'S456', 'S457', 'S458', 'S459', 'S460',\n",
+       "       'S461', 'S462', 'S463', 'S464', 'S465', 'S466', 'S467', 'S468',\n",
+       "       'S469', 'S470', 'S471', 'S472', 'S473', 'S474', 'S475', 'S476',\n",
+       "       'S477', 'S478', 'S479', 'S480', 'S481', 'S482', 'S483', 'S484',\n",
+       "       'S485', 'S486', 'S487', 'S488', 'S489', 'S490', 'S491', 'S492',\n",
+       "       'S493', 'S494', 'S495', 'S496', 'S497', 'S498', 'S499', 'S500',\n",
+       "       'S501', 'S502', 'S503', 'S504', 'S505', 'S506', 'S507', 'S508',\n",
+       "       'S509', 'S510', 'S511', 'S512', 'S513', 'S514', 'S515', 'S516',\n",
+       "       'S517', 'S518', 'S519', 'S520', 'S521', 'S522', 'S523', 'S524',\n",
+       "       'S525', 'S526', 'S527', 'S528', 'S529', 'S530', 'S531', 'S532',\n",
+       "       'S533', 'S534', 'S535', 'S536', 'S537', 'S538', 'S539', 'S540',\n",
+       "       'S541', 'S542', 'S543', 'S544', 'S545', 'S546', 'S547', 'S548',\n",
+       "       'S549', 'S550', 'S551', 'S552', 'S553', 'S554', 'S555', 'S556',\n",
+       "       'S557', 'S558', 'S559', 'S560', 'S561', 'S562', 'S563', 'S564',\n",
+       "       'S565', 'S566', 'S567', 'S568', 'S569', 'S570', 'S571', 'S572',\n",
+       "       'S573', 'S574', 'S575', 'S576', 'S577', 'S578', 'S579', 'S580',\n",
+       "       'S581', 'S582', 'S583', 'S584', 'S585', 'S586', 'S587', 'S588',\n",
+       "       'S589', 'S590', 'S591', 'S592', 'S593', 'S594', 'S595', 'S596',\n",
+       "       'S597', 'S598', 'S599', 'S600', 'S601', 'S602', 'S603', 'S604',\n",
+       "       'S605', 'S606', 'S607', 'S608', 'S609', 'S610', 'S611', 'S612',\n",
+       "       'S613', 'S614', 'S615', 'S616', 'S617', 'S618', 'S619', 'S620',\n",
+       "       'S621', 'S622', 'S623', 'S624', 'S625', 'S626', 'S627', 'S628',\n",
+       "       'S629', 'S630', 'S631', 'S632', 'S633', 'S634', 'S635', 'S636',\n",
+       "       'S637', 'S638', 'S639', 'S640', 'S641', 'S642', 'S643', 'S644',\n",
+       "       'S645', 'S646', 'S647', 'S648', 'S649', 'S650', 'S651', 'S652',\n",
+       "       'S653', 'S654', 'S655', 'S656', 'S657', 'S658', 'S659', 'S660',\n",
+       "       'S661', 'S662', 'S663', 'S664', 'S665', 'S666', 'S667', 'S668',\n",
+       "       'S669', 'S670', 'S671', 'S672', 'S673', 'S674', 'S675', 'S676',\n",
+       "       'S677', 'S678', 'S679', 'S680', 'S681', 'S682', 'S683', 'S684',\n",
+       "       'S685', 'S686', 'S687', 'S688', 'S689', 'S690', 'S691', 'S692',\n",
+       "       'S693', 'S694', 'S695', 'S696', 'S697', 'S698', 'S699', 'S700',\n",
+       "       'S701', 'S702', 'S703', 'S704', 'S705', 'S706', 'S707', 'S708',\n",
+       "       'S709', 'S710', 'S711', 'S712', 'S713', 'S714', 'S715', 'S716',\n",
+       "       'S717', 'S718', 'S719', 'S720', 'S721', 'S722', 'S723', 'S724',\n",
+       "       'S725', 'S726', 'S727', 'S728', 'S729', 'S730', 'S731', 'S732',\n",
+       "       'S733', 'S734', 'S735', 'S736', 'S737', 'S738', 'S739', 'S740',\n",
+       "       'S741', 'S742', 'S743', 'S744', 'S745', 'S746', 'S747', 'S748',\n",
+       "       'S749', 'S750', 'S751', 'S752', 'S753', 'S754', 'S755', 'S756',\n",
+       "       'S757', 'S758', 'S759', 'S760', 'S761', 'S762', 'S763', 'S764',\n",
+       "       'S765', 'S766', 'S767', 'S768', 'S769', 'S770', 'S771', 'S772',\n",
+       "       'S773', 'S774', 'S775', 'S776', 'S777', 'S778', 'S779', 'S780',\n",
+       "       'S781', 'S782', 'S783', 'S784', 'S785', 'S786', 'S787', 'S788',\n",
+       "       'S789', 'S790', 'S791', 'S792', 'S793', 'S794', 'S795', 'S796',\n",
+       "       'S797', 'S798', 'S799', 'S800', 'S801', 'S802', 'S803', 'S804',\n",
+       "       'S805', 'S806', 'S807', 'S808', 'S809', 'S810', 'S811', 'S812',\n",
+       "       'S813', 'S814', 'S815', 'S816', 'S817', 'S818', 'S819', 'S820',\n",
+       "       'S821', 'S822', 'S823', 'S824', 'S825', 'S826', 'S827', 'S828',\n",
+       "       'S829', 'S830', 'S831', 'S832', 'S833', 'S834', 'S835', 'S836',\n",
+       "       'S837', 'S838', 'S839', 'S840', 'S841', 'S842', 'S843', 'S844',\n",
+       "       'S845', 'S846', 'S847', 'S848', 'S849', 'S850', 'S851', 'S852',\n",
+       "       'S853', 'S854', 'S855', 'S856', 'S857', 'S858', 'S859', 'S860',\n",
+       "       'S861', 'S862', 'S863', 'S864', 'S865', 'S866', 'S867', 'S868',\n",
+       "       'S869', 'S870', 'S871', 'S872', 'S873', 'S874', 'S875', 'S876',\n",
+       "       'S877', 'S878', 'S879', 'S880', 'S881', 'S882', 'S883', 'S884',\n",
+       "       'S885', 'S886', 'S887', 'S888', 'S889', 'S890', 'S891', 'S892',\n",
+       "       'S893', 'S894', 'S895', 'S896', 'S897', 'S898', 'S899', 'S900',\n",
+       "       'S901', 'S902', 'S903', 'S904', 'S905', 'S906', 'S907', 'S908',\n",
+       "       'S909', 'S910', 'S911', 'S912', 'S913', 'S914', 'S915', 'S916',\n",
+       "       'S917', 'S918', 'S919', 'S920', 'S921', 'S922', 'S923', 'S924',\n",
+       "       'S925', 'S926', 'S927', 'S928', 'S929', 'S930', 'S931', 'S932',\n",
+       "       'S933', 'S934', 'S935', 'S936', 'S937', 'S938', 'S939', 'S940',\n",
+       "       'S941', 'S942', 'S943', 'S944', 'S945', 'S946', 'S947', 'S948',\n",
+       "       'S949', 'S950', 'S951', 'S952', 'S953', 'S954', 'S955', 'S956',\n",
+       "       'S957', 'S958', 'S959', 'S960', 'S961', 'S962', 'S963', 'S964',\n",
+       "       'S965', 'S966', 'S967', 'S968', 'S969', 'S970', 'S971', 'S972',\n",
+       "       'S973', 'S974', 'S975', 'S976', 'S977', 'S978', 'S979', 'S980',\n",
+       "       'S981', 'S982', 'S983', 'S984', 'S985', 'S986', 'S987', 'S988',\n",
+       "       'S989', 'S990', 'S991', 'S992', 'S993', 'S994', 'S995', 'S996',\n",
+       "       'S997', 'S998', 'S999'], dtype='<U4')
      • call/genotype
        (variants, samples, ploidy)
        int64
        1 0 1 0 0 0 1 0 ... 1 0 1 1 0 0 1 1
        array([[[1, 0],\n",
+       "        [1, 0],\n",
+       "        [0, 0],\n",
+       "        ...,\n",
+       "        [1, 0],\n",
+       "        [1, 0],\n",
+       "        [1, 1]],\n",
+       "\n",
+       "       [[1, 0],\n",
+       "        [1, 1],\n",
+       "        [1, 1],\n",
+       "        ...,\n",
+       "        [1, 1],\n",
+       "        [0, 0],\n",
+       "        [1, 1]],\n",
+       "\n",
+       "       [[1, 0],\n",
+       "        [1, 1],\n",
+       "        [1, 0],\n",
+       "        ...,\n",
+       "        [1, 0],\n",
+       "        [0, 0],\n",
+       "        [1, 0]],\n",
+       "\n",
+       "       ...,\n",
+       "\n",
+       "       [[1, 1],\n",
+       "        [1, 0],\n",
+       "        [1, 0],\n",
+       "        ...,\n",
+       "        [0, 0],\n",
+       "        [0, 0],\n",
+       "        [1, 0]],\n",
+       "\n",
+       "       [[1, 0],\n",
+       "        [1, 0],\n",
+       "        [1, 1],\n",
+       "        ...,\n",
+       "        [0, 0],\n",
+       "        [1, 0],\n",
+       "        [1, 0]],\n",
+       "\n",
+       "       [[0, 0],\n",
+       "        [0, 0],\n",
+       "        [0, 0],\n",
+       "        ...,\n",
+       "        [1, 1],\n",
+       "        [0, 0],\n",
+       "        [1, 1]]])
      • call/genotype_mask
        (variants, samples, ploidy)
        bool
        False False False ... False False
        array([[[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]],\n",
+       "\n",
+       "       [[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]],\n",
+       "\n",
+       "       [[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]],\n",
+       "\n",
+       "       ...,\n",
+       "\n",
+       "       [[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]],\n",
+       "\n",
+       "       [[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]],\n",
+       "\n",
+       "       [[False, False],\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        ...,\n",
+       "        [False, False],\n",
+       "        [False, False],\n",
+       "        [False, False]]])
      • variant/hwe_p_value
        (variants)
        float64
        dask.array<chunksize=(50,), meta=np.ndarray>
        \n",
+       "\n",
+       "\n",
+       "\n",
+       "\n",
+       "
        \n",
+       "\n",
+       "  \n",
+       "    \n",
+       "  \n",
+       "  \n",
+       "    \n",
+       "    \n",
+       "    \n",
+       "    \n",
+       "  \n",
+       "
          Array  Chunk 
         Bytes  400 B   400 B 
         Shape  (50,)   (50,) 
         Count  4 Tasks  1 Chunks 
         Type  float64  numpy.ndarray 
        \n",
+       "        		\n",
+       "\n",
+       "\n",
+       "  \n",
+       "  \n",
+       "  \n",
+       "\n",
+       "  \n",
+       "  \n",
+       "  \n",
+       "\n",
+       "  \n",
+       "  \n",
+       "\n",
+       "  \n",
+       "  50\n",
+       "  1\n",
+       "\n",
+       "
    • contigs :
      [0]
    " + ], + "text/plain": [ + "\n", + "Dimensions: (alleles: 2, ploidy: 2, samples: 1000, variants: 50)\n", + "Dimensions without coordinates: alleles, ploidy, samples, variants\n", + "Data variables:\n", + " variant/contig (variants) int64 0 0 0 0 0 0 0 0 0 ... 0 0 0 0 0 0 0 0\n", + " variant/position (variants) int64 0 1 2 3 4 5 6 ... 43 44 45 46 47 48 49\n", + " variant/alleles (variants, alleles) |S1 b'T' b'C' b'T' ... b'T' b'C'\n", + " sample/id (samples) \n", + "Attributes:\n", + " contigs: [0]" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from typing import Optional\n", + "from sgkit.api import create_genotype_call_dataset\n", + "from sgkit.stats.hwe import hardy_weinberg_test\n", + "from sgkit.tests.test_hwe import to_genotype_call_dataset, simulate_genotype_calls\n", + "gt_dist = [0.25, 0.5, 0.25]\n", + "call_genotype = simulate_genotype_calls(50, 1000, p=gt_dist)\n", + "ds = to_genotype_call_dataset(call_genotype)\n", + "ds = ds.merge(hardy_weinberg_test(ds))\n", + "ds" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0.84965934, 0.16354278, 0.07601487, 0.5275052 , 0.6128462 ,\n", + " 1. , 0.70306802, 0.94962855, 0.41039482, 0.56905156,\n", + " 0.01141049, 0.04312044, 0.16446281, 0.34292471, 0.65827469,\n", + " 0.65828459, 0.7046632 , 0.41143553, 1. , 0.56913128,\n", + " 0.80049101, 0.75171363, 0.18377625, 0.34334841, 0.16443965,\n", + " 0.94961921, 0.41048285, 0.11418474, 0.89944443, 0.48613015,\n", + " 0.84946805, 0.31078525, 0.03701711, 0.61322632, 0.48702148,\n", + " 0.94963516, 0.89943975, 0.48450522, 0.56908112, 0.28101448,\n", + " 0.00194671, 0.80023881, 0.44766596, 0.23003134, 0.44829854,\n", + " 0.44830607, 0.84946805, 0.4094545 , 0.37569028, 0.22925408])" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p = ds['variant/hwe_p_value'].values\n", + "p" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
    \n", + "\n", + "\n", + "Show/Hide data repr\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Show/Hide attributes\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
    xarray.Dataset
      • ploidy: 3
        • x
          (ploidy)
          float64
          0.0 0.0 0.0
          array([0., 0., 0.])
      " + ], + "text/plain": [ + "\n", + "Dimensions: (ploidy: 3)\n", + "Dimensions without coordinates: ploidy\n", + "Data variables:\n", + " x (ploidy) float64 0.0 0.0 0.0" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import xarray as xr\n", + "xr.Dataset({'x': ('ploidy', np.zeros(3))})" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(0.0019467095399049773, 1.0)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "p.min(), p.max()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/validation/gwas/method/hwe/invoke.yaml b/validation/gwas/method/hwe/invoke.yaml new file mode 100644 index 000000000..ff787a8ca --- /dev/null +++ b/validation/gwas/method/hwe/invoke.yaml @@ -0,0 +1,2 @@ +tasks: + auto_dash_names: false \ No newline at end of file diff --git a/validation/gwas/method/hwe/logging.ini b/validation/gwas/method/hwe/logging.ini new file mode 100644 index 000000000..604feb7d0 --- /dev/null +++ b/validation/gwas/method/hwe/logging.ini @@ -0,0 +1,21 @@ +[loggers] +keys=root + +[handlers] +keys=console + +[formatters] +keys=console_formatter + +[logger_root] +level=INFO +handlers=console + +[handler_console] +level=INFO +class=StreamHandler +formatter=console_formatter +args=(sys.stdout,) + +[formatter_console_formatter] +format=%(asctime)s|%(levelname)s|%(name)s.%(funcName)s:%(lineno)d| %(message)s \ No newline at end of file diff --git a/validation/gwas/method/hwe/tasks.py b/validation/gwas/method/hwe/tasks.py new file mode 100644 index 000000000..f4da9666b --- /dev/null +++ b/validation/gwas/method/hwe/tasks.py @@ -0,0 +1,80 @@ +import ctypes +import glob +import logging +import logging.config +import os +import shutil +from pathlib import Path + +import numpy as np +import pandas as pd +from invoke import task + +logging.config.fileConfig("logging.ini") +logger = logging.getLogger(__name__) + +DEFAULT_SIM_DATADIR = os.getenv("SIM_DATADIR", "data") +DEFAULT_TEST_DATADIR = os.getenv("TEST_DATADIR", "../../../../sgkit/tests/test_hwe") + + +@task +def compile(ctx): + """Build reference implementation C library""" + logger.info("Building reference C library") + ctx.run("make") + logger.info("Build complete") + + +def get_genotype_counts(): + """Generate genotype counts for testing.""" + rs = np.random.RandomState(0) + n, s = 10_000, 50 + n_het = np.expand_dims(np.arange(n, step=s) + 1, -1) + frac = rs.uniform(0.3, 0.7, size=(n // s, 2)) + n_hom = frac * n_het + n_hom = n_hom.astype(int) + return pd.DataFrame( + np.concatenate((n_het, n_hom), axis=1), columns=["n_het", "n_hom_1", "n_hom_2"] + ) + + +@task +def simulate(ctx, sim_datadir=DEFAULT_SIM_DATADIR): + """Create inputs and outputs for unit tests.""" + logger.info("Generating unit test data") + libc = ctypes.CDLL("./libchwe.so") + chwep = libc.hwep + chwep.restype = ctypes.c_double + df = get_genotype_counts() + df["p"] = df.apply( + lambda r: chwep(int(r["n_het"]), int(r["n_hom_1"]), int(r["n_hom_2"])), axis=1 + ) + output_dir = Path(sim_datadir) + if not output_dir.exists(): + output_dir.mkdir(parents=True, exist_ok=True) + path = output_dir / "sim_01.csv" + df.to_csv(path, index=False) + logger.info(f"Unit test data written to {path}") + + +@task +def export( + ctx, + sim_datadir=DEFAULT_SIM_DATADIR, + test_datadir=DEFAULT_TEST_DATADIR, + clear=True, + runs=None, +): + sim_datadir = Path(sim_datadir) + test_datadir = Path(test_datadir).resolve() + logger.info(f"Exporting test data to {test_datadir}") + if clear and test_datadir.exists(): + logger.info(f"Clearing test datadir at {test_datadir}") + shutil.rmtree(test_datadir) + test_datadir.mkdir(exist_ok=True) + for f in glob.glob(str(sim_datadir / "*.csv")): + src = f + dst = test_datadir / Path(f).name + logger.info(f"Copying {src} to {dst}") + shutil.copy(src, dst) + logger.info("Export complete")