Difference between revisions of "830 - Haplotyper::SampleChromosomes()"
From Genome Analysis Wiki
Jump to navigationJump to search (Created page with " void Haplotyper::SampleChromosomes(Random * rand) { RewindMemoryPool(); RetrieveMemoryBlock(markers - 1); int first = 0, second = 0; SampleHaplotypes(leftMat...") |
|||
Line 1: | Line 1: | ||
− | + | <source lang="cpp"> | |
void Haplotyper::SampleChromosomes(Random * rand) | void Haplotyper::SampleChromosomes(Random * rand) | ||
Line 44: | Line 44: | ||
RecordHaplotypes(0, first, second); | RecordHaplotypes(0, first, second); | ||
} | } | ||
+ | |||
+ | </source> |
Latest revision as of 14:11, 2 October 2013
void Haplotyper::SampleChromosomes(Random * rand)
{
RewindMemoryPool();
RetrieveMemoryBlock(markers - 1);
int first = 0, second = 0;
SampleHaplotypes(leftMatrices[markers - 1], first, second, rand);
float left[2][2], right[2][2];
for (int j = markers - 2; j >= 0; j--)
{
RecordHaplotype(j + 1, first, second);
RetrieveMemoryBlock(j);
float sum = SummarizeOptions(leftMatrices[j], left, right, thetas[j], first, second);
float choice = rand->Uniform(0, sum);
// The most likely outcome is that no changes occur ...
choice -= left[0][0] * right[0][0];
if (choice <= 0.0) continue;
// Otherwise, check if the second haplotype changed ...
choice -= left[0][1] * right[0][1];
if (choice <= 0.0)
{
SampleOneHaplotype(leftMatrices[j], first, second, -choice / right[0][1]);
continue;
}
// Otherwise, check if the first haplotype changed ...
choice -= left[1][0] * right[1][0];
if (choice <= 0.0)
{
SampleOneHaplotype(leftMatrices[j], second, first, -choice / right[1][0]);
continue;
}
// And sometimes, both haplotypes recombine and we start fresh
SampleHaplotypes(leftMatrices[j], first, second, rand);
}
RecordHaplotypes(0, first, second);
}