From Genome Analysis Wiki
Jump to navigationJump to search
1,257 bytes added
, 15:01, 25 September 2013
<source lang="cpp">
void MarkovModel::CountErrorsAndRecombinants(char * observed, char ** haplotypes, float ** freq)
{
float * swap;
float * vector = new float [states];
float * extra = new float [states];
// Initialize likelihoods at first position
for (int i = 0; i < states; i++)
vector[i] = 1.;
// Scan along chromosome
for (int i = markers - 1; i > 0; i--)
{
for (int j = 0; j < states; j++)
extra[j] = vector[j] * matrix[i][j];
if (observed[i])
{
errorCount[i] += CountErrors(extra, haplotypes, i, observed[i], E[i], freq[observed[i]][i]);
Condition(vector, haplotypes, i, observed[i], E[i], freq[observed[i]][i]);
}
else
errorCount[i] += E[i];
recombinantCount[i-1] += CountRecombinants(vector, matrix[i - 1], R[i-1]);
Transpose(vector, extra, R[i - 1]);
swap = vector; vector = extra; extra = swap;
}
if (observed[0])
{
Condition(vector, haplotypes, 0, observed[0], E[0], freq[observed[0]][0]);
errorCount[0] += CountErrors(vector, haplotypes, 0, observed[0], E[0], freq[observed[0]][0]);
}
else
errorCount[0] += E[0];
delete [] vector;
delete [] extra;
}
</source>