-static double mc_gtfreq_iter(double g[3], const bcf_p1aux_t *ma, int beg, int end)
-{
- double err, gg[3];
- int i;
- gg[0] = gg[1] = gg[2] = 0.;
- for (i = beg; i < end; ++i) {
- double *pdg, sum, tmp[3];
- pdg = ma->pdg + i * 3;
- tmp[0] = pdg[0] * g[0]; tmp[1] = pdg[1] * g[1]; tmp[2] = pdg[2] * g[2];
- sum = (tmp[0] + tmp[1] + tmp[2]) * (end - beg);
- gg[0] += tmp[0] / sum; gg[1] += tmp[1] / sum; gg[2] += tmp[2] / sum;
- }
- err = fabs(gg[0] - g[0]) > fabs(gg[1] - g[1])? fabs(gg[0] - g[0]) : fabs(gg[1] - g[1]);
- err = err > fabs(gg[2] - g[2])? err : fabs(gg[2] - g[2]);
- g[0] = gg[0]; g[1] = gg[1]; g[2] = gg[2];
- return err;
+ double norm=lk[0];
+ for (i=1; i<n; i++) if (lk[i]<norm) norm=lk[i];
+ #if DBG
+ for (i=0,j=0; i<b->n_alleles; i++)
+ {
+ int k; for (k=0; k<=i; k++) printf("%.0f\t", lk[j++]);
+ printf("\n");
+ }
+ #endif
+ for (i=0; i<n; i++) lk[i] = pow(10,-0.1*(lk[i]-norm));
+
+ // Find out the most likely alleles. In contrast to the original version,
+ // ALT alleles may not be printed when they are more likely than REF but
+ // significantly less likely than the most likely ALT. The only criterion
+ // is the LK ratio now. To obtain behaviour similar to the original one,
+ // use the pref variant below.
+ double pmax=0; //,pref=0;
+ n = ma->is_indel ? b->n_alleles : b->n_alleles-1;
+ for (i=0; i<n; i++)
+ {
+ double pr=0;
+ int k=i*(i+1)/2;
+ for (j=0; j<=i; j++) { pr+=lk[k]; k++; }
+ for (j=i+1; j<b->n_alleles; j++) { k=j*(j+1)/2+i; pr+=lk[k]; }
+ #if DBG
+ printf("%d\t%e\n", i,pr);
+ #endif
+ if (pmax<pr) pmax=pr;
+ // if (i==0) pref=pr;
+ // if (pr<pref && pr/pmax < 1e-4) break;
+ if (pr/pmax < 1e-4) break; // Assuming the alleles are sorted by the lk
+ }
+ return i-1;