1,555
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to search830 - Haplotyper::SampleHaplotypes() (view source)
Revision as of 14:05, 2 October 2013
, 14:05, 2 October 2013Created page with "<source lang="cpp"> void Haplotyper::SampleChromosomes(Random * rand) { RewindMemoryPool(); RetrieveMemoryBlock(markers - 1); int first = 0, second = 0; Sampl..."
<source lang="cpp">
void Haplotyper::SampleChromosomes(Random * rand)
{
RewindMemoryPool();
RetrieveMemoryBlock(markers - 1);
int first = 0, second = 0;
SampleHaplotypes(leftMatrices[markers - 1], first, second);
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(probability, first, second);
}
RecordHaplotypes(0, first, second);
}
</source>
void Haplotyper::SampleChromosomes(Random * rand)
{
RewindMemoryPool();
RetrieveMemoryBlock(markers - 1);
int first = 0, second = 0;
SampleHaplotypes(leftMatrices[markers - 1], first, second);
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(probability, first, second);
}
RecordHaplotypes(0, first, second);
}
</source>
