- // generate esum and fsum
- memset(&aux, 0, sizeof(auxaux_t));
- for (j = n - 1, r->sum_Q2 = 0; j >= 0; --j) { // calculate esum and fsum
- uint32_t info = bca->info[j];
- int tmp;
- if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
- k = info>>16&7;
- if (info>>24 > 0) {
- aux.esum[k&3] += bca->fk[aux.w[k]] * (info>>24);
- aux.fsum[k&3] += bca->fk[aux.w[k]];
- if (aux.w[k] + 1 < CALL_MAX) ++aux.w[k];
- ++aux.c[k&3];
- }
- tmp = (int)(info&0xff) < bca->capQ? (int)(info&0xff) : bca->capQ;
- r->sum_Q2 += tmp * tmp;
- }
- memcpy(r->esum, aux.esum, 5 * sizeof(float));
- // rescale ->c[]
- for (j = c = 0; j != 4; ++j) c += aux.c[j];
- if (c > 255) {
- for (j = 0; j != 4; ++j) aux.c[j] = (int)(254.0 * aux.c[j] / c + 0.499);
- for (j = c = 0; j != 4; ++j) c += aux.c[j];
- }
- // generate likelihood
- for (j = 0; j != 5; ++j) {
- float tmp;
- // homozygous
- for (k = 0, tmp = 0.0; k != 5; ++k)
- if (j != k) tmp += aux.esum[k];
- p[j*5+j] = tmp; // anything that is not j
- // heterozygous
- for (k = j + 1; k < 5; ++k) {
- for (i = 0, tmp = 0.0; i != 5; ++i)
- if (i != j && i != k) tmp += aux.esum[i];
- p[j*5+k] = p[k*5+j] = 3.01 * (aux.c[j] + aux.c[k]) + tmp;
- }
- }