- for (i = 0; i < b->n_alleles; ++i) p[i] = p[i]<<4 | i;
- for (i = 1; i < b->n_alleles; ++i) // insertion sort
- for (j = i; j > 0 && p[j] < p[j-1]; --j)
- tmp = p[j], p[j] = p[j-1], p[j-1] = tmp;
- for (i = b->n_alleles - 1; i >= 0; --i)
- if ((p[i]&0xf) == 0) break;
- return i;
-}
-// f0 is the reference allele frequency
-static double mc_freq_iter(double f0, const bcf_p1aux_t *ma)
-{
- double f, f3[3];
- int i;
- f3[0] = (1.-f0)*(1.-f0); f3[1] = 2.*f0*(1.-f0); f3[2] = f0*f0;
- for (i = 0, f = 0.; i < ma->n; ++i) {
- double *pdg;
- pdg = ma->pdg + i * 3;
- f += (pdg[1] * f3[1] + 2. * pdg[2] * f3[2])
- / (pdg[0] * f3[0] + pdg[1] * f3[1] + pdg[2] * f3[2]);
- }
- f /= ma->n * 2.;
- return f;