text: //
  // freebayes
  //
  // A bayesian genetic variant detector.
  //
  
  // standard includes
  //#include <cstdio>
  #include <iostream>
  #include <fstream>
  #include <string>
  #include <vector>
  #include <map>
  #include <iterator>
  #include <algorithm>
  #include <cmath>
  #include <time.h>
  #include <float.h>
  #include <stdlib.h>
  
  // private libraries
  #ifdef HAVE_BAMTOOLS
  #include "api/BamReader.h"
  #endif
  #include "TryCatch.h"
  #include "Parameters.h"
  #include "Allele.h"
  #include "Sample.h"
  #include "AlleleParser.h"
  #include "Utility.h"
  #include "SegfaultHandler.h"
  
  #include "multichoose.h"
  #include "multipermute.h"
  
  #include "Genotype.h"
  #include "DataLikelihood.h"
  #include "Marginals.h"
  #include "ResultData.h"
  
  #include "Bias.h"
  #include "Contamination.h"
  #include "NonCall.h"
  #include "Logging.h"
  
  using namespace std;
  
  // todo
  // generalize the main function to take the parameters as input
  // so that we can invoke the entire algorithm on different regions
  // when requested to run in parallel
  // take the targets (or whole genome) and make small jobs
  // run the main function for each region in an omp parallel for loop
  // only do this if the --parallel flag is set > 1
  
  // freebayes main
  int main (int argc, char *argv[]) {
  
      // install segfault handler
      signal(SIGSEGV, segfaultHandler);
  
      AlleleParser* parser = new AlleleParser(argc, argv);
      Parameters& parameters = parser->parameters;
      list<Allele*> alleles;
  
      Samples samples;
      NonCalls nonCalls;
  
      ostream& out = *(parser->output);
  
      Bias observationBias;
      if (!parameters.alleleObservationBiasFile.empty()) {
          observationBias.open(parameters.alleleObservationBiasFile);
      }
  
      Contamination contaminationEstimates(0.5+parameters.probContamination, parameters.probContamination);
      if (!parameters.contaminationEstimateFile.empty()) {
          contaminationEstimates.open(parameters.contaminationEstimateFile);
      }
  
      // this can be uncommented to force operation on a specific set of genotypes
      vector<Allele> allGenotypeAlleles;
      allGenotypeAlleles.push_back(genotypeAllele(ALLELE_GENOTYPE, "A", 1));
      allGenotypeAlleles.push_back(genotypeAllele(ALLELE_GENOTYPE, "T", 1));
      allGenotypeAlleles.push_back(genotypeAllele(ALLELE_GENOTYPE, "G", 1));
      allGenotypeAlleles.push_back(genotypeAllele(ALLELE_GENOTYPE, "C", 1));
  
      int allowedAlleleTypes = ALLELE_REFERENCE;
      if (parameters.allowSNPs) {
          allowedAlleleTypes |= ALLELE_SNP;
      }
      if (parameters.allowIndels) {
          allowedAlleleTypes |= ALLELE_INSERTION;
          allowedAlleleTypes |= ALLELE_DELETION;
      }
      if (parameters.allowMNPs) {
          allowedAlleleTypes |= ALLELE_MNP;
      }
      if (parameters.allowComplex) {
          allowedAlleleTypes |= ALLELE_COMPLEX;
      }
  
      // output VCF header
      if (parameters.output == "vcf" && !parameters.noOutput) {
          out << parser->variantCallFile.header << endl;
      }
  
      if (0 < parameters.limitCoverage) {
          srand(13);
      }
  
      Allele nullAllele = genotypeAllele(ALLELE_NULL, "N", 1, "1N");
  
      unsigned long total_sites = 0;
      unsigned long processed_sites = 0;
  
      while (parser->getNextAlleles(samples, allowedAlleleTypes)) {
  
          ++total_sites;
  
          DEBUG2("at start of main loop");
  
          // did we switch chromosomes or exceed our gVCF chunk size, or do we not want to use chunks?
          // if so, we may need to output a gVCF record
          Results results;
          if (parameters.gVCFout
                 &&  !(nonCalls.empty())
                 &&  (  (parameters.gVCFNoChunk)
                     || (nonCalls.begin()->first != parser->currentSequenceName)
                     || (parameters.gVCFchunk
                         && nonCalls.lastPos().second - nonCalls.firstPos().second >= parameters.gVCFchunk
                        )
                    )
              ){
              vcflib::Variant var(parser->variantCallFile);
              out << results.gvcf(var, nonCalls, parser) << endl;
              nonCalls.clear();
          }
  
          // don't process non-ATGCN's in the reference
          string cb = parser->currentReferenceBaseString();
          if (cb != "A" && cb != "T" && cb != "C" && cb != "G" && cb != "N") {
              DEBUG2("current reference base is not in { A T G C N }");
              continue;
          }
  
          int coverage = countAlleles(samples);
  
          DEBUG("position: " << parser->currentSequenceName << ":" << (long unsigned int) parser->currentPosition + 1 << " coverage: " << coverage);
  
          bool skip = false;
          if (!parser->hasInputVariantAllelesAtCurrentPosition()) {
              // skips 0-coverage regions
              if (coverage == 0) {
                  DEBUG("no alleles left at this site after filtering");
                  skip = true;
              } else if (coverage < parameters.minCoverage) {
                  DEBUG("post-filtering coverage of " << coverage << " is less than --min-coverage of " << parameters.minCoverage);
                  skip = true;
              } else if (parameters.onlyUseInputAlleles) {
                  DEBUG("no input alleles, but using only input alleles for analysis, skipping position");
                  skip = true;
              } else if (0 < parameters.limitCoverage) {
                  // go through each sample
                  for (Samples::iterator s = samples.begin(); s != samples.end(); ++s) {
                      string sampleName = s->first;
                      Sample& sample = s->second;
                      // get the coverage for this sample
                      int sampleCoverage = 0;
                      for (Sample::iterator sg = sample.begin(); sg != sample.end(); ++sg) {
                          sampleCoverage += sg->second.size();
                      }
                      if (sampleCoverage <= parameters.limitCoverage) {
                          continue;
                      }
  
                      DEBUG("coverage " << sampleCoverage << " for sample " << sampleName << " was > " << parameters.limitCoverage << ", so we will remove " << (sampleCoverage - parameters.limitCoverage) << " genotypes");
                      vector<string> genotypesToErase;
                      do {
                          double probRemove = (sampleCoverage - parameters.limitCoverage) / (double)sampleCoverage;
                          vector<string> genotypesToErase;
                          // iterate through the genotypes
                          for (Sample::iterator sg = sample.begin(); sg != sample.end(); ++sg) {
                              vector<Allele*> allelesToKeep;
                              // iterate through each allele
                              for (int alleleIndex = 0; alleleIndex < sg->second.size(); alleleIndex++) {
                                  // only if we have more alleles to remove
                                  if (parameters.limitCoverage < sampleCoverage) {
                                      double r = rand() / (double)RAND_MAX;
                                      if (r < probRemove) { // skip over this allele
                                          sampleCoverage--;
                                          continue;
                                      }
                                  }
                                  // keep it
                                  allelesToKeep.push_back(sg->second[alleleIndex]);
                              }
                              // re-assign the alleles to this genotype
                              if (allelesToKeep.size() < sg->second.size()) {
                                  sg->second.assign(allelesToKeep.begin(), allelesToKeep.end());
                              }
                              // if no more alleles for this genotype, remove it later
                              if (sg->second.empty()) {
                                  genotypesToErase.push_back(sg->first);
                              }
                          }
                          // remove empty genotypes
                          for (vector<string>::iterator gt = genotypesToErase.begin(); gt != genotypesToErase.end(); ++gt) {
                              sample.erase(*gt);
                          }
                      } while (parameters.limitCoverage < sampleCoverage);
                      sampleCoverage = 0;
                      for (Sample::iterator sg = sample.begin(); sg != sample.end(); ++sg) {
                          sampleCoverage += sg->second.size();
                      }
                      DEBUG("coverage for sample " << sampleName << " is now " << sampleCoverage);
                  }
                  // update coverage
                  coverage = countAlleles(samples);
              }
  
              DEBUG2("coverage " << parser->currentSequenceName << ":" << parser->currentPosition << " == " << coverage);
  
              // establish a set of possible alternate alleles to evaluate at this location
  
              if (!parameters.reportMonomorphic
                  && !sufficientAlternateObservations(samples, parameters.minAltCount, parameters.minAltFraction)) {
                  DEBUG("insufficient alternate observations");
                  skip = true;
              }
              if (parameters.reportMonomorphic) {
                  DEBUG("calling at site even though there are no alternate observations");
              }
          } else {
              /*
              cerr << "has input variants at " << parser->currentSequenceName << ":" << parser->currentPosition << endl;
              vector<Allele>& inputs = parser->inputVariantAlleles[parser->currentSequenceName][parser->currentPosition];
              for (vector<Allele>::iterator a = inputs.begin(); a != inputs.end(); ++a) {
                  cerr << *a << endl;
              }
              */
          }
  
          if (skip) {
              // record data for gVCF
              if (parameters.gVCFout) {
                  nonCalls.record(parser->currentSequenceName, parser->currentPosition, samples);
              }
              // and step ahead
              continue;
          }
  
          // to ensure proper ordering of output stream
          vector<string> sampleListPlusRef;
  
          for (vector<string>::iterator s = parser->sampleList.begin(); s != parser->sampleList.end(); ++s) {
              sampleListPlusRef.push_back(*s);
          }
          if (parameters.useRefAllele) {
              sampleListPlusRef.push_back(parser->currentSequenceName);
          }
  
          // establish genotype alleles using input filters
          map<string, vector<Allele*> > alleleGroups;
          groupAlleles(samples, alleleGroups);
          DEBUG2("grouped alleles by equivalence");
  
          vector<Allele> genotypeAlleles = parser->genotypeAlleles(alleleGroups, samples, parameters.onlyUseInputAlleles);
  
          // always include the reference allele as a possible genotype, even when we don't include it by default
          if (!parameters.useRefAllele) {
              vector<Allele> refAlleleVector;
              refAlleleVector.push_back(genotypeAllele(ALLELE_REFERENCE, string(1, parser->currentReferenceBase), 1, "1M"));
              genotypeAlleles = alleleUnion(genotypeAlleles, refAlleleVector);
          }
  
          map<string, vector<Allele*> > partialObservationGroups;
          map<Allele*, set<Allele*> > partialObservationSupport;
  
          // build haplotype alleles matching the current longest allele (often will do nothing)
          // this will adjust genotypeAlleles if changes are made
          DEBUG("building haplotype alleles, currently there are " << genotypeAlleles.size() << " genotype alleles");
          DEBUG(genotypeAlleles);
          parser->buildHaplotypeAlleles(genotypeAlleles,
                                        samples,
                                        alleleGroups,
                                        partialObservationGroups,
                                        partialObservationSupport,
                                        allowedAlleleTypes);
          DEBUG("built haplotype alleles, now there are " << genotypeAlleles.size() << " genotype alleles");
          DEBUG(genotypeAlleles);
  
          string referenceBase = parser->currentReferenceHaplotype();
  
  
          // for debugging
          /*
          for (Samples::iterator s = samples.begin(); s != samples.end(); ++s) {
              string sampleName = s->first;
              Sample& sample = s->second;
              cerr << sampleName << ": " << sample << endl;
          }
          */
  
          // re-calculate coverage, as this could change now that we've built haplotype alleles
          coverage = countAlleles(samples);
  
          // estimate theta using the haplotype length
          long double theta = parameters.TH * parser->lastHaplotypeLength;
  
          // if we have only one viable allele, we don't have evidence for variation at this site
          if (!parser->hasInputVariantAllelesAtCurrentPosition() && !parameters.reportMonomorphic && genotypeAlleles.size() <= 1 && genotypeAlleles.front().isReference()) {
              DEBUG("no alternate genotype alleles passed filters at " << parser->currentSequenceName << ":" << parser->currentPosition);
              continue;
          }
          DEBUG("genotype alleles: " << genotypeAlleles);
  
          // add the null genotype
          bool usingNull = false;
          if (parameters.excludeUnobservedGenotypes && genotypeAlleles.size() > 2) {
              genotypeAlleles.push_back(nullAllele);
              usingNull = true;
          }
  
          if (parameters.noOutput) {
              continue;
          }
  
          ++processed_sites;
  
          // generate possible genotypes
  
          // for each possible ploidy in the dataset, generate all possible genotypes
          vector<int> ploidies = parser->currentPloidies(samples);
          map<int, vector<Genotype> > genotypesByPloidy = getGenotypesByPloidy(ploidies, genotypeAlleles);
          int numCopiesOfLocus = parser->copiesOfLocus(samples);
  
  
          DEBUG2("generated all possible genotypes:");
          if (parameters.debug2) {
              for (map<int, vector<Genotype> >::iterator s = genotypesByPloidy.begin(); s != genotypesByPloidy.end(); ++s) {
                  vector<Genotype>& genotypes = s->second;
                  for (vector<Genotype>::iterator g = genotypes.begin(); g != genotypes.end(); ++g) {
                      DEBUG2(*g);
                  }
              }
          }
  
          // get estimated allele frequencies using sum of estimated qualities
          map<string, double> estimatedAlleleFrequencies = samples.estimatedAlleleFrequencies();
          double estimatedMaxAlleleFrequency = 0;
          double estimatedMaxAlleleCount = 0;
          double estimatedMajorFrequency = estimatedAlleleFrequencies[referenceBase];
          if (estimatedMajorFrequency < 0.5) estimatedMajorFrequency = 1-estimatedMajorFrequency;
          double estimatedMinorFrequency = 1-estimatedMajorFrequency;
          //cerr << "num copies of locus " << numCopiesOfLocus << endl;
          int estimatedMinorAllelesAtLocus = max(1, (int) ceil((double) numCopiesOfLocus * estimatedMinorFrequency));
          //cerr << "estimated minor frequency " << estimatedMinorFrequency << endl;
          //cerr << "estimated minor count " << estimatedMinorAllelesAtLocus << endl;
  
  
          map<string, vector<vector<SampleDataLikelihood> > > sampleDataLikelihoodsByPopulation;
          map<string, vector<vector<SampleDataLikelihood> > > variantSampleDataLikelihoodsByPopulation;
          map<string, vector<vector<SampleDataLikelihood> > > invariantSampleDataLikelihoodsByPopulation;
  
          map<string, int> inputAlleleCounts;
          int inputLikelihoodCount = 0;
  
          DEBUG2("calculating data likelihoods");
          calculateSampleDataLikelihoods(
              samples,
              results,
              parser,
              genotypesByPloidy,
              parameters,
              usingNull,
              observationBias,
              genotypeAlleles,
              contaminationEstimates,
              estimatedAlleleFrequencies,
              sampleDataLikelihoodsByPopulation,
              variantSampleDataLikelihoodsByPopulation,
              invariantSampleDataLikelihoodsByPopulation);
  
          DEBUG2("finished calculating data likelihoods");
  
  
          // if somehow we get here without any possible sample genotype likelihoods, bail out
          bool hasSampleLikelihoods = false;
          for (map<string, vector<vector<SampleDataLikelihood> > >::iterator s = sampleDataLikelihoodsByPopulation.begin();
               s != sampleDataLikelihoodsByPopulation.end(); ++s) {
              if (!s->second.empty()) {
                  hasSampleLikelihoods = true;
                  break;
              }
          }
          if (!hasSampleLikelihoods) {
              continue;
          }
  
          DEBUG2("calulating combo posteriors over " << parser->populationSamples.size() << " populations");
  
          // XXX
          // TODO skip these steps in the case that there is only one population?
  
  
          // we provide p(var|data), or the probability that the location has
          // variation between individuals relative to the probability that it
          // has no variation
          //
          // in other words:
          // p(var|d) = 1 - p(AA|d) - p(TT|d) - P(GG|d) - P(CC|d)
          //
          // the approach is go through all the homozygous combos
          // and then subtract this from 1... resolving p(var|d)
  
          BigFloat pVar = 1.0;
          BigFloat pHom = 0.0;
  
          long double bestComboOddsRatio = 0;
  
          bool bestOverallComboIsHet = false;
          GenotypeCombo bestCombo; // = NULL;
  
          GenotypeCombo bestGenotypeComboByMarginals;
          vector<vector<SampleDataLikelihood> > allSampleDataLikelihoods;
  
          DEBUG("searching genotype space");
  
          // resample the posterior, this time without bounds on the
          // samples we vary, ensuring that we can generate marginals for
          // all sample/genotype combinations
  
          //SampleDataLikelihoods marginalLikelihoods = sampleDataLikelihoods;  // heavyweight copy...
          map<string, list<GenotypeCombo> > genotypeCombosByPopulation;
          int genotypingTotalIterations = 0; // tally total iterations required to reach convergence
          map<string, list<GenotypeCombo> > glMaxCombos;
  
          for (map<string, SampleDataLikelihoods>::iterator p = sampleDataLikelihoodsByPopulation.begin(); p != sampleDataLikelihoodsByPopulation.end(); ++p) {
  
              const string& population = p->first;
              SampleDataLikelihoods& sampleDataLikelihoods = p->second;
              list<GenotypeCombo>& populationGenotypeCombos = genotypeCombosByPopulation[population];
  
              DEBUG2("genqerating banded genotype combinations from " << sampleDataLikelihoods.size() << " sample genotypes in population " << population);
  
              // cap the number of iterations at 2 x the number of alternate alleles
              // max it at parameters.genotypingMaxIterations iterations, min at 10
              int itermax = min(max(10, 2 * estimatedMinorAllelesAtLocus), parameters.genotypingMaxIterations);
              //int itermax = parameters.genotypingMaxIterations;
  
              // XXX HACK
              // passing 0 for bandwidth and banddepth means "exhaustive local search"
              // this produces properly normalized GQ's at polyallelic sites
              int adjustedBandwidth = 0;
              int adjustedBanddepth = 0;
              // however, this can lead to huge performance problems at complex sites,
              // so we implement this hack...
              if (parameters.genotypingMaxBandDepth > 0 &&
                  genotypeAlleles.size() > parameters.genotypingMaxBandDepth) {
                  adjustedBandwidth = 1;
                  adjustedBanddepth = parameters.genotypingMaxBandDepth;
              }
  
              GenotypeCombo nullCombo;
              SampleDataLikelihoods nullSampleDataLikelihoods;
  
              // this is the genotype-likelihood maximum
              if (parameters.reportGenotypeLikelihoodMax) {
                  GenotypeCombo comboKing;
                  vector<int> initialPosition;
                  initialPosition.assign(sampleDataLikelihoods.size(), 0);
                  SampleDataLikelihoods nullDataLikelihoods; // dummy variable
                  makeComboByDatalLikelihoodRank(comboKing,
                                                 initialPosition,
                                                 sampleDataLikelihoods,
                                                 nullDataLikelihoods,
                                                 inputAlleleCounts,
                                                 theta,
                                                 parameters.pooledDiscrete,
                                                 parameters.ewensPriors,
                                                 parameters.permute,
                                                 parameters.hwePriors,
                                                 parameters.obsBinomialPriors,
                                                 parameters.alleleBalancePriors,
                                                 parameters.diffusionPriorScalar);
  
                  glMaxCombos[population].push_back(comboKing);
              }
  
              // search much longer for convergence
              convergentGenotypeComboSearch(
                  populationGenotypeCombos,
                  nullCombo,
                  sampleDataLikelihoods, // vary everything
                  sampleDataLikelihoods,
                  nullSampleDataLikelihoods,
                  samples,
                  genotypeAlleles,
                  inputAlleleCounts,
                  adjustedBandwidth,
                  adjustedBanddepth,
                  theta,
                  parameters.pooledDiscrete,
                  parameters.ewensPriors,
                  parameters.permute,
                  parameters.hwePriors,
                  parameters.obsBinomialPriors,
                  parameters.alleleBalancePriors,
                  parameters.diffusionPriorScalar,
                  itermax,
                  genotypingTotalIterations,
                  true); // add homozygous combos
                  // ^^ combo results are sorted by default
          }
  
          // generate the GL max combo
          GenotypeCombo glMax;
          if (parameters.reportGenotypeLikelihoodMax) {
              list<GenotypeCombo> glMaxGenotypeCombos;
              combinePopulationCombos(glMaxGenotypeCombos, glMaxCombos);
              glMax = glMaxGenotypeCombos.front();
          }
  
          // accumulate combos from independently-calculated populations into the list of combos
          list<GenotypeCombo> genotypeCombos; // build new combos into this list
          combinePopulationCombos(genotypeCombos, genotypeCombosByPopulation);
          // TODO factor out the following blocks as they are repeated from above
  
          // re-get posterior normalizer
          vector<long double> comboProbs;
          for (list<GenotypeCombo>::iterator gc = genotypeCombos.begin(); gc != genotypeCombos.end(); ++gc) {
              comboProbs.push_back(gc->posteriorProb);
          }
          long double posteriorNormalizer = logsumexp_probs(comboProbs);
  
          // recalculate posterior normalizer
          pVar = 1.0;
          pHom = 0.0;
          // calculates pvar and gets the best het combo
          list<GenotypeCombo>::iterator gc = genotypeCombos.begin();
          bestCombo = *gc;
          for ( ; gc != genotypeCombos.end(); ++gc) {
              if (gc->isHomozygous() && gc->alleles().front() == referenceBase) {
                  pVar -= big_exp(gc->posteriorProb - posteriorNormalizer);
                  pHom += big_exp(gc->posteriorProb - posteriorNormalizer);
              } else if (gc == genotypeCombos.begin()) {
                  bestOverallComboIsHet = true;
              }
          }
  
          // odds ratio between the first and second-best combinations
          if (genotypeCombos.size() > 1) {
              bestComboOddsRatio = genotypeCombos.front().posteriorProb - (++genotypeCombos.begin())->posteriorProb;
          }
  
          if (parameters.calculateMarginals) {
              // make a combined, all-populations sample data likelihoods vector to accumulate marginals
              SampleDataLikelihoods allSampleDataLikelihoods;
              for (map<string, SampleDataLikelihoods>::iterator p = sampleDataLikelihoodsByPopulation.begin(); p != sampleDataLikelihoodsByPopulation.end(); ++p) {
                  SampleDataLikelihoods& sdls = p->second;
                  allSampleDataLikelihoods.reserve(allSampleDataLikelihoods.size() + distance(sdls.begin(), sdls.end()));
                  allSampleDataLikelihoods.insert(allSampleDataLikelihoods.end(), sdls.begin(), sdls.end());
              }
              // calculate the marginal likelihoods for this population
              marginalGenotypeLikelihoods(genotypeCombos, allSampleDataLikelihoods);
              // store the marginal data likelihoods in the results, for easy parsing
              // like a vector -> map conversion...
              results.update(allSampleDataLikelihoods);
          }
  
          map<string, int> repeats;
          if (parameters.showReferenceRepeats) {
              repeats = parser->repeatCounts(parser->currentSequencePosition(), parser->currentSequence, 12);
          }
  
          vector<Allele> alts;
          if (parameters.onlyUseInputAlleles
              || parameters.reportAllHaplotypeAlleles
              || parameters.pooledContinuous) {
              //alts = genotypeAlleles;
              for (vector<Allele>::iterator a = genotypeAlleles.begin(); a != genotypeAlleles.end(); ++a) {
                  if (!a->isReference()) {
                      alts.push_back(*a);
                  }
              }
          } else {
              // get the unique alternate alleles in this combo, sorted by frequency in the combo
              vector<pair<Allele, int> > alternates = alternateAlleles(bestCombo, referenceBase);
              for (vector<pair<Allele, int> >::iterator a = alternates.begin(); a != alternates.end(); ++a) {
                  Allele& alt = a->first;
                  if (!alt.isNull() && !alt.isReference())
                      alts.push_back(alt);
              }
              // if there are no alternate alleles in the best combo, use the genotype alleles
              // XXX ...
              if (alts.empty()) {
                  for (vector<Allele>::iterator a = genotypeAlleles.begin(); a != genotypeAlleles.end(); ++a) {
                      if (!a->isReference()) {
                          alts.push_back(*a);
                      }
                  }
              }
          }
  
          // reporting the GL maximum *over all alleles*
          if (parameters.reportGenotypeLikelihoodMax) {
              bestCombo = glMax;
          } else {
              // the default behavior is to report the GL maximum genotyping over the alleles in the best posterior genotyping
  
              // select the maximum-likelihood GL given the alternates we have
              // this is not the same thing as the GL max over all alleles!
              // it is the GL max over the selected alleles at this point
  
              vector<Allele> alleles = alts;
              for (vector<Allele>::iterator a = genotypeAlleles.begin(); a != genotypeAlleles.end(); ++a) {
                  if (a->isReference()) {
                      alleles.push_back(*a);
                  }
              }
              map<string, list<GenotypeCombo> > glMaxComboBasedOnAltsByPop;
              for (map<string, SampleDataLikelihoods>::iterator p = sampleDataLikelihoodsByPopulation.begin(); p != sampleDataLikelihoodsByPopulation.end(); ++p) {
                  const string& population = p->first;
                  SampleDataLikelihoods& sampleDataLikelihoods = p->second;
                  GenotypeCombo glMaxBasedOnAlts;
                  for (SampleDataLikelihoods::iterator v = sampleDataLikelihoods.begin(); v != sampleDataLikelihoods.end(); ++v) {
                      SampleDataLikelihood* m = NULL;
                      for (vector<SampleDataLikelihood>::iterator d = v->begin(); d != v->end(); ++d) {
                          if (d->genotype->matchesAlleles(alleles)) {
                              m = &*d;
                              break;
                          }
                      }
                      assert(m != NULL);
                      glMaxBasedOnAlts.push_back(m);
                  }
                  glMaxComboBasedOnAltsByPop[population].push_back(glMaxBasedOnAlts);
              }
              list<GenotypeCombo> glMaxBasedOnAltsGenotypeCombos; // build new combos into this list
              combinePopulationCombos(glMaxBasedOnAltsGenotypeCombos, glMaxComboBasedOnAltsByPop);
              bestCombo = glMaxBasedOnAltsGenotypeCombos.front();
          }
  
          DEBUG("best combo: " << bestCombo);
  
          // output
  
          if ((!alts.empty() && (1 - pHom.ToDouble()) >= parameters.PVL) || parameters.PVL == 0){
  
              // write the last gVCF record(s)
              if (parameters.gVCFout && !nonCalls.empty()) {
                  vcflib::Variant var(parser->variantCallFile);
                  out << results.gvcf(var, nonCalls, parser) << endl;
                  nonCalls.clear();
              }
  
              vcflib::Variant var(parser->variantCallFile);
  
              out << results.vcf(
                  var,
                  pHom,
                  bestComboOddsRatio,
                  samples,
                  referenceBase,
                  alts,
                  repeats,
                  genotypingTotalIterations,
                  parser->sampleList,
                  coverage,
                  bestCombo,
                  alleleGroups,
                  partialObservationGroups,
                  partialObservationSupport,
                  genotypesByPloidy,
                  parser->sequencingTechnologies,
                  parser)
                  << endl;
  
          } else if (parameters.gVCFout) {
              // record statistics for gVCF output
              nonCalls.record(parser->currentSequenceName, parser->currentPosition, samples);
          }
          DEBUG2("finished position");
  
      }
  
      // write the last gVCF record
      // NOTE: for some resion this is only needed if we are using gVCF chunks, if minimal chunking it is not requird, in fact it breaks....
      if (parameters.gVCFout && !nonCalls.empty() && !parameters.gVCFNoChunk) {
          Results results;
          vcflib::Variant var(parser->variantCallFile);
          out << results.gvcf(var, nonCalls, parser) << endl;
          nonCalls.clear();
      }
  
      DEBUG("total sites: " << total_sites << endl
            << "processed sites: " << processed_sites << endl
            << "ratio: " << (float) processed_sites / (float) total_sites);
  
      delete parser;
  
      return 0;
  
  }
  