830 - Haplotyper::SampleChromosomes()

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);
  }