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