text: // Author: Mark Chaisson
  #pragma once
  
  #include "BlasrAlign.hpp"
  
  //----------------------MODIFY ALIGNMENTS--------------------------//
  void AssignRefContigLocation(T_AlignmentCandidate &alignment,
                               SequenceIndexDatabase<FASTQSequence> &seqdb, DNASequence &genome)
  {
      //
      // If the sequence database is used, the start position of
      // the alignment is relative to the start of the chromosome,
      // not the entire index.  Subtract off the start position of
      // the chromosome to get the true position.
      //
      DNALength forwardTPos;
      int seqDBIndex;
      if (alignment.tStrand == 0) {
          forwardTPos = alignment.tAlignedSeqPos;
          seqDBIndex = seqdb.SearchForIndex(forwardTPos);
          alignment.tAlignedSeqPos -= seqdb.seqStartPos[seqDBIndex];
      } else {
          //
          // Flip coordinates into forward strand in order to find the boundaries
          // of the contig, then reverse them in order to find offset.
          //
  
          // Find the reverse complement coordinate of the index of the last aligned base.
          assert(alignment.tAlignedSeqLength > 0);
          forwardTPos =
              genome.MakeRCCoordinate(alignment.tAlignedSeqPos + alignment.tAlignedSeqLength - 1);
          seqDBIndex = seqdb.SearchForIndex(forwardTPos);
  
          //
          // Find the reverse comlement coordinate of the last base of this
          // sequence.  This would normally be the start of the next contig
          // -1 to get the length, but since an 'N' is added between every
          // pair of sequences, this is -2.
          //
          DNALength reverseTOffset;
          reverseTOffset = genome.MakeRCCoordinate(seqdb.seqStartPos[seqDBIndex + 1] - 2);
          alignment.tAlignedSeqPos -= reverseTOffset;
      }
  }
  
  void AssignRefContigLocations(std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                SequenceIndexDatabase<FASTQSequence> &seqdb, DNASequence &genome)
  {
  
      UInt i;
      for (i = 0; i < alignmentPtrs.size(); i++) {
          T_AlignmentCandidate *aref = alignmentPtrs[i];
          AssignRefContigLocation(*aref, seqdb, genome);
      }
  }
  
  template <typename T_RefSequence>
  void AssignGenericRefContigName(std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                  T_RefSequence &genome)
  {
      UInt i;
      for (i = 0; i < alignmentPtrs.size(); i++) {
          T_AlignmentCandidate *aref = alignmentPtrs[i];
          aref->tName = genome.title;
      }
  }
  
  void StoreRankingStats(std::vector<T_AlignmentCandidate *> &alignments,
                         VarianceAccumulator<float> &accumPValue,
                         VarianceAccumulator<float> &accumWeight)
  {
      int i;
      for (i = 0; i < int(alignments.size()); i++) {
          alignments[i]->pvalVariance = accumPValue.GetVariance();
          alignments[i]->pvalNStdDev = accumPValue.GetNStdDev(alignments[i]->clusterScore);
          alignments[i]->weightVariance = accumWeight.GetVariance();
          alignments[i]->weightNStdDev = accumWeight.GetNStdDev(alignments[i]->clusterWeight);
      }
  }
  
  void AssignMapQV(std::vector<T_AlignmentCandidate *> &alignmentPtrs)
  {
      int i;
      int mapQV = 1;
      if (alignmentPtrs.size() > 1 and alignmentPtrs[0]->score == alignmentPtrs[1]->score) {
          // the top two alignments have the same score, don't consider them as mapped.
          mapQV = 0;
      }
  
      for (i = 0; i < int(alignmentPtrs.size()); i++) {
          alignmentPtrs[i]->mapQV = mapQV;
      }
  }
  
  void ScaleMapQVByClusterSize(T_AlignmentCandidate &alignment, MappingParameters &params)
  {
      if (alignment.numSignificantClusters > int(params.nCandidates)) {
          alignment.mapQV = Phred((1 - InversePhred(alignment.mapQV)) *
                                  ((float)params.nCandidates / alignment.numSignificantClusters));
      } else if (alignment.numSignificantClusters == 0) {
          alignment.mapQV = 0;
      }
  }
  
  void StoreMapQVs(SMRTSequence &read, std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                   MappingParameters &params)
  {
      //
      // Only weight alignments for mapqv against eachother if they are overlapping.
      //
      int a;
      std::vector<std::set<int> >
          partitions;  // Each set contains alignments that overlap on the read.
      DistanceMatrixScoreFunction<DNASequence, FASTQSequence> distScoreFn;
      distScoreFn.del = params.deletion;
      distScoreFn.ins = params.insertion;
      // bug 24363, set affineOpen and affineExtend for distScoreFn
      distScoreFn.affineOpen = params.affineOpen;
      distScoreFn.affineExtend = params.affineExtend;
      distScoreFn.InitializeScoreMatrix(SMRTLogProbMatrix);
      IDSScoreFunction<DNASequence, FASTQSequence> idsScoreFn;
      idsScoreFn.ins = params.insertion;
      idsScoreFn.del = params.deletion;
      idsScoreFn.affineExtend = params.affineExtend;
      idsScoreFn.affineOpen = params.affineOpen;
      idsScoreFn.substitutionPrior = params.substitutionPrior;
      idsScoreFn.globalDeletionPrior = params.globalDeletionPrior;
  
      //
      // Rescore the alignment so that it uses probabilities.
      //
      for (a = 0; a < int(alignmentPtrs.size()); a++) {
          if (params.ignoreQualities == false) {
              // bug 24363, pass -affineAlign to compute correct alignment score.
              alignmentPtrs[a]->probScore =
                  -ComputeAlignmentScore(*alignmentPtrs[a], alignmentPtrs[a]->qAlignedSeq,
                                         alignmentPtrs[a]->tAlignedSeq, idsScoreFn,
                                         params.affineAlign) /
                  10.0;
          } else {
              alignmentPtrs[a]->probScore =
                  -ComputeAlignmentScore(*alignmentPtrs[a], alignmentPtrs[a]->qAlignedSeq,
                                         alignmentPtrs[a]->tAlignedSeq, distScoreFn,
                                         params.affineAlign) /
                  10.0;
          }
      }
      PartitionOverlappingAlignments(alignmentPtrs, partitions,
                                     params.minFractionToBeConsideredOverlapping);
  
      int p;
      std::set<int>::iterator partIt, partEnd;
  
      //
      // For each partition, store where on the read it begins, and where
      // it ends.
      //
      std::vector<int> partitionBeginPos, partitionEndPos;
      partitionBeginPos.resize(partitions.size());
      partitionEndPos.resize(partitions.size());
      std::fill(partitionBeginPos.begin(), partitionBeginPos.end(), -1);
      std::fill(partitionEndPos.begin(), partitionEndPos.end(), -1);
      std::vector<char> assigned;
      assigned.resize(alignmentPtrs.size());
      std::fill(assigned.begin(), assigned.end(), false);
  
      for (p = 0; p < int(partitions.size()); p++) {
          partEnd = partitions[p].end();
          int alnStart, alnEnd;
  
          if (partitions[p].size() > 0) {
              partIt = partitions[p].begin();
              alignmentPtrs[*partIt]->GetQInterval(alnStart, alnEnd);
              partitionBeginPos[p] = alnStart;
              partitionEndPos[p] = alnEnd;
              ++partIt;
              partEnd = partitions[p].end();
              for (; partIt != partEnd; ++partIt) {
                  //  Comment out because all reads are now in the forward strand.
                  //  alignmentPtrs[*partIt]->GetQInterval(alnStart, alnEnd, convertToForwardStrand);
                  alignmentPtrs[*partIt]->GetQInterval(alnStart, alnEnd);
                  if (alnEnd - alnStart > partitionEndPos[p] - partitionBeginPos[p]) {
                      partitionBeginPos[p] = alnStart;
                      partitionEndPos[p] = alnEnd;
                  }
              }
          }
      }
  
      //
      // For each partition, determine the widest parts of the read that
      // are aligned in the partition.  All alignments will be extended to
      // the end of the widest parts of the partition.
      //
      const static bool convertToForwardStrand = true;
  
      UInt i;
  
      //
      // For now, just use the alignment score as the probability score.
      // Although it is possible to use the full forward probability, for
      // the most part it is pretty much the same as the Vitterbi
      // probability, but it takes a lot longer to compute.
      //
  
      //
      // Now estimate what the alignment scores would be if they were
      // extended past the ends of their current alignment.
      //
  
      for (p = 0; p < int(partitions.size()); p++) {
          partEnd = partitions[p].end();
          int alnStart, alnEnd;
          for (partIt = partitions[p].begin(); partitions[p].size() > 0 and partIt != partEnd;
               ++partIt) {
              int mismatchSum = 0;
              alignmentPtrs[*partIt]->GetQInterval(alnStart, alnEnd, convertToForwardStrand);
              if (alnStart - partitionBeginPos[p] > MAPQV_END_ALIGN_WIGGLE or
                  partitionEndPos[p] - alnEnd > MAPQV_END_ALIGN_WIGGLE) {
                  // bug 24363, use updated SumMismatches to compute mismatch score when
                  // no QV is available.
                  SumMismatches(read, *alignmentPtrs[*partIt], 15, partitionBeginPos[p],
                                partitionEndPos[p], params, mismatchSum);
              }
              //
              // Random sequence can be aligned with about 50% similarity due
              // to optimization, so weight the qv sum
              //
              alignmentPtrs[*partIt]->probScore += -(mismatchSum)*0.5;
          }
      }
  
      //
      // Determine mapqv by summing qvscores in partitions
  
      float mapQVDenominator = 0;
      for (p = 0; p < int(partitions.size()); p++) {
          std::set<int>::iterator nextIt;
          if (partitions[p].size() == 0) {
              continue;
          }
          int index = *partitions[p].begin();
  
          mapQVDenominator = alignmentPtrs[index]->probScore;
  
          if (partitions[p].size() > 1) {
              partIt = partitions[p].begin();
              partEnd = partitions[p].end();
              ++partIt;
  
              for (; partIt != partEnd; ++partIt) {
                  index = *partIt;
                  mapQVDenominator = LogSumOfTwo(mapQVDenominator, alignmentPtrs[index]->probScore);
              }
          }
  
          for (partIt = partitions[p].begin(); partIt != partitions[p].end(); ++partIt) {
              //
              // If only one alignment is found, assume maximum mapqv.
              //
              assigned[*partIt] = true;
              if (partitions[p].size() == 1) {
                  alignmentPtrs[*partIt]->mapQV = MAX_PHRED_SCORE;
              }
  
              //
              // Look for overflow.
              //
              else if (alignmentPtrs[*partIt]->probScore - mapQVDenominator < -20) {
                  alignmentPtrs[*partIt]->mapQV = 0;
              } else {
                  double log10 = log(10);
                  double sub = alignmentPtrs[*partIt]->probScore - mapQVDenominator;
                  double expo = exp(log10 * sub);
                  double diff = 1.0 - expo;
                  int phredValue;
  
                  if (expo == 0) {
                      phredValue = 0;
                  } else if (diff == 0) {
                      phredValue = MAX_PHRED_SCORE;
                  } else {
                      phredValue = Phred(diff);
                  }
                  if (phredValue > MAX_PHRED_SCORE) {
                      phredValue = MAX_PHRED_SCORE;
                  }
  
                  alignmentPtrs[*partIt]->mapQV = phredValue;
                  assigned[*partIt] = true;
              }
  
              if (params.scaleMapQVByNumSignificantClusters) {
                  ScaleMapQVByClusterSize(*alignmentPtrs[*partIt], params);
              }
          }
      }
  
      for (i = 0; i < assigned.size(); i++) {
          assert(assigned[i]);
      }
  }
  
  //--------------------SEARCH & CHECK ALIGNMENTS-------------------//
  template <typename T_Sequence>
  bool CheckForSufficientMatch(T_Sequence &read, std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                               MappingParameters &params)
  {
      (void)(read);
      if (alignmentPtrs.size() > 0 and alignmentPtrs[0]->score < params.maxScore) {
          return true;
      } else {
          return false;
      }
  }
  
  int FindMaxLengthAlignment(std::vector<T_AlignmentCandidate *> alignmentPtrs, int &maxLengthIndex)
  {
      int i;
      int maxLength = 0;
      maxLengthIndex = -1;
  
      for (i = 0; i < int(alignmentPtrs.size()); i++) {
          int qStart, qEnd;
          alignmentPtrs[i]->GetQInterval(qStart, qEnd);
          if (qEnd - qStart > maxLength) {
              maxLengthIndex = i;
              maxLength = qEnd - qStart;
          }
      }
      return (maxLength != -1);
  }
  
  void SumMismatches(SMRTSequence &read, T_AlignmentCandidate &alignment, int mismatchScore,
                     int fullIntvStart, int fullIntvEnd, MappingParameters &params, int &sum)
  {
      int alnStart, alnEnd;
      alignment.GetQIntervalOnForwardStrand(alnStart, alnEnd);
      int p;
      sum = 0;
      if (not params.ignoreQualities and read.substitutionQV.Empty() == false) {
          for (p = fullIntvStart; p < alnStart; p++) {
              sum += read.substitutionQV[p];
          }
          for (p = alnEnd; p < fullIntvEnd; p++) {
              sum += read.substitutionQV[p];
          }
      } else {
          // bug 24363, compute mismatch score when QV is not available.
          sum += mismatchScore * ((alnStart - fullIntvStart) + (fullIntvEnd - alnEnd));
      }
  }
  
  bool AlignmentsOverlap(T_AlignmentCandidate &alnA, T_AlignmentCandidate &alnB,
                         float minPercentOverlap)
  {
      int alnAStart, alnAEnd, alnBStart, alnBEnd;
      bool useForwardStrand = true;
      alnA.GetQInterval(alnAStart, alnAEnd, useForwardStrand);
      alnB.GetQInterval(alnBStart, alnBEnd, useForwardStrand);
      // Look if one alignment encompasses the other
      int ovp = 0;
      if (alnAStart <= alnBStart and alnAEnd >= alnBEnd) {
          return true;
      } else if (alnBStart <= alnAStart and alnBEnd >= alnAEnd) {
          return true;
          //ovp = alnAEnd - alnAStart;
      } else {
          //
          // Look to see if the alignments overlap
          //
  
          if (alnAEnd >= alnBStart and alnAEnd <= alnBEnd) {
              ovp = alnAEnd - alnBStart;
          } else if (alnAStart >= alnBStart and alnAStart <= alnBEnd) {
              ovp = alnBEnd - alnAStart;
          }
      }
  
      // float ovpPercent = (2.0*ovp) / ((alnAEnd - alnAStart) + (alnBEnd - alnBStart));
      float ovpPercent = 0;
      if (alnAEnd - alnAStart > 0 and alnBEnd - alnBStart > 0) {
          // overlap percentage: maximum overlap percent in A and B.
          ovpPercent = std::max(float(ovp) / float(alnAEnd - alnAStart),
                                float(ovp) / float(alnBEnd - alnBStart));
      }
  
      // returns true when an overlap is found.
      return (ovpPercent > minPercentOverlap);
  }
  
  void PartitionOverlappingAlignments(std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                      std::vector<std::set<int> > &partitions, float minOverlap)
  {
      if (alignmentPtrs.size() == 0) {
          partitions.clear();
          return;
      }
  
      std::set<int>::iterator setIt, setEnd;
      int i, p;
      bool overlapFound = false;
      for (i = 0; i < int(alignmentPtrs.size()); i++) {
          overlapFound = false;
          for (p = 0; p < int(partitions.size()) and overlapFound == false; p++) {
              setEnd = partitions[p].end();
              for (setIt = partitions[p].begin();
                   setIt != partitions[p].end() and overlapFound == false; ++setIt) {
                  if (AlignmentsOverlap(*alignmentPtrs[i], *alignmentPtrs[*setIt], minOverlap) or
                      ((alignmentPtrs[i]->QAlignStart() <= alignmentPtrs[*setIt]->QAlignStart()) and
                       (alignmentPtrs[i]->QAlignEnd() > alignmentPtrs[*setIt]->QAlignEnd()))) {
                      partitions[p].insert(i);
                      overlapFound = true;
                  }
              }
          }
          //
          // If this alignment does not overlap any other, create a
          // partition with it as the first element.
          //
          if (overlapFound == false) {
              partitions.push_back(std::set<int>());
              partitions[partitions.size() - 1].insert(i);
          }
      }
  }
  
  //--------------------FILTER ALIGNMENTS---------------------------//
  int RemoveLowQualitySDPAlignments(int readLength,
                                    std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                    MappingParameters &params)
  {
      // Just a hack.  For now, assume there is at least 1 match per 50 bases.
      int totalBasesMatched = 0;
      int a;
      for (a = 0; a < int(alignmentPtrs.size()); a++) {
          int b;
          for (b = 0; b < int(alignmentPtrs[a]->blocks.size()); b++) {
              totalBasesMatched += alignmentPtrs[a]->blocks[b].length;
          }
          int expectedMatches = params.sdpTupleSize / 50.0 * readLength;
          if (totalBasesMatched < expectedMatches) {
              delete alignmentPtrs[a];
              alignmentPtrs[a] = NULL;
          }
      }
      int packedAlignmentIndex = 0;
      for (a = 0; a < int(alignmentPtrs.size()); a++) {
          if (alignmentPtrs[a] != NULL) {
              alignmentPtrs[packedAlignmentIndex] = alignmentPtrs[a];
              packedAlignmentIndex++;
          }
      }
      alignmentPtrs.resize(packedAlignmentIndex);
      return packedAlignmentIndex;
  }
  
  template <typename T_Sequence>
  int RemoveLowQualityAlignments(T_Sequence &read, std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                 MappingParameters &params)
  {
      PB_UNUSED(read);
      if (params.verbosity > 0) {
          std::cout << "checking at least " << alignmentPtrs.size()
                    << " alignments to see if they are accurate." << std::endl;
      }
      for (size_t i = 0; i < MIN(static_cast<size_t>(params.nCandidates), alignmentPtrs.size());
           i++) {
          if (params.verbosity > 0) {
              std::cout << "Quality check  " << i << " " << alignmentPtrs[i]->score << std::endl;
          }
          if (alignmentPtrs[i]->blocks.size() == 0 or alignmentPtrs[i]->score > params.maxScore) {
              //
              // Since the alignments are sorted according to alignment
              // score, once one of the alignments is too low of a score,
              // all remaining alignments are also too low, and should be
              // removed as well.  Do that all at once.
              //
              if (alignmentPtrs[i]->blocks.size() == 0 and params.verbosity > 0) {
                  std::cout << "Removing empty alignment " << alignmentPtrs[i]->qName << std::endl;
              }
              if (params.verbosity > 0) {
                  std::cout << alignmentPtrs[i]->qName << " alignment " << i
                            << " is too low of a score." << alignmentPtrs[i]->score << std::endl;
              }
              for (size_t deletedIndex = i; deletedIndex < alignmentPtrs.size(); deletedIndex++) {
                  delete alignmentPtrs[deletedIndex];
                  alignmentPtrs[deletedIndex] = NULL;
              }
              alignmentPtrs.erase(i + alignmentPtrs.begin(), alignmentPtrs.end());
              break;
          } else {
              if (params.verbosity > 0) {
                  std::cout << "Keeping alignment " << i << " " << alignmentPtrs[i]->qPos << " "
                            << alignmentPtrs[i]->qLength << " " << alignmentPtrs[i]->tName << " "
                            << alignmentPtrs[i]->tPos << " " << alignmentPtrs[i]->tLength
                            << " from score: " << alignmentPtrs[i]->score << std::endl;
              }
          }
      }
      return alignmentPtrs.size();
  }
  
  //FIXME: move to class ReadAlignments
  int RemoveOverlappingAlignments(std::vector<T_AlignmentCandidate *> &alignmentPtrs,
                                  MappingParameters &params)
  {
      std::vector<unsigned char> alignmentIsContained;
      alignmentIsContained.resize(alignmentPtrs.size());
      std::fill(alignmentIsContained.begin(), alignmentIsContained.end(), false);
  
      int j;
      int numContained = 0;
      int curNotContained = 0;
  
      if (alignmentPtrs.size() > 0) {
          UInt i;
          for (i = 0; i < alignmentPtrs.size() - 1; i++) {
              T_AlignmentCandidate *aref = alignmentPtrs[i];
              if (aref->pctSimilarity < params.minPctSimilarity) {
                  continue;
              }
              for (j = i + 1; j < int(alignmentPtrs.size()); j++) {
                  //
                  // Make sure this alignment isn't already removed.
                  //
                  if (alignmentIsContained[j]) {
                      continue;
                  }
  
                  //
                  // Only check for containment if the two sequences are from the same contig.
                  //
                  if (alignmentPtrs[i]->tIndex != alignmentPtrs[j]->tIndex) {
                      continue;
                  }
  
                  //
                  // Check for an alignment that is fully overlapping another
                  // alignment.
                  if (aref->GenomicTBegin() <= alignmentPtrs[j]->GenomicTBegin() and
                      aref->GenomicTEnd() >= alignmentPtrs[j]->GenomicTEnd() and
                      alignmentPtrs[i]->tIndex == alignmentPtrs[j]->tIndex) {
                      //
                      // Alignment i is contained in j is only true if it has a worse score.
                      //
                      if (aref->score <= alignmentPtrs[j]->score) {
                          alignmentIsContained[j] = true;
                      }
                      if (params.verbosity >= 2) {
                          std::cout << "alignment " << i << " is contained in " << j << std::endl;
                          std::cout << aref->tAlignedSeqPos << " " << alignmentPtrs[j]->tAlignedSeqPos
                                    << " " << aref->tAlignedSeqPos + aref->tAlignedSeqLength << " "
                                    << alignmentPtrs[j]->tAlignedSeqPos +
                                           alignmentPtrs[j]->tAlignedSeqLength
                                    << std::endl;
                      }
                  } else if (alignmentPtrs[j]->GenomicTBegin() <= aref->GenomicTBegin() and
                             alignmentPtrs[j]->GenomicTEnd() >= aref->GenomicTEnd() and
                             alignmentPtrs[i]->tIndex == alignmentPtrs[j]->tIndex) {
                      if (params.verbosity >= 2) {
                          std::cout << "ALIGNMENT " << j << " is contained in " << i << std::endl;
                          std::cout << alignmentPtrs[j]->tAlignedSeqPos << " " << aref->tAlignedSeqPos
                                    << " "
                                    << alignmentPtrs[j]->tAlignedSeqPos +
                                           alignmentPtrs[j]->tAlignedSeqLength
                                    << " " << aref->tAlignedSeqPos + aref->tAlignedSeqLength
                                    << std::endl;
                      }
                      if (alignmentPtrs[j]->score <= aref->score) {
                          alignmentIsContained[i] = true;
                      }
                  }
              }
          }
          for (i = 0; i < alignmentPtrs.size(); i++) {
              T_AlignmentCandidate *aref = alignmentPtrs[i];
              if (alignmentIsContained[i]) {
                  delete alignmentPtrs[i];
                  alignmentPtrs[i] = NULL;
                  numContained++;
              } else {
                  alignmentPtrs[curNotContained] = aref;
                  ++curNotContained;
              }
          }
          alignmentPtrs.resize(alignmentPtrs.size() - numContained);
      }
      return alignmentPtrs.size();
  }
  
  // Delete all alignments from index startIndex in vector, inclusive.
  void DeleteAlignments(std::vector<T_AlignmentCandidate *> &alignmentPtrs, int startIndex)
  {
      int i;
      for (i = startIndex; i < int(alignmentPtrs.size()); i++) {
          delete alignmentPtrs[i];
      }
      alignmentPtrs.resize(0);
  }
  
  //--------------------REFINE ALIGNMENTS---------------------------//
  template <typename T_RefSequence, typename T_Sequence>
  void RefineAlignment(std::vector<T_Sequence *> &bothQueryStrands, T_RefSequence &genome,
                       T_AlignmentCandidate &alignmentCandidate, MappingParameters &params,
                       MappingBuffers &mappingBuffers)
  {
      (void)(genome);
      FASTQSequence qSeq;
      DNASequence tSeq;
      DistanceMatrixScoreFunction<DNASequence, FASTQSequence> distScoreFn(
          SMRTDistanceMatrix, params.deletion, params.insertion);
  
      DistanceMatrixScoreFunction<DNASequence, FASTQSequence> distScoreFn2(
          SMRTDistanceMatrix, params.indel, params.indel);
  
      QualityValueScoreFunction<DNASequence, FASTQSequence> scoreFn;
      IDSScoreFunction<DNASequence, FASTQSequence> idsScoreFn;
      idsScoreFn.InitializeScoreMatrix(SMRTDistanceMatrix);
      scoreFn.del = params.indel;
      scoreFn.ins = params.indel;
      idsScoreFn.ins = params.insertion;
      idsScoreFn.del = params.deletion;
      idsScoreFn.affineExtend = params.affineExtend;
      idsScoreFn.affineOpen = params.affineOpen;
      idsScoreFn.substitutionPrior = params.substitutionPrior;
      idsScoreFn.globalDeletionPrior = params.globalDeletionPrior;
  
      if (params.doGlobalAlignment) {
          // global and placeGapConsistently can not set at the same time
          assert(not params.placeGapConsistently);
          SMRTSequence subread;
          subread.ReferenceSubstring(*bothQueryStrands[0], bothQueryStrands[0]->SubreadStart(),
                                     (bothQueryStrands[0]->SubreadLength()));
  
          int drift = ComputeDrift(alignmentCandidate);
          T_AlignmentCandidate refinedAlignment;
  
          KBandAlign(subread, alignmentCandidate.tAlignedSeq, SMRTDistanceMatrix, params.insertion,
                     params.deletion, drift, mappingBuffers.scoreMat, mappingBuffers.pathMat,
                     refinedAlignment, idsScoreFn, Global);
          refinedAlignment.RemoveEndGaps();
          ComputeAlignmentStats(refinedAlignment, subread.seq, alignmentCandidate.tAlignedSeq.seq,
                                distScoreFn2);
          //idsScoreFn);
  
          alignmentCandidate.blocks = refinedAlignment.blocks;
          alignmentCandidate.gaps = refinedAlignment.gaps;
          alignmentCandidate.tPos = refinedAlignment.tPos;
          alignmentCandidate.qPos = refinedAlignment.qPos + bothQueryStrands[0]->SubreadStart();
          alignmentCandidate.score = refinedAlignment.score;
          subread.Free();
      } else if (params.useGuidedAlign) {
          T_AlignmentCandidate refinedAlignment;
          int lastBlock = alignmentCandidate.blocks.size() - 1;
  
          if (alignmentCandidate.blocks.size() > 0) {
  
              /*
               * Refine the alignment without expanding past the current
               * boundaries of the sequences that are already aligned.
               */
  
              //
              // NOTE** this only makes sense when
              // alignmentCandidate.blocks[0].tPos == 0. Otherwise the length
              // of the sequence is not correct.
              //
              tSeq.Copy(
                  alignmentCandidate.tAlignedSeq, alignmentCandidate.tPos,
                  (alignmentCandidate.blocks[lastBlock].tPos +
                   alignmentCandidate.blocks[lastBlock].length - alignmentCandidate.blocks[0].tPos));
  
              //      qSeq.ReferenceSubstring(alignmentCandidate.qAlignedSeq,
              qSeq.ReferenceSubstring(*bothQueryStrands[static_cast<size_t>(
                                          (alignmentCandidate.qStrand == Forward) ? 0 : 1)],
                                      alignmentCandidate.qAlignedSeqPos + alignmentCandidate.qPos,
                                      (alignmentCandidate.blocks[lastBlock].qPos +
                                       alignmentCandidate.blocks[lastBlock].length));
  
              assert(not(params.affineAlign and params.placeGapConsistently));
  
              if (!params.ignoreQualities &&
                  ReadHasMeaningfulQualityValues(alignmentCandidate.qAlignedSeq)) {
                  if (params.affineAlign) {
                      AffineGuidedAlign(qSeq, tSeq, alignmentCandidate, idsScoreFn, params.bandSize,
                                        mappingBuffers, refinedAlignment, Global, false);
                  } else {
                      GuidedAlign(qSeq, tSeq, alignmentCandidate, idsScoreFn,
                                  params.guidedAlignBandSize, mappingBuffers, refinedAlignment,
                                  Global, false);
                  }
              } else {
                  if (params.affineAlign) {
                      AffineGuidedAlign(qSeq, tSeq, alignmentCandidate, distScoreFn, params.bandSize,
                                        mappingBuffers, refinedAlignment, Global, false);
                  } else {
                      GuidedAlign(qSeq, tSeq, alignmentCandidate, distScoreFn,
                                  params.guidedAlignBandSize, mappingBuffers, refinedAlignment,
                                  Global, false);
                  }
              }
              ComputeAlignmentStats(refinedAlignment, qSeq.seq, tSeq.seq, distScoreFn2,
                                    params.affineAlign);
              //
              // Copy the refine alignment, which may be a subsequence of the
              // alignmentCandidate into the alignment candidate.
              //
  
              // First copy the alignment block and gap (the description of
              // the base by base alignment).
  
              alignmentCandidate.blocks.clear();
              alignmentCandidate.blocks = refinedAlignment.blocks;
  
              alignmentCandidate.CopyStats(refinedAlignment);
  
              alignmentCandidate.gaps = refinedAlignment.gaps;
              alignmentCandidate.score = refinedAlignment.score;
              alignmentCandidate.nCells = refinedAlignment.nCells;
  
              // Next copy the information that describes what interval was
              // aligned.  Since the reference sequences of the alignment
              // candidate have been modified, they are reassigned.
              alignmentCandidate.tAlignedSeq.Free();
              alignmentCandidate.tAlignedSeq.TakeOwnership(tSeq);
              alignmentCandidate.ReassignQSequence(qSeq);
              alignmentCandidate.tAlignedSeqPos += alignmentCandidate.tPos;
              alignmentCandidate.qAlignedSeqPos += alignmentCandidate.qPos;
  
              //
              // tPos and qPos are the positions within the interval where the
              // alignment begins. The refined alignment has adifferent tPos
              // and qPos from the alignment candidate.
              alignmentCandidate.tPos = refinedAlignment.tPos;
              alignmentCandidate.qPos = refinedAlignment.qPos;
  
              // The lengths of the newly aligned sequences may differ, update those.
              alignmentCandidate.tAlignedSeqLength = tSeq.length;
              alignmentCandidate.qAlignedSeqLength = qSeq.length;
          }
      } else {
          assert(not params.placeGapConsistently);
  
          //
          // This assumes an SDP alignment has been performed to create 'alignmentCandidate'.
  
          //
          // Recompute the alignment using a banded smith waterman to
          // get rid of any spurious effects of usign the seeded gaps.
          //
  
          //
          // The k-banded alignment is over a subsequence of the first
          // (sparse dynamic programming, SDP) alignment.  The SDP
          // alignment is over a large window that may contain the
          // candidate sequence.  The k-band alignment is over a tighter
          // region.
  
          int drift = ComputeDrift(alignmentCandidate);
  
          //
          // Rescore the alignment with a banded alignment that has a
          // better model of sequencing error.
          //
  
          if (alignmentCandidate.blocks.size() == 0) {
              alignmentCandidate.score = 0;
              return;
          }
          int lastBlock = alignmentCandidate.blocks.size() - 1;
  
          //
          // Assign the sequences that are going to be realigned using
          // banded alignment.  The SDP alignment does not give that great
          // of a score, but it does do a good job at finding a backbone
          // alignment that closely defines the sequence that is aligned.
          // Reassign the subsequences for alignment with a tight bound
          // around the beginning and ending of each sequence, so that
          // global banded alignment may be performed.
          //
  
          //
          // This section needs to be cleaned up substantially.  Right now it
          // copies a substring from the ref to a temp, then from the temp
          // back to the ref.  It may be possible to just keep one pointer per
          // read to the memory that was allocated, then allow the seq
          // parameter to float around.  The reason for all the copying is
          // that in case there is a compressed version of the genome the
          // seqences must be transformed before alignment.
          //
  
          if (alignmentCandidate.qIsSubstring) {
              qSeq.ReferenceSubstring(*bothQueryStrands[0],  // the original sequence
                                      alignmentCandidate.qPos + alignmentCandidate.qAlignedSeqPos,
                                      alignmentCandidate.blocks[lastBlock].qPos +
                                          alignmentCandidate.blocks[lastBlock].length);
          } else {
              qSeq.ReferenceSubstring(
                  alignmentCandidate.qAlignedSeq,  // the subsequence that the alignment points to
                  alignmentCandidate.qPos + alignmentCandidate.qAlignedSeqPos,
                  alignmentCandidate.blocks[lastBlock].qPos +
                      alignmentCandidate.blocks[lastBlock].length -
                      alignmentCandidate.blocks[0].qPos);
          }
  
          tSeq.Copy(alignmentCandidate.tAlignedSeq,  // the subsequence the alignment points to
                    alignmentCandidate.tPos,         // ofset into the subsequence
                    alignmentCandidate.blocks[lastBlock].tPos +
                        alignmentCandidate.blocks[lastBlock].length -
                        alignmentCandidate.blocks[0].tPos);
  
          T_AlignmentCandidate refinedAlignment;
  
          //
          // When the parameter bandSize is 0, set the alignment band size
          // to the drift off the diagonal, plus a little more for wiggle
          // room.  When the parameteris nonzero, use that as a fixed band.
          //
          int k;
          if (params.bandSize == 0) {
              k = std::abs(drift) * 1.5;
          } else {
              k = params.bandSize;
          }
          if (params.verbosity > 0) {
              std::cout << "drift: " << drift << " qlen: " << alignmentCandidate.qAlignedSeq.length
                        << " tlen: " << alignmentCandidate.tAlignedSeq.length << " k: " << k
                        << std::endl;
              std::cout << "aligning in " << k << " * " << alignmentCandidate.tAlignedSeq.length
                        << " " << k * alignmentCandidate.tAlignedSeq.length << std::endl;
          }
          if (k < 10) {
              k = 10;
          }
  
          alignmentCandidate.tAlignedSeqPos += alignmentCandidate.tPos;
  
          VectorIndex lastSDPBlock = alignmentCandidate.blocks.size() - 1;
  
          if (alignmentCandidate.blocks.size() > 0) {
              alignmentCandidate.tAlignedSeqLength = (alignmentCandidate.blocks[lastSDPBlock].tPos +
                                                      alignmentCandidate.blocks[lastSDPBlock].length -
                                                      alignmentCandidate.blocks[0].tPos);
          } else {
              alignmentCandidate.tAlignedSeqLength = 0;
          }
  
          alignmentCandidate.tPos = 0;
          alignmentCandidate.qAlignedSeqPos += alignmentCandidate.qPos;
  
          if (alignmentCandidate.blocks.size() > 0) {
              alignmentCandidate.qAlignedSeqLength = (alignmentCandidate.blocks[lastSDPBlock].qPos +
                                                      alignmentCandidate.blocks[lastSDPBlock].length -
                                                      alignmentCandidate.blocks[0].qPos);
          } else {
              alignmentCandidate.qAlignedSeqLength = 0;
          }
          alignmentCandidate.qPos = 0;
  
          alignmentCandidate.blocks.clear();
          alignmentCandidate.tAlignedSeq.Free();
          alignmentCandidate.tAlignedSeq.TakeOwnership(tSeq);
          alignmentCandidate.ReassignQSequence(qSeq);
  
          if (params.verbosity >= 2) {
              std::cout << "refining target: " << std::endl;
              alignmentCandidate.tAlignedSeq.PrintSeq(std::cout);
              std::cout << "refining query: " << std::endl;
              static_cast<DNASequence *>(&alignmentCandidate.qAlignedSeq)->PrintSeq(std::cout);
              std::cout << std::endl;
          }
          PairwiseLocalAlign(qSeq, tSeq, k, params, alignmentCandidate, mappingBuffers, Fit);
      }
  }
  
  template <typename T_RefSequence, typename T_Sequence>
  void RefineAlignments(std::vector<T_Sequence *> &bothQueryStrands, T_RefSequence &genome,
                        std::vector<T_AlignmentCandidate *> &alignmentPtrs, MappingParameters &params,
                        MappingBuffers &mappingBuffers)
  {
  
      UInt i;
      for (i = 0; i < alignmentPtrs.size(); i++) {
          RefineAlignment(bothQueryStrands, genome, *alignmentPtrs[i], params, mappingBuffers);
      }
      //
      // It's possible the alignment references change their order after running
      // the local alignments.  This is made into a parameter rather than resorting
      // every time so that the performance gain by resorting may be measured.
      //
      if (params.sortRefinedAlignments) {
          std::sort(alignmentPtrs.begin(), alignmentPtrs.end(), SortAlignmentPointersByScore());
      }
  }
  
  std::vector<T_AlignmentCandidate *> SelectAlignmentsToPrint(
      std::vector<T_AlignmentCandidate *> alignmentPtrs, MappingParameters &params,
      const int &associatedRandInt)
  {
      if (params.placeRandomly) {
          assert(params.hitPolicy.IsRandombest());
      }
  
      if (alignmentPtrs.size() == 0) {
          return std::vector<T_AlignmentCandidate *>({});
      }
  
      std::sort(alignmentPtrs.begin(), alignmentPtrs.end(), SortAlignmentPointersByScore());
  
      // Apply filter criteria and hit policy.
      // Shallow copy AlignmentCandidate pointers.
      std::vector<T_AlignmentCandidate *> filtered;
      for (auto ptr : alignmentPtrs) {
          if (params.filterCriteria.Satisfy(ptr)) {
              filtered.push_back(ptr);
              if (int(filtered.size()) == params.nBest) break;
          }
      }
  
      return params.hitPolicy.Apply(filtered, false, associatedRandInt);
  }
  
  // The full read is not the subread, and does not have masked off characters.
  void PrintAlignment(T_AlignmentCandidate &alignment, SMRTSequence &fullRead,
                      MappingParameters &params, AlignmentContext &alignmentContext,
                      std::ostream &outFile
  #ifdef USE_PBBAM
                      ,
                      SMRTSequence &subread, PacBio::BAM::IRecordWriter *bamWriterPtr
  #endif
                      )
  {
      try {
          // Before printing alignments, make sure query is always Forward.
          // If query was Reverse, reverse complement both query and target
          // and recomputes all coordinates.
          alignment.MakeQueryForward();
  
          if (params.printFormat == StickPrint) {
              PrintAlignmentStats(alignment, outFile);
              StickPrintAlignment(alignment, (DNASequence &)alignment.qAlignedSeq,
                                  (DNASequence &)alignment.tAlignedSeq, outFile,
                                  alignment.qAlignedSeqPos, alignment.tAlignedSeqPos);
          } else if (params.printFormat == SAM) {
              SAMOutput::PrintAlignment(alignment, fullRead, outFile, alignmentContext,
                                        params.samQVList, params.clipping, params.cigarUseSeqMatch,
                                        params.allowAdjacentIndels);
          } else if (params.printFormat == BAM) {
  #ifdef USE_PBBAM
              BAMOutput::PrintAlignment(alignment, fullRead, subread, *bamWriterPtr, alignmentContext,
                                        params.samQVList, params.clipping, params.cigarUseSeqMatch,
                                        params.allowAdjacentIndels);
  #else
              REQUIRE_PBBAM_ERROR();
  #endif
          } else if (params.printFormat == CompareXML) {
              XMLOutput::Print(alignment, (DNASequence &)alignment.qAlignedSeq,
                               (DNASequence &)alignment.tAlignedSeq, outFile,
                               alignment.qAlignedSeqPos, alignment.tAlignedSeqPos);
          } else if (params.printFormat == Vulgar) {
              PrintAlignmentStats(alignment, outFile);
              VulgarOutput::Print(alignment, outFile);
          } else if (params.printFormat == CompareSequencesParsable) {
              CompareSequencesOutput::Print(alignment, alignment.qAlignedSeq, alignment.tAlignedSeq,
                                            outFile);
          } else if (params.printFormat == Interval) {
              if (alignment.blocks.size() > 0) {
                  IntervalOutput::Print(alignment, outFile);
              }
          } else if (params.printFormat == SummaryPrint) {
              if (alignment.blocks.size() > 0) {
                  SummaryOutput::Print(alignment, outFile);
              }
          }
      } catch (std::ostream::failure f) {
          std::cout << "ERROR writing to output file. The output drive may be full, or you  "
                    << std::endl;
          std::cout << "may not have proper write permissions." << std::endl;
          std::exit(EXIT_FAILURE);
      }
  }
  
  // Print all alignments in std::vector<T_AlignmentCandidate*> alignmentPtrs
  void PrintAlignments(std::vector<T_AlignmentCandidate *> alignmentPtrs, SMRTSequence &read,
                       MappingParameters &params, std::ostream &outFile,
                       AlignmentContext alignmentContext,
  #ifdef USE_PBBAM
                       SMRTSequence &subread, PacBio::BAM::IRecordWriter *bamWriterPtr,
  #endif
                       MappingSemaphores &semaphores)
  {
      if (params.nProc > 1) {
  #ifdef __APPLE__
          sem_wait(semaphores.writer);
  #else
          sem_wait(&semaphores.writer);
  #endif
      }
      for (int i = 0; i < int(alignmentPtrs.size()); i++) {
          T_AlignmentCandidate *aref = alignmentPtrs[i];
  
          if (aref->blocks.size() == 0) {
  
              //
              // If the SDP alignment finds nothing, there will be no
              // blocks.  This may happen if the sdp block size is larger
              // than the anchor size found with the suffix array.  When no
              // blocks are found there is no alignment, so zero-out the
              // score and continue.
              //
              aref->score = 0;
              if (params.verbosity > 0) {
                  std::cout << "Zero blocks found for " << aref->qName << " " << aref->qAlignedSeqPos
                            << " " << aref->tAlignedSeqPos << std::endl;
              }
              continue;
          }
  
          //
          // Configure some of the alignment context before printing.
          //
          if (i > 0 and params.placeRandomly == false) {
              alignmentContext.isPrimary = false;
          } else {
              alignmentContext.isPrimary = true;
          }
  
          if (params.printSAM or params.printBAM) {
              DistanceMatrixScoreFunction<DNASequence, FASTASequence> editdistScoreFn(
                  EditDistanceMatrix, 1, 1);
              T_AlignmentCandidate &alignment = *alignmentPtrs[i];
              alignmentContext.editDist = ComputeAlignmentScore(
                  alignment, alignment.qAlignedSeq, alignment.tAlignedSeq, editdistScoreFn);
          }
  
          PrintAlignment(*alignmentPtrs[i], read, params, alignmentContext, outFile
  #ifdef USE_PBBAM
                         ,
                         subread, bamWriterPtr
  #endif
                         );
      }
  
      if (params.nProc > 1) {
  #ifdef __APPLE__
          sem_post(semaphores.writer);
  #else
          sem_post(&semaphores.writer);
  #endif
      }
  }
  
  void PrintAlignmentPtrs(std::vector<T_AlignmentCandidate *> &alignmentPtrs, std::ostream &out)
  {
      for (int alignmentIndex = 0; alignmentIndex < int(alignmentPtrs.size()); alignmentIndex++) {
          out << "[" << alignmentIndex << "/" << alignmentPtrs.size() << "]" << std::endl;
          T_AlignmentCandidate *alignment = alignmentPtrs[alignmentIndex];
          alignment->Print(out);
      }
      out << std::endl;
  }
  
  void PrintUnaligned(const SMRTSequence &unalignedRead, std::ostream &unalignedFilePtr,
                      const bool noPrintUnalignedSeqs)
  {
      if (noPrintUnalignedSeqs) {
          std::string s = unalignedRead.GetTitle();
          SMRTTitle st(s);
          if (st.isSMRTTitle)
              unalignedFilePtr << st.ToString() << std::endl;
          else
              //size_t pos = s.rfind("/");
              //if (pos != string::npos)
              //    unalignedFilePtr << s.substr(0, pos) << std::endl;
              //else
              unalignedFilePtr << s << std::endl;
      } else
          unalignedRead.PrintSeq(unalignedFilePtr);
  }
  
  // Print all alignments for subreads in allReadAlignments.
  // Input:
  //   allReadAlignments - contains a set of subreads, each of which
  //                       is associated with a group of alignments.
  //   alignmentContext  - an alignment context of each subread used
  //                       for printing in SAM format.
  //   params            - mapping parameters.
  // Output:
  //   outFilePtr        - where to print alignments for subreads.
  //   unalignedFilePtr  - where to print sequences for unaligned subreads.
  void PrintAllReadAlignments(ReadAlignments &allReadAlignments, AlignmentContext &alignmentContext,
                              std::ostream &outFilePtr, std::ostream &unalignedFilePtr,
                              MappingParameters &params, std::vector<SMRTSequence> &subreads,
  #ifdef USE_PBBAM
                              PacBio::BAM::IRecordWriter *bamWriterPtr,
  #endif
                              MappingSemaphores &semaphores)
  {
      int subreadIndex;
      int nAlignedSubreads = allReadAlignments.GetNAlignedSeq();
  
      //
      // Initialize the alignemnt context with information applicable to SAM output.
      //
      alignmentContext.alignMode = allReadAlignments.alignMode;
      for (subreadIndex = 0; subreadIndex < nAlignedSubreads; subreadIndex++) {
          if (allReadAlignments.subreadAlignments[subreadIndex].size() > 0) {
              alignmentContext.numProperlyAlignedSubreads++;
          }
      }
  
      if (alignmentContext.numProperlyAlignedSubreads ==
          int(allReadAlignments.subreadAlignments.size())) {
          alignmentContext.allSubreadsProperlyAligned = true;
      }
      alignmentContext.nSubreads = nAlignedSubreads;
  
      for (subreadIndex = 0; subreadIndex < nAlignedSubreads; subreadIndex++) {
          alignmentContext.subreadIndex = subreadIndex;
          if (subreadIndex < nAlignedSubreads - 1 and
              allReadAlignments.subreadAlignments[subreadIndex + 1].size() > 0) {
              alignmentContext.nextSubreadPos =
                  allReadAlignments.subreadAlignments[subreadIndex + 1][0]->QAlignStart();
              alignmentContext.nextSubreadDir =
                  allReadAlignments.subreadAlignments[subreadIndex + 1][0]->qStrand;
              alignmentContext.rNext =
                  allReadAlignments.subreadAlignments[subreadIndex + 1][0]->tName;
              alignmentContext.hasNextSubreadPos = true;
          } else {
              alignmentContext.nextSubreadPos = 0;
              alignmentContext.nextSubreadDir = 0;
              alignmentContext.rNext = "";
              alignmentContext.hasNextSubreadPos = false;
          }
          SMRTSequence *sourceSubread = &(allReadAlignments.subreads[subreadIndex]);
          if (subreads.size() == allReadAlignments.subreads.size()) {
              sourceSubread = &subreads[subreadIndex];
          }
          if (allReadAlignments.subreadAlignments[subreadIndex].size() > 0) {
              PrintAlignments(allReadAlignments.subreadAlignments[subreadIndex],
                              allReadAlignments.subreads[subreadIndex],
                              // for these alignments
                              params, outFilePtr,  //*mapData->outFilePtr,
                              alignmentContext,
  #ifdef USE_PBBAM
                              *sourceSubread, bamWriterPtr,
  #endif
                              semaphores);
          } else {
              //
              // Print the unaligned sequences.
              //
              if (params.printUnaligned == true) {
                  if (params.nProc == 1) {
                      PrintUnaligned(*sourceSubread, unalignedFilePtr, params.noPrintUnalignedSeqs);
                  } else {
  #ifdef __APPLE__
                      sem_wait(semaphores.unaligned);
  #else
                      sem_wait(&semaphores.unaligned);
  #endif
                      PrintUnaligned(*sourceSubread,  //subreads[subreadIndex],
                                     unalignedFilePtr, params.noPrintUnalignedSeqs);
  #ifdef __APPLE__
                      sem_post(semaphores.unaligned);
  #else
                      sem_post(&semaphores.unaligned);
  #endif
                  }  // End of nproc > 1.
              }      // End of printing  unaligned sequences.
          }          // End of finding no alignments for the subread with subreadIndex.
      }              // End of printing and processing alignmentContext for each subread.
  }
  