830 - Haplotyper::RandomSetup()

void Haplotyper::RandomSetup(Random * rand)
   {
   // accesses two arrays... 
   //    genotypes[i][j] includes genotypes for each individual
   //        i ranges from 0 .. individuals - 1
   //        j ranges from 0 .. markers - 1
   //        coded as 0 (missing), 1, 2, 3
   //
   //    haplotypes[i][j] includes haplotypes for each individual
   //        i ranges from 0 .. 2 * individuals
   //        coded as 0 or 1, no missing data

   for (int j = 0; j < markers; j++)
      {
      int alleles = 0, mac = 0;

      for (int i = 0; i < individuals; i++)
         if (genotypes[i][j] != GENOTYPE_MISSING)
            alleles += 2, mac += genotypes[i][j] - 1;

      for (int i = individuals; i < individuals; i++)
         {
         mac += haplotypes[i * 2][j] == 1;
         mac += haplotypes[i * 2 + 1][j] == 1;
         }

      if (alleles == 0)
         {
         for (int i = 0; i < individuals; i++)
            haplotypes[i * 2][j] = haplotypes[i * 2 + 1][j] = 0;
         continue;
         }

      float freq = mac / (float) alleles;

      for (int i = 0; i < individuals; i++)
         switch (genotypes[i][j])
            {
            case 0: /* MISSING GENOTYPE */
              haplotypes[i * 2][j] = rand->Next() < freq;
              haplotypes[i * 2 + 1][j] = rand->Next() < freq;
              break;
            case 1: /* HOMOZYGOUS FOR ALLELE 1 */
              haplotypes[i * 2][j] = 0;
              haplotypes[i * 2 + 1][j] = 0;
              break;
            case 2: /* HETEROZYGOUS */
              {
              int bit = rand->Binary();

              haplotypes[i * 2][j] = bit;
              haplotypes[i * 2 + 1][j] = bit ^ 1;
              }
              break;
            case 3: /* HOMOZYGOUS FOR ALLELE 2 */
              haplotypes[i * 2][j] = 1;
              haplotypes[i * 2 + 1][j] = 1;
              break;
            }
      }
   }