/*****************************************************************************
#   Copyright (C) 1994-2008 by David Gordon.
#   All rights reserved.                           
#                                                                           
#   This software is part of a beta-test version of the Consed/Autofinish
#   package.  It should not be redistributed or
#   used for any commercial purpose, including commercially funded
#   sequencing, without written permission from the author and the
#   University of Washington.
#   
#   This software is provided ``AS IS'' and any express or implied
#   warranties, including, but not limited to, the implied warranties of
#   merchantability and fitness for a particular purpose, are disclaimed.
#   In no event shall the authors or the University of Washington be
#   liable for any direct, indirect, incidental, special, exemplary, or
#   consequential damages (including, but not limited to, procurement of
#   substitute goods or services; loss of use, data, or profits; or
#   business interruption) however caused and on any theory of liability,
#   whether in contract, strict liability, or tort (including negligence
#   or otherwise) arising in any way out of the use of this software, even
#   if advised of the possibility of such damage.
#
#   Building Consed from source is error prone and not simple which is
#   why I provide executables.  Due to time limitations I cannot
#   provide any assistance in building Consed.  Even if you do not
#   modify the source, you may introduce errors due to using a
#   different version of the compiler, a different version of motif,
#   different versions of other libraries than I used, etc.  For this
#   reason, if you discover Consed bugs, I can only offer help with
#   those bugs if you first reproduce those bugs with an executable
#   provided by me--not an executable you have built.
# 
#   Modifying Consed is also difficult.  Although Consed is modular,
#   some modules are used by many other modules.  Thus making a change
#   in one place can have unforeseen effects on many other features.
#   It may takes months for you to notice these other side-effects
#   which may not seen connected at all.  It is not feasable for me to
#   provide help with modifying Consed sources because of the
#   potentially huge amount of time involved.
#
#*****************************************************************************/
#include    "contig.h"
#include    "primerType.h"
#include    "experiment.h"
#include    "locatedFragment.h"
#include    "mbtValVector.h"
#include    "consedParameters.h"
#include    "min.h"
#include    "max.h"
#include    "subcloneTTemplate.h"
#include    "consed.h"
#include    "fileDefines.h"



static int nMaxSubcloneCoverage = 10000;  // some big number


// This file contains the parts of autofinish for covering single 
// subclone regions.


void Contig :: coverSingleSubcloneRegions() {

   // now this is already done in autoFinish.cpp
   // so that evalExpCluster can take advantage of where the
   // single subclone regions are.
   //   autoFinishFindSingleSubcloneBases();


   if ( nCountSingleSubcloneBases() == 0 ) {
      fprintf( pAO, "No single subclone regions\n" );
      return;
   }


   autoFinishFindSingleSubcloneRegionsFromSuggestedExperiments();


   autoFinishSetCumulativeSingleSubcloneArray();


   autoFinishCreateCandidateExperimentsToCoverSingleSubcloneRegions();

   autoFinishChooseSingleSubcloneExperiments();
}




bool Contig :: bGetLikelyAlignedRegionOnContig( experiment* pExp,
                                                int& nConsPosLeftUnpadded,
                                                int& nConsPosRightUnpadded ) {

   if ( pExp->nReadType_ == nWalk )
      return( bGetLikelyAlignedRegionOnContig2( 
               pExp->nFirstBaseOfReadUnpadded_,
               pExp->bTopStrandNotBottomStrand_,
               nConsPosLeftUnpadded,
               nConsPosRightUnpadded ) );
   else {
      // clearly vector bases are not going to fix single subclone bases
      // if the read is a universal primer read, 
      // nUnpaddedLeft_ and nUnpaddedRight_ have already excluded
      // vector bases.  The average aligned portion of reads
      // may not do this, hence the intersection.

      if ( bGetLikelyAlignedRegionOnContig2(                
               pExp->nFirstBaseOfReadUnpadded_,
               pExp->bTopStrandNotBottomStrand_,
               nConsPosLeftUnpadded,
               nConsPosRightUnpadded ) ) {
         // now must exclude vector bases

         return( bIntersect( nConsPosLeftUnpadded,
                          nConsPosRightUnpadded,
                          pExp->nUnpaddedLeft_,
                          pExp->nUnpaddedRight_,
                          nConsPosLeftUnpadded,
                          nConsPosRightUnpadded ) );
      }
      else
         return( false );
   }
}



// I think this is only valid for walks since it doesn't
// consider vector bases at the beginning
bool Contig :: bGetLikelyAlignedRegionOnContig2( 
                        const int nFirstBaseOfReadUnpadded,
                        const bool bTopStrandNotBottomStrand,
                        int& nConsPosLeftUnpadded,
                        int& nConsPosRightUnpadded ) {

   if ( bTopStrandNotBottomStrand ) {
      nConsPosLeftUnpadded = nFirstBaseOfReadUnpadded +
         ConsEd::pGetAssembly()->nModelReadAlignedSegmentStart_ - 1;
      nConsPosRightUnpadded = nFirstBaseOfReadUnpadded +
         ConsEd::pGetAssembly()->nModelReadAlignedSegmentEnd_ - 1;
   }
   else {
      nConsPosRightUnpadded = nFirstBaseOfReadUnpadded -
         ConsEd::pGetAssembly()->nModelReadAlignedSegmentStart_ + 1;
      nConsPosLeftUnpadded = nFirstBaseOfReadUnpadded -
         ConsEd::pGetAssembly()->nModelReadAlignedSegmentEnd_ + 1;
   }

   if ( bIntersect( nConsPosLeftUnpadded,
                    nConsPosRightUnpadded,
                    nGetUnpaddedStartIndex(),
                    nGetUnpaddedEndIndex(),
                    nConsPosLeftUnpadded,
                    nConsPosRightUnpadded ) ) {

      return( true );
   }
   else
      return( false );
}




void Contig :: sortTemplatesByHowManySingleSubcloneBasesEachFixes( 
                              primerType* pPrimer,
                              int nSingleSubcloneBasesFixedByThisTemplate[ nNUMBER_OF_TEMPLATES] ) {

   int nTemplate;
   for( nTemplate = 0; nTemplate < ( nNUMBER_OF_TEMPLATES - 1); ++nTemplate ) {
      for( int nTemplate2 = nTemplate + 1; nTemplate2 < nNUMBER_OF_TEMPLATES;
           ++nTemplate2 ) {

         if ( nSingleSubcloneBasesFixedByThisTemplate[ nTemplate2 ] >   
              nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] ) {
            int nTemp = nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ];
            subcloneTTemplate* pSubTemp = pPrimer->pSubcloneTemplate_[ nTemplate ];
            nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] =
               nSingleSubcloneBasesFixedByThisTemplate[ nTemplate2 ];
            
            pPrimer->pSubcloneTemplate_[ nTemplate ] =
               pPrimer->pSubcloneTemplate_[ nTemplate2 ];

            
            nSingleSubcloneBasesFixedByThisTemplate[ nTemplate2 ] =
               nTemp;

            pPrimer->pSubcloneTemplate_[ nTemplate2 ] = pSubTemp;
         }
      }
   }

   // check above sort


   for( nTemplate = 1; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
      assert( nSingleSubcloneBasesFixedByThisTemplate[ nTemplate - 1 ]
              >=
              nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] );
   }






}



void Contig :: eliminateTemplatesThatDoNotSolveAnySingleSubcloneBases(
           primerType* pPrimer,
           int nSingleSubcloneBasesFixedByThisTemplate[ nNUMBER_OF_TEMPLATES] ) {

   for( int nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
      if ( pPrimer->pSubcloneTemplate_[ nTemplate ] &&
           nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] == 0 ) {

         pPrimer->pSubcloneTemplate_[ nTemplate ] = NULL;
      }
   }
}






void Contig :: autoFinishSetSingleSubcloneArrays( const bool bObserveDoNotFinishTags ) {
  
   int nUnpaddedLength = nGetUnpaddedSeqLength();

   pWhichSingleSubcloneIsCoveringEachBase_ = 
      new mbtPtrVector<LocatedFragment>( nUnpaddedLength,  nGetUnpaddedStartIndex() );
   
   pNumberOfSubclonesCoveringEachBase_ = 
      new mbtValVector<int>( nUnpaddedLength, nGetUnpaddedStartIndex(), 0 );


   pNumberOfSubclonesCoveringEachBaseWithoutSuggestedReads_ =
      new mbtValVector<int>( nUnpaddedLength, nGetUnpaddedStartIndex(), 0 );

   autoFinishFindSingleSubcloneRegionsFromExistingReads();

   if ( bObserveDoNotFinishTags )
      autoFinishExamineDoNotFinishTagsForSingleSubcloneRegions();

   saveCopyOfSingleSubcloneArray();
}



void Contig :: saveCopyOfSingleSubcloneArray() {
   for( int nUnpadded = pNumberOfSubclonesCoveringEachBase_->nGetStartIndex();
        nUnpadded <= pNumberOfSubclonesCoveringEachBase_->nGetEndIndex();
        ++nUnpadded ) {
      (*pNumberOfSubclonesCoveringEachBaseWithoutSuggestedReads_)[ nUnpadded ] =
         (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ];
         
   }
}




void Contig :: autoFinishSetCumulativeSingleSubcloneArray() {

   // set cumulative subclone array.  This will just be used 
   // once--to decide which experiment candidates to create

   int nUnpaddedLength = nGetUnpaddedSeqLength();


   // the offset is nGetUnpaddedStartIndex() - 1 instead of
   // nGetUnpaddedStartIndex()
   // This allows us to take the part of the experiment that is on the
   // consensus (nLeft, nRight), and then just use nLeft - 1 and nRight
   // as subscripts to the array below without worrying about subscript 
   // errors.  If we had used nGetUnpaddedStartIndex() as the offset,
   // then if the experiment hung out over the left edge of the consensus,
   // nLeft - 1 would give a subscript error.  Furthermore, we need to
   // include the number of single subclone bases at left-most base of the
   // consensus.

   pCumNumberOfSingleSubcloneBases_ = new mbtValVector<int>( nUnpaddedLength + 1, 
                                                          nGetUnpaddedStartIndex() - 1,
                                                          0 );
   
   // set the element at nGetUnpaddedStartIndex() - 1 to 0 so that 
   // every other element is indeed cumulative from this base

   (*pCumNumberOfSingleSubcloneBases_)[ 
            pCumNumberOfSingleSubcloneBases_->nGetStartIndex()
                                       ]  = 0;

   
   for( int nUnpaddedConsPos = pCumNumberOfSingleSubcloneBases_->nGetStartIndex() + 1;
        nUnpaddedConsPos <=  pCumNumberOfSingleSubcloneBases_->nGetEndIndex();
        ++nUnpaddedConsPos ) {

      int nSingleSubcloneOnThisBase =
         ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpaddedConsPos ] < 2 ) ?
         1 : 0;

      (*pCumNumberOfSingleSubcloneBases_)[ nUnpaddedConsPos ] =
         (*pCumNumberOfSingleSubcloneBases_)[ nUnpaddedConsPos - 1 ]
         + 
         nSingleSubcloneOnThisBase;
   }

}




void Contig :: autoFinishCreateCandidateExperimentsToCoverSingleSubcloneRegions() {


   // just throw away previous list of experiment candidates
   // (If consed were to stay up, rather than terminate, this would
   // be a memory leak.  In that case, we would need to find the
   // experiment candidates that are not also on the experiments list
   // and delete them.)

   pCandidateExperiments_ = new mbtPtrOrderedVector<experiment>;
   
   autoFinishCreateTopStrandCandidateExperimentsForSingleSubcloneRegions();

   autoFinishCreateBottomStrandCandidateExperimentsForSingleSubcloneRegions();


}


void Contig :: autoFinishCreateTopStrandCandidateExperimentsForSingleSubcloneRegions() {

   pCP->setParametersForSequencingPrimersNotPCRPrimers( true );
   

   int nFirstExperimentPrimer3PrimeEnd = nGetUnpaddedStartIndex() +  
      consedParameters::pGetConsedParameters()->nPrimersMaximumLengthOfAPrimerToUse_ 
      - 1;

   int nLastExperimentPrimer3PrimeEnd = nGetUnpaddedEndIndex();


   for( int nPrimer3PrimeEnd = nFirstExperimentPrimer3PrimeEnd;
        nPrimer3PrimeEnd <= nLastExperimentPrimer3PrimeEnd;
        ++nPrimer3PrimeEnd ) {

      int nFirstBaseOfReadUnpadded = nPrimer3PrimeEnd + 1;

      int nLikelyAlignedLeftEndUnpadded;
      int nLikelyAlignedRightEndUnpadded;

      // this likely aligned portion of the read may not be completely
      // on the contig.  In fact, it might not be on at all.

      // case in which the high quality portion is completely off
      // the contig so no single subclone bases are fixed
      if ( !bGetLikelyAlignedRegionOnContig2( nFirstBaseOfReadUnpadded,
                                            true, // top strand not bottom
                                            nLikelyAlignedLeftEndUnpadded,
                                            nLikelyAlignedRightEndUnpadded ) )
         continue;

      

      // This will give a rough idea of how many single subclone
      // bases will be fixed by this experiment
      // (A precise number must take into account which template
      // is already covering each base.)
         
      int nMaxNumberOfSingleSubcloneBasesFixedByThisRead = 
         (*pCumNumberOfSingleSubcloneBases_)[ nLikelyAlignedRightEndUnpadded ]
         -
         (*pCumNumberOfSingleSubcloneBases_)[ nLikelyAlignedLeftEndUnpadded  - 1 ];

      if ( nMaxNumberOfSingleSubcloneBasesFixedByThisRead 
           <
           consedParameters::pGetConsedParameters()->nAutoFinishMinNumberOfSingleSubcloneBasesFixedByAnExp_ )
         continue;

      // check that the primer is high enough quality
      if ( !bCustomPrimerQualityTest(
                                     nFirstBaseOfReadUnpadded,
                                     true ) ) continue;

         
      if ( consedParameters::pGetConsedParameters()->bAutoFinishAllowWholeCloneReads_ ) {
         maybeCreateExperimentWithCustomPrimer(
                                            TOP_STRAND_TERMINATOR_READS,
                                            true, // bTopStrandNotBottonStrand,
                                            nFirstBaseOfReadUnpadded,
                                            true, // bCloneTemplateNotSubcloneTemplate
                                            true ); // for single subclone regions
      }



      if ( pCP->bAutoFinishAllowCustomPrimerSubcloneReads_ ) {
         maybeCreateExperimentWithCustomPrimer(
                                            TOP_STRAND_TERMINATOR_READS,
                                            true, // bTopStrandNotBottomStrand,
                                            nFirstBaseOfReadUnpadded,
                                            false, // bCloneTemplateNotSubcloneTemplate
                                            true ); // for single subclone regions
      }

   }
}




void Contig :: autoFinishCreateBottomStrandCandidateExperimentsForSingleSubcloneRegions() {


   pCP->setParametersForSequencingPrimersNotPCRPrimers( true );


   int nFirstExperimentPrimer3PrimeEnd = nGetUnpaddedEndIndex() - 
      consedParameters::pGetConsedParameters()->nPrimersMaximumLengthOfAPrimerToUse_ 
      + 1;

   int nLastExperimentPrimer3PrimeEnd = nGetUnpaddedStartIndex();

   for( int nPrimer3PrimeEnd = nFirstExperimentPrimer3PrimeEnd;
        nPrimer3PrimeEnd >= nLastExperimentPrimer3PrimeEnd;
        --nPrimer3PrimeEnd ) {

      int nFirstBaseOfReadUnpadded = nPrimer3PrimeEnd - 1;

      int nLikelyAlignedLeftEndUnpadded;
      int nLikelyAlignedRightEndUnpadded;

      if ( !bGetLikelyAlignedRegionOnContig2( nFirstBaseOfReadUnpadded,
                                            false, // top strand not bottom
                                            nLikelyAlignedLeftEndUnpadded,
                                            nLikelyAlignedRightEndUnpadded ) )
         continue;
      

      // This will give a rough idea of how many single subclone
      // bases will be fixed by this experiment.  (A precise number
      // must take into account which template is already covering 
      // each base.)

      int nMaxNumberOfSingleSubcloneBasesFixedByThisRead =
         (*pCumNumberOfSingleSubcloneBases_)[ nLikelyAlignedRightEndUnpadded ]
         -
         (*pCumNumberOfSingleSubcloneBases_)[ nLikelyAlignedLeftEndUnpadded - 1 ];

      if ( nMaxNumberOfSingleSubcloneBasesFixedByThisRead 
           <
           consedParameters::pGetConsedParameters()->nAutoFinishMinNumberOfSingleSubcloneBasesFixedByAnExp_ )
         continue;

      // check that the primer is high enough quality to be made
      if ( !bCustomPrimerQualityTest(
                                     nFirstBaseOfReadUnpadded,
                                     false ) // bTopStrandNotBottomStrand
           ) continue;

      if ( consedParameters::pGetConsedParameters()->bAutoFinishAllowWholeCloneReads_ ) {
         maybeCreateExperimentWithCustomPrimer(
                                            BOTTOM_STRAND_TERMINATOR_READS,
                                            false, // bTopStrandNotBottomStrand,
                                            nFirstBaseOfReadUnpadded,
                                            true, // bCloneTemplateNotSubcloneTemplate
                                            true ); // for single subclone regions
      }

      if ( pCP->bAutoFinishAllowCustomPrimerSubcloneReads_ ) {
         maybeCreateExperimentWithCustomPrimer(
                                            BOTTOM_STRAND_TERMINATOR_READS,
                                            false, // bTopStrandNotBottomStrand,
                                            nFirstBaseOfReadUnpadded,
                                            false, // bCloneTemplateNotSubcloneTemplate
                                            true ); // for single subclone regions
      }
   } //   for( int nPrimer3PrimeEnd = nFirstExperiment3PrimeEnd;
}







void Contig :: autoFinishFindSingleSubcloneRegionsFromExistingReads() {


   for(  int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
      LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );

      if ( pLocFrag->bIsWholeReadUnaligned() )
         continue;

      int nUnpaddedConsPosStart = nUnpaddedIndex( pLocFrag->nGetAlignClipStart() );
      int nUnpaddedConsPosEnd = nUnpaddedIndex( pLocFrag->nGetAlignClipEnd() );

      nUnpaddedConsPosStart = MAX( nUnpaddedConsPosStart, nGetUnpaddedStartIndex() );
      nUnpaddedConsPosEnd = MIN( nUnpaddedConsPosEnd, nGetUnpaddedEndIndex() );


      // this takes care of reads that are totally off the contig
      if ( nUnpaddedConsPosEnd < nUnpaddedConsPosStart ) continue;

      RWCString soSubcloneName = pLocFrag->soGetTemplateName();

      for( int nConsPosUnpadded = nUnpaddedConsPosStart; 
           nConsPosUnpadded <= nUnpaddedConsPosEnd; 
           ++nConsPosUnpadded ) {

         // if *any* of the reads covering
         // this region are whole clone reads, then I should
         // consider the region to be covered by more than one
         // subclone since, if there is a deletion in the whole clone,
         // additional reads through the region won't help

         if ( pLocFrag->bIsWholeCloneTemplateNotSubcloneTemplate() ) {
            (*pNumberOfSubclonesCoveringEachBase_)[ nConsPosUnpadded ] = 
               nMaxSubcloneCoverage;
            // in case someone wants more than double subclone coverage
         }
         else {
            if ( !(*pNumberOfSubclonesCoveringEachBase_)[ nConsPosUnpadded ] ) {
               // there is no read so far covering this position
               (*pNumberOfSubclonesCoveringEachBase_)[ nConsPosUnpadded ] = 1;
               (*pWhichSingleSubcloneIsCoveringEachBase_)[ nConsPosUnpadded ] = pLocFrag;
            }
            else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nConsPosUnpadded ] == 1 ) {
               // there is some other read covering this base or 
               // this base *will* be covered by some read that
               // will be suggested by autofinish

               // Let's see if it is the same subclone.

               LocatedFragment* pLocFrag2 = (*pWhichSingleSubcloneIsCoveringEachBase_)[ nConsPosUnpadded ];

               if ( pLocFrag2->soGetTemplateName() != soSubcloneName ) {
                  (*pNumberOfSubclonesCoveringEachBase_)[ nConsPosUnpadded ] = 2;
               }
            } // else if ( (*pNumberOfSubclonesCoveringEachBase_)...

            // if there are already 2 subclones covering this region,
            // don't bother checking for more
         } // if ( pLocFrag->bIsWholeCloneTemplate() ) ... else
      } //  for( int nConsPosUnpadded = nConsPosStart...
   }  // for(  int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead )

   
   // Now we have found all the regions covered by existing reads.
}



void Contig :: autoFinishFindSingleSubcloneRegionsFromSuggestedExperiments() {


   // Next, add to that the regions that will be covered by the
   // potential high quality portion of suggested experiments

   // Suggested experiments are of the following 3 types:
   //    1) universal primer reads, which won't help at all
   //    2) custom primer whole clone reads, which will solve every base
   //           regardless of what else covers those bases
   //    3) custom primer subclone reads, which may or may not solve
   //           a base, depending on what other subclone covers that base


   for( int nExp = 0; nExp < pChosenExperiments_->length(); ++nExp ) {
      experiment* pExp = (*pChosenExperiments_)[ nExp ];

      //      if ( !pExp->bCustomPrimerNotUniversalPrimer_ ) continue;
      if ( pExp->bIsResequencingRead() ) continue;

      int nConsPosLeftUnpadded;
      int nConsPosRightUnpadded;


      // handle case in which the high quality portion of the
      // read is completely off the contig
      if (!bGetLikelyAlignedRegionOnContig( pExp, 
                                          nConsPosLeftUnpadded,
                                          nConsPosRightUnpadded ) ) continue;


      if ( pExp->bIsUniversalPrimerRead() ) {

         subcloneTTemplate* pSub = pExp->pSubcloneTemplate_;
         LocatedFragment* pLocFragOfSub = pSub->aReads_[0];
         
         for( int nUnpadded = nConsPosLeftUnpadded; 
              nUnpadded <= nConsPosRightUnpadded; ++nUnpadded ) {
            if ( !(*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] ) {
               // there is no read so far covering this position
               (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 1;
               (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ] =
                  pLocFragOfSub;
            }
            else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 
                      1 ) {
               // there is some other read covering this base or
               // this base is covered by some other read that has been
               // suggested by Autofinish in this round.
               // is it the same subclone?
               
               LocatedFragment* pLocFrag2 = 
                  (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ];
               if ( pLocFrag2->soGetTemplateName() != pSub->soTemplateName_ ) {
                  (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 2;
               }
            } // if  ( !(*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] ) else
         } //  for( int nUnpadded  ...
      } //   if ( pExp->bIsUniveralPrimerRead() ) {
      else if ( pExp->bCloneTemplateNotSubcloneTemplate_ ) {
         for( int nUnpadded = nConsPosLeftUnpadded; 
              nUnpadded <= nConsPosRightUnpadded; ++nUnpadded ) {
            (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 
               nMaxSubcloneCoverage;
         }
      }
      else {
         // subclone template
         // Let's find which subclone it is from

         primerType* pPrimer = pExp->pCustomPrimer_;
         for( int nTemplate = 0; ( nTemplate < nNUMBER_OF_TEMPLATES ) &&
                 ( nTemplate < 
                 consedParameters::pGetConsedParameters()->nAutoFinishHowManyTemplatesYouIntendToUseForCustomPrimerSubcloneReactions_ ) &&
                 ( pPrimer->pSubcloneTemplate_[ nTemplate ] != NULL );
              ++nTemplate ) {
            subcloneTTemplate* pSubcloneTemplate = pPrimer->pSubcloneTemplate_[ nTemplate ];

            RWCString soSubcloneName = pSubcloneTemplate->soTemplateName_;

	    // fix bug found by Darren Grafham, Sanger:
	    //  ssss TTTTTTTTTTTT  where sss is single subclone and TTT is
	    // a template.  In making a bottom strand read, Autofinish used
	    // to ignore the left end of the TTT and just assume that a read
	    // length of bases would be made, thus fixing this single
	    // subclone region
	 
	    int nUnpaddedLeftOnTemplate;
	    int nUnpaddedRightOnTemplate;

	    if (!bIntersect( nConsPosLeftUnpadded,
			    nConsPosRightUnpadded,
			  pSubcloneTemplate->nUnpaddedStart_,
			  pSubcloneTemplate->nUnpaddedEnd_,
			  nUnpaddedLeftOnTemplate,
			  nUnpaddedRightOnTemplate ) ) continue;
         
	    for( int nUnpadded = nUnpaddedLeftOnTemplate; 
		nUnpadded <= nUnpaddedRightOnTemplate;
		++nUnpadded ) {

               if ( !(*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] ) {
                  // there is no read so far covering this position
                  (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 1;
                  (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ] = 
                     pSubcloneTemplate->aReads_[0];
                  // I just picked an arbitrary read from the subclone
                  // template as a stand-in for it
               }
               else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 1 ) {
                  // there is some other read covering this base
                  // or this base is covered by some read that has
                  // been suggested by autofinish

                  // Let's see if it is the same subclone

                  LocatedFragment* pLocFrag2 = (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ];
                  if ( pLocFrag2->soGetTemplateName() != soSubcloneName ) {
                     (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 2;
                  }
               }
               
               // ignore case in which there are already 2 or more
               // subclones covering this base
               
            }// for( nUnpadded ...
         } //  for( int nTemplate = 0; ( nTemplate < nNUMBER_OF_TEMPLATES ) &&
      } //  if ( pExp->bCloneTemplateNotSubcloneTemplate_ )  else
   } // for( int nExp ...
}



void Contig :: autoFinishChooseSingleSubcloneExperiments() {

   // should continue to find experiments until there are 
   // absolutely no single subclone regions left in this contig

   bool bMoreSingleSubcloneRegions = true;

   
   while( bMoreSingleSubcloneRegions ) {
      experiment* pBestExp = pGetBestSingleSubcloneExperiment();
      if ( !pBestExp ) {
         bMoreSingleSubcloneRegions = false;
      }
      else {
         pBestExp->nPurpose_ = nCoverSingleSubcloneBases;
         findWindowOfMostSingleSubcloneBasesFixed( pBestExp );
         pBestExp->chooseExperiment();
         pBestExp->print();
         pChosenExperiments_->insert( pBestExp );
         ConsEd::pGetAssembly()->pAllChosenExperiments_->insert( pBestExp );
         adjustSingleSubcloneRegionsInRemainingExperiments( pBestExp );
         // the reason this should be called is so that, when we are done
         // with single subclone regions, we can find any
         // remaining weak regions and print them out.

         adjustConsensusDistributionsAndArrays( pBestExp );
      }
   }

   int nRemainingSingleSubcloneBases = nCountSingleSubcloneBases();
   fprintf( pAO, "Number of remaining single subclone bases: %d\n", 
           nCountSingleSubcloneBases() );

   if ( nRemainingSingleSubcloneBases > 0 ) {
      fprintf( pAO, "Could find no experiment to resolve these remaining regions {\n" );
      autoFinishPrintSingleSubcloneRegions();
      fprintf( pAO, "}\n" );
   }
}


void Contig :: adjustSingleSubcloneRegionsInRemainingExperiments( 
                                  experiment* pChosenExperiment ) {

   autoFinishAdjustSingleSubcloneArrays( pChosenExperiment );

   for( int nExp = 0; nExp < pCandidateExperiments_->length(); ++nExp ) {
      experiment* pExpToBeAdjusted = (*pCandidateExperiments_)[ nExp ];

      // universal primer reads will almost never help
      if ( pExpToBeAdjusted->nReadType_ != nWalk ) continue;
      // if ( !pExpToBeAdjusted->bCustomPrimerNotUniversalPrimer_ ) continue;

      if ( pExpToBeAdjusted->nState_ != nAvailableExperiment ) continue;
      
      if ( ! pChosenExperiment->bOverlaps( pExpToBeAdjusted ) ) continue;
      
      if ( pExpToBeAdjusted == pChosenExperiment ) continue;

      // if reached here, the pExpToBeAdjusted overlaps the experiment
      // just chosen.  Recompute the number of single subclone bases
      // it will fix.  

      // If the pExpToBeAdjusted already doesn't fix any single subclone bases,
      // don't bother recomputing anything
      if ( pExpToBeAdjusted->nSingleSubcloneBasesFixed_ == 0 ) continue;

      adjustSingleSubcloneBasesFixedByExperiment( pExpToBeAdjusted );
   }
}
      

void Contig :: adjustSingleSubcloneBasesFixedByExperiment( 
                                    experiment* pExpToBeAdjusted ) {
   
   int nConsPosLeftUnpadded;
   int nConsPosRightUnpadded;

   if ( !bGetLikelyAlignedRegionOnContig( pExpToBeAdjusted,
                                       nConsPosLeftUnpadded,
                                       nConsPosRightUnpadded ) )
      return;

   if ( pExpToBeAdjusted->bCloneTemplateNotSubcloneTemplate_ ) {
      // every single subclone base will be fixed by a whole clone 
      // experiment

      adjustSingleSubcloneBasesFixedByWholeCloneExperiment( pExpToBeAdjusted,
                                                            nConsPosLeftUnpadded,
                                                            nConsPosRightUnpadded );
   } // if ( pExpToBeAdjusted->bCloneTemplateNotSubcloneTemplate_ ) {
   else {

      // Case of subclone template for primer
      adjustSingleSubcloneBasesFixedBySubcloneExperiment( pExpToBeAdjusted,
                                                          nConsPosLeftUnpadded,
                                                          nConsPosRightUnpadded );
   }
}


void Contig :: adjustSingleSubcloneBasesFixedBySubcloneExperiment(
                                          experiment* pExpToBeAdjusted,
                                          const int nConsPosLeftUnpadded,
                                          const int nConsPosRightUnpadded ) {


   primerType* pPrimer = pExpToBeAdjusted->pCustomPrimer_;
   int nSingleSubcloneBasesFixedByThisTemplate[ nNUMBER_OF_TEMPLATES ];

   int nTemplate;
   for( nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
      nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] = 0;
   }

   for( nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES;
        ++nTemplate ) {

      // There is a problem here--each different template may fix a
      // different number of errors.  In effect, each template is a
      // different experiment.  More and more I am thinking that
      // single subclone coverage should be done in a separate
      // pass.  The way I would do it differently is that I would
      // split out the different templates.  I am concerned about
      // the case in which the top most quality templates do no
      // good for single subclone coverage.  Thus it is a waste of
      // money to do them.  Perhaps we could handle this as a
      // special case--if a template does no good, it isn't done.
      // But that would still leave the possibility of doing some
      // templates that only fix a few single subclone regions due
      // to some template of lower quality that fixes a great
      // number of them.  Perhaps for single subclone regions I
      // could split the experiment list into a much longer list
      // (up to 5x as long) in which each experiment only has 1
      // template.  This has the disadvantage in that several
      // templates may each work to solve a single subclone region
      // and you might want to pick more than one in order to be
      // sure of having one that works.  Alternatively, I could
      // score the toal experiment based on how well the best
      // template works.  This has the problem that the group might
      // not do that experiment, especially if it is near the end
      // of the list.  How about if I score the experiment based on
      // the best (or all) templates within the list of templates
      // that are in the top range of templates that the user
      // specified as making?  What about if the last template is
      // the one that will be most useful in solving a single
      // subclone region?  Perhaps the best thing is to asterisk
      // the template that is the one that is best.  Alternatively,
      // we could keep track of a score for each of the templates.
      // Each experiment could have up to 5 different counts of
      // single subclone bases fixed.  And thus there should be up
      // to 5 different single subclone bases per cost.  There are
      // competing objectives with picking templates--the one with
      // the most single subclone bases fixed vs the one that has
      // the highest quality.  I think that we should go with Pat's
      // quality cut-off idea--just don't use a template at all if
      // it has too poor a quality.  If 2 templates are both above
      // the quality cut-off, then use the one that fixes the most
      // single subclone bases.  OK, so we could reorder the
      // subcloneTemplates array.  (We may need to reorder it more
      // than once.  Shall we remove a template entirely if it
      // fixes less than a minimum number of subclone bases?
      // Certainly if it fixes *no* single subclone bases, it
      // should be removed.  Shall the finishers be informed if the
      // 2nd fixes al


      subcloneTTemplate* pSubcloneTemplate = 
         pPrimer->pSubcloneTemplate_[ nTemplate ];

      if ( ! pSubcloneTemplate ) continue;

      RWCString soSubcloneName = pSubcloneTemplate->soTemplateName_;

      // fix bug found by Darren Grafham, Sanger:
      //  ssss TTTTTTTTTTTT  where sss is single subclone and TTT is
      // a template.  In making a bottom strand read, Autofinish used
      // to ignore the left end of the TTT and just assume that a read
      // length of bases would be made, thus fixing this single
      // subclone region

      
      int nUnpaddedLeftOnTemplate;
      int nUnpaddedRightOnTemplate;

      if (!bIntersect( nConsPosLeftUnpadded,
		      nConsPosRightUnpadded,
		      pSubcloneTemplate->nUnpaddedStart_,
		      pSubcloneTemplate->nUnpaddedEnd_,
		      nUnpaddedLeftOnTemplate,
		      nUnpaddedRightOnTemplate ) ) continue;
      
      for( int nUnpadded = nUnpaddedLeftOnTemplate; 
	  nUnpadded <= nUnpaddedRightOnTemplate;
	  ++nUnpadded ) {
            
         if ( !(*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] )
            ++nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ];
         else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 1 ) {
            LocatedFragment* pLocFrag = 
               (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ];

            if ( pLocFrag->soGetTemplateName() != soSubcloneName ) {
               ++nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ];
            }
         }
      }
   } // for( int nTemplate = ...

   // sort the templates based on how many single subclone bases
   // each solves


   sortTemplatesByHowManySingleSubcloneBasesEachFixes(
           pPrimer,
           nSingleSubcloneBasesFixedByThisTemplate );

   // make a pass through the array eliminating any template
   // that doesn't solve any single subclone bases

   eliminateTemplatesThatDoNotSolveAnySingleSubcloneBases(
           pPrimer,
           nSingleSubcloneBasesFixedByThisTemplate );

   // the # of ssb's fixed by the experiment is defined to be the
   // greatest # of ssb's fixed by any template of the experiment
   pExpToBeAdjusted->nSingleSubcloneBasesFixed_ = 
      nSingleSubcloneBasesFixedByThisTemplate[ 0 ];

   pExpToBeAdjusted->dSingleSubcloneBasesFixedPerDollar_ =
      ( (double) pExpToBeAdjusted->nSingleSubcloneBasesFixed_ ) 
      /
      ( (double) pExpToBeAdjusted->fCost_ );

} // if ( pExpToBeAdjusted->bCloneTemplateNotSubcloneTemplate_ )  else






void Contig :: adjustSingleSubcloneBasesFixedByWholeCloneExperiment(
                                       experiment* pExpToBeAdjusted,
                                       const int nConsPosLeftUnpadded,
                                       const int nConsPosRightUnpadded ) {
      
   int nNumberOfSingleSubcloneBasesFixed = 0;
   for( int nUnpadded = nConsPosLeftUnpadded; nUnpadded <= nConsPosRightUnpadded; ++nUnpadded ) {
      if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] < 2 ) 
         ++nNumberOfSingleSubcloneBasesFixed;
   }

   pExpToBeAdjusted->nSingleSubcloneBasesFixed_ = 
      nNumberOfSingleSubcloneBasesFixed;

   pExpToBeAdjusted->dSingleSubcloneBasesFixedPerDollar_ =
      ( (double) pExpToBeAdjusted->nSingleSubcloneBasesFixed_ )
      /
      ( (double) pExpToBeAdjusted->fCost_ );
}




experiment* Contig :: pGetBestSingleSubcloneExperiment() {
   experiment* pBestExperiment = NULL;

   bool bTryAgain = false;
   do {
      pBestExperiment = pGetBestSingleSubcloneExperimentButMightHaveBeenTriedBefore();

      if ( pBestExperiment ) {
         if ( pBestExperiment->bHasThisExperimentOrOneLikeItBeenTriedBefore() ) {
            pBestExperiment->nState_ = nRejectedExperiment;
            bTryAgain = true;
         }
         else 
            bTryAgain = false;
         
      }
      else
         // case in which there are no more available experiments
         bTryAgain = false;
   } while ( bTryAgain );


   return( pBestExperiment );
}








experiment* Contig :: pGetBestSingleSubcloneExperimentButMightHaveBeenTriedBefore() {

   experiment* pBestSoFar = NULL;
   double dSingleSubcloneBasesFixedPerDollarBestSoFar = 0.0;

   for( int nCan = 0; nCan < pCandidateExperiments_->length(); ++nCan ) {
      experiment* pExp = (*pCandidateExperiments_)[ nCan ];

      if ( pExp->nState_ == nChosenExperiment ) continue;
      if ( pExp->nState_ == nRejectedExperiment ) continue;

      if ( pExp->nSingleSubcloneBasesFixed_ == 0 ) continue;

      if ( pExp->dSingleSubcloneBasesFixedPerDollar_ >
           dSingleSubcloneBasesFixedPerDollarBestSoFar ) {

	 // check that, if the read is a subclone custom primer read,
	 // it isn't too close to another subclone custom primer read
	 // on the same strand that has already been chosen

	 if (pCP->bAutoFinishDoNotAllowSubcloneCustomPrimerReadsCloseTogether_ ) {
	    if ( pExp->bIsCustomPrimerSubcloneRead() &&
	      bThereIsAnotherSubcloneCustomPrimerReadTooCloseOnTheSameStrand( pExp ) ) {
	       // this experiment can never be considered again in this
	       // round of autofinish, so reject it
	       pExp->nState_ = nRejectedExperiment;
	       continue;
	    }
	 } // if (pCP->bAutoFinishDoNotAllowSubcloneCustomPrimer...

         
         pBestSoFar = pExp;
         dSingleSubcloneBasesFixedPerDollarBestSoFar =
            pExp->dSingleSubcloneBasesFixedPerDollar_;
      }
   }

   //   if ( pBestSoFar == NULL ) 
   // printf( "no best single subclone experiment\n" );


   // Note:  This might return NULL
   return( pBestSoFar );
}




void Contig :: autoFinishAdjustSingleSubcloneArrays( 
                                   experiment* pChosenExperiment ) {

   
   int nConsPosLeftUnpadded;
   int nConsPosRightUnpadded;
   
   if ( !bGetLikelyAlignedRegionOnContig( pChosenExperiment, 
                                       nConsPosLeftUnpadded,
                                       nConsPosRightUnpadded ) ) 
      return; // case in which the high quality region is not on the contig


   if ( pChosenExperiment->bCloneTemplateNotSubcloneTemplate_ )
      autoFinishAdjustSingleSubcloneArraysDueToPickingWholeCloneExperiment( pChosenExperiment );
   else
      autoFinishAdjustSingleSubcloneArraysDueToPickingSubcloneExperiment( pChosenExperiment );

}


void Contig :: autoFinishAdjustSingleSubcloneArraysDueToPickingSubcloneExperiment( experiment* pChosenExperiment ) {
   
   int nConsPosLeftUnpadded;
   int nConsPosRightUnpadded;
   
   if (!bGetLikelyAlignedRegionOnContig( pChosenExperiment,
                                       nConsPosLeftUnpadded,
                                       nConsPosRightUnpadded ) ) return;

   
   // I'm going to assume that the lab will only do some of the templates

   primerType* pPrimer = pChosenExperiment->pCustomPrimer_;
   for( int nTemplate = 0; 
        ( nTemplate < nNUMBER_OF_TEMPLATES ) &&
           ( nTemplate < consedParameters::pGetConsedParameters()->nAutoFinishHowManyTemplatesYouIntendToUseForCustomPrimerSubcloneReactions_ ) &&
           ( pPrimer->pSubcloneTemplate_[ nTemplate ] != NULL );
        ++nTemplate ) {

      subcloneTTemplate* pSubcloneTemplate = pPrimer->pSubcloneTemplate_[ nTemplate ];

      RWCString soSubcloneName = pSubcloneTemplate->soTemplateName_;


      // fix bug found by Darren Grafham, Sanger:
      //  ssss TTTTTTTTTTTT  where sss is single subclone and TTT is
      // a template.  In making a bottom strand read, Autofinish used
      // to ignore the left end of the TTT and just assume that a read
      // length of bases would be made, thus fixing this single
      // subclone region

      
      int nUnpaddedLeftOnTemplate;
      int nUnpaddedRightOnTemplate;

      if (!bIntersect( nConsPosLeftUnpadded,
		      nConsPosRightUnpadded,
		      pSubcloneTemplate->nUnpaddedStart_,
		      pSubcloneTemplate->nUnpaddedEnd_,
		      nUnpaddedLeftOnTemplate,
		      nUnpaddedRightOnTemplate ) ) continue;
      
      for( int nUnpadded = nUnpaddedLeftOnTemplate; 
	  nUnpadded <= nUnpaddedRightOnTemplate;
	  ++nUnpadded ) {

         if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 0 ) {
            // there is no aligned read so far covering this position
            (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 1;
            (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ] = 
               pSubcloneTemplate->aReads_[0];
            // just picked an arbitrary read to represent the template
         }
         else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 1 ) {
            // there is some other read covering this base
            // or this base is covered by some read that has
            // been suggested by autofinish

            // Let's see if it is the same subclone

            LocatedFragment* pLocFrag2 = 
               (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ];

            if ( pLocFrag2->soGetTemplateName() != soSubcloneName ) {
               (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 2;
            }
         }

         // ignore case in which there are already 2 or more subclones 
         // covering this base
      } // for( int nUnpadded ...
   } // for( int nTemplate...
}
         




void Contig :: autoFinishAdjustSingleSubcloneArraysDueToPickingWholeCloneExperiment( experiment* pChosenExperiment ) {

   int nConsPosLeftUnpadded;
   int nConsPosRightUnpadded;

   if (!bGetLikelyAlignedRegionOnContig( pChosenExperiment,
                                       nConsPosLeftUnpadded,
                                       nConsPosRightUnpadded ) ) return;


   for( int nUnpadded = nConsPosLeftUnpadded; nUnpadded <= nConsPosRightUnpadded; ++nUnpadded ) {
      (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = nMaxSubcloneCoverage;
   }
}







int Contig :: nCountSingleSubcloneBases() {
   
   int nNumberOfSingleSubcloneBases = 0;
   for( int nUnpadded = nGetUnpaddedStartIndex(); nUnpadded <= nGetUnpaddedEndIndex(); ++nUnpadded ) {
      if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] < 2 )
         ++nNumberOfSingleSubcloneBases;
   }

   return( nNumberOfSingleSubcloneBases );
}



void Contig :: autoFinishPrintSingleSubcloneRegions() {


   bool bInSingleSubcloneRegion = false;
   int nSingleSubcloneRegionStart;
   for( int nUnpadded = nGetUnpaddedStartIndex(); nUnpadded <= nGetUnpaddedEndIndex(); ++nUnpadded ) {
      if ( bInSingleSubcloneRegion ) {
         if  ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] >= 2 ) {
            bInSingleSubcloneRegion = false;

            int nSingleSubcloneRegionEnd = nUnpadded - 1;
            
            fprintf( pAO, "  Single subclone region from %d to %d (%d bp)\n",
               nSingleSubcloneRegionStart,
               nSingleSubcloneRegionEnd,
                    nSingleSubcloneRegionEnd -
                    nSingleSubcloneRegionStart + 1
                    );
         }
      }
      else {
         if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] < 2 ) {
            bInSingleSubcloneRegion = true;
            nSingleSubcloneRegionStart = nUnpadded;
         }
      }
   }

   // handle case in which single subclone region goes to the end

   if ( bInSingleSubcloneRegion ) {
      int nSingleSubcloneRegionEnd = nGetUnpaddedEndIndex();
      

      fprintf( pAO, "  Single subclone region from %d to %d (%d bp)\n",
         nSingleSubcloneRegionStart,
         nSingleSubcloneRegionEnd,
              nSingleSubcloneRegionEnd -
              nSingleSubcloneRegionStart + 1
              );
   }
}




void Contig :: orderTemplatesAndSetSingleSubcloneValues( experiment* pExp ) {

   int nConsPosLeftUnpadded;
   int nConsPosRightUnpadded;

   if ( !bGetLikelyAlignedRegionOnContig( pExp, nConsPosLeftUnpadded, nConsPosRightUnpadded ) ) {
      pExp->nSingleSubcloneBasesFixed_ = 0;
      pExp->dSingleSubcloneBasesFixedPerDollar_ = 0.0;
      return;
   }


   if ( pExp->bCloneTemplateNotSubcloneTemplate_ ) {
      // note that when counting single subclone bases this way,
      // if there is a single subclone base that has no aligned
      // read covering it (Phil on 12/23/98 said this is possible),
      // then this is counted as fixing a single base even though
      // it would take 2 M13 reads to cover the same base.
      // Thus there is a slight bias in favor of M13 reads.

      pExp->nSingleSubcloneBasesFixed_ = 
         (*pCumNumberOfSingleSubcloneBases_)[ nConsPosRightUnpadded ] 
         -
         (*pCumNumberOfSingleSubcloneBases_)[ nConsPosLeftUnpadded - 1 ];

      pExp->dSingleSubcloneBasesFixedPerDollar_ =
         ( (double) pExp->nSingleSubcloneBasesFixed_ )
         /
         ( (double) pExp->fCost_ );
      
   }
   else {

      primerType* pPrimer = pExp->pCustomPrimer_;
      int nSingleSubcloneBasesFixedByThisTemplate[ nNUMBER_OF_TEMPLATES ];

      int nTemplate;
      for( nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
         nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ] = 0;
      }

      for( nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
         subcloneTTemplate* pSubcloneTemplate = 
            pPrimer->pSubcloneTemplate_[ nTemplate ];

         if ( ! pSubcloneTemplate ) continue;

	 RWCString soSubcloneName = pSubcloneTemplate->soTemplateName_;

	 // fix bug found by Darren Grafham, Sanger:
	 //  ssss TTTTTTTTTTTT  where sss is single subclone and TTT is
	 // a template.  In making a bottom strand read, Autofinish used
	 // to ignore the left end of the TTT and just assume that a read
	 // length of bases would be made, thus fixing this single
	 // subclone region

	 
	 int nUnpaddedLeftOnTemplate;
	 int nUnpaddedRightOnTemplate;

	 if (!bIntersect( nConsPosLeftUnpadded,
	                  nConsPosRightUnpadded,
			  pSubcloneTemplate->nUnpaddedStart_,
			  pSubcloneTemplate->nUnpaddedEnd_,
			  nUnpaddedLeftOnTemplate,
			  nUnpaddedRightOnTemplate ) ) continue;
         
         for( int nUnpadded = nUnpaddedLeftOnTemplate; 
	     nUnpadded <= nUnpaddedRightOnTemplate;
              ++nUnpadded ) {

            if ( !(*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] ) {
               // just an approximation since it will take 2 M13 reads to
               // completely fix the single subclone base
               if ( nUnpadded % 2 == 0 )
                  ++nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ];
            }
            else if ( (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 1 ) {
               LocatedFragment* pLocFrag = 
                  (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ];

               if ( pLocFrag->soGetTemplateName() != soSubcloneName ) {
                  ++nSingleSubcloneBasesFixedByThisTemplate[ nTemplate ];
               }
            }
         } // for( int nUnpadded ...
      } // for( int nTemplate ...

      sortTemplatesByHowManySingleSubcloneBasesEachFixes( 
                         pPrimer,
                         nSingleSubcloneBasesFixedByThisTemplate
                                                          );

      eliminateTemplatesThatDoNotSolveAnySingleSubcloneBases(
                         pPrimer,
                         nSingleSubcloneBasesFixedByThisTemplate
                                                          );

      // use the template that fixes the most single subclone bases
      pExp->nSingleSubcloneBasesFixed_ = 
         nSingleSubcloneBasesFixedByThisTemplate[0];
      
      pExp->dSingleSubcloneBasesFixedPerDollar_ =
         ( (double) pExp->nSingleSubcloneBasesFixed_ )
         /
         ( (double) pExp->fCost_ );

   } //    if ( pExp->bCloneTemplateNotSubcloneTemplate_ ) { else

}


void Contig :: autoFinishExamineDoNotFinishTagsForSingleSubcloneRegions() {
   
   // look for do-not-finish tags in both reads and in the consensus

   // First look in the reads.

   for( int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
      LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );

      for( int nTag = 0; nTag < pLocFrag->nGetNumberOfTags(); ++nTag ) {
         tag* pTag = pLocFrag->pGetTag( nTag );

         if ( pTag->soGetTagType() == "cloneEnd" ) 
            dealWithCloneEndTagForSingleSubcloneRegions( pTag );
         else if ( consedParameters::pGetConsedParameters()->aAutoFinishTagsToNotFinish_.index( pTag->soGetTagType() ) == RW_NPOS ) continue;

         int nUnpaddedConsPosStart = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
         int nUnpaddedConsPosEnd = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );

         // read tags might be on a read that protrudes beyond the
         // end of the consensus

         int nOverlapLeft;
         int nOverlapRight;
         if ( bIntersect( nUnpaddedConsPosStart,
                          nUnpaddedConsPosEnd,
                          pNumberOfSubclonesCoveringEachBase_->nGetStartIndex(),
                          pNumberOfSubclonesCoveringEachBase_->nGetEndIndex(),
                          nOverlapLeft,
                          nOverlapRight ) ) {
            
            for( int nUnpadded = nOverlapLeft;
                 nUnpadded <= nOverlapRight; 
                 ++nUnpadded ) {
               (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 
                  nMaxSubcloneCoverage;            
            }
         }
      }
   }


   // now look for consensus tags

   for( int nConsensusTag = 0; nConsensusTag < nGetNumberOfTags(); 
        ++nConsensusTag ) {

      tag* pTag = pGetTag( nConsensusTag );

      if ( pTag->soGetTagType() == "cloneEnd" )
         dealWithCloneEndTagForSingleSubcloneRegions( pTag );
      else if ( consedParameters::pGetConsedParameters()->aAutoFinishTagsToNotFinish_.index( pTag->soGetTagType() ) == RW_NPOS ) continue;

      int nUnpaddedConsPosStart = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
      int nUnpaddedConsPosEnd = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );

      for( int nUnpadded = nUnpaddedConsPosStart;
           nUnpadded <= nUnpaddedConsPosEnd;
           ++nUnpadded ) {
         (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] =
            nMaxSubcloneCoverage;
      }
   }
}



void Contig :: dealWithCloneEndTagForSingleSubcloneRegions( tag* pTag ) {

   int nUnpaddedLeft;
   int nUnpaddedRight;
   if ( pTag->soMiscData_ == "->" ) {
      nUnpaddedLeft = pNumberOfSubclonesCoveringEachBase_->nGetStartIndex();
      nUnpaddedRight = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
   }
   else if ( pTag->soMiscData_ == "<-" ) {
      nUnpaddedLeft = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
      nUnpaddedRight = pNumberOfSubclonesCoveringEachBase_->nGetEndIndex();
   }
   else 
      return; // corrupted cloneEnd tag

   for( int nUnpadded = nUnpaddedLeft; nUnpadded <= nUnpaddedRight; 
        ++nUnpadded ) {
      (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] = 
         nMaxSubcloneCoverage;
   }
         


}






static int nWindowOfSingleSubcloneBasesFixed[10];


void Contig :: findWindowOfMostSingleSubcloneBasesFixed( experiment* pExp ) {

   int nSingleSubcloneBasesFixedInWindowOfMostFixed = 0;
   int nSingleSubcloneBasesFixedInCurrentWindow = 0;
   int nWindowOfMostSingleSubcloneBasesFixedUnpaddedRight;


   int nPointerToLastAddedBase = -1;

   int nUnpaddedLeftOnContig;
   int nUnpaddedRightOnContig;

   if ( !bIntersect( pNumberOfSubclonesCoveringEachBase_->nGetStartIndex(),
                     pNumberOfSubclonesCoveringEachBase_->nGetEndIndex(),
                     pExp->nUnpaddedLeft_,
                     pExp->nUnpaddedRight_,
                     nUnpaddedLeftOnContig,
                     nUnpaddedRightOnContig ) )
      return;

   // for the case of using a subclone walk, we will just consider the
   // number of single subclone bases fixed by the first template.  I'm
   // being a little sloppy here, but this window is only used for people
   // looking at the results, and has no effect on which reads are chosen.


   RWCString soSubcloneName;

   if ( pExp->bIsCustomPrimerSubcloneRead() ) {
      
      subcloneTTemplate* pSub = pExp->pCustomPrimer_->pSubcloneTemplate_[0];

      if (!bIntersect( 
                nUnpaddedLeftOnContig,
                nUnpaddedRightOnContig,
                pSub->nUnpaddedStart_,
                pSub->nUnpaddedEnd_,
                nUnpaddedLeftOnContig,
                nUnpaddedRightOnContig ) )
         return;

      soSubcloneName = pSub->soTemplateName_;
   }
   else
      soSubcloneName = "";
                       


   int n10thBaseUnpadded = 
      nUnpaddedLeftOnContig + 9;


   for( int nUnpadded = nUnpaddedLeftOnContig; 
        nUnpadded <= nUnpaddedRightOnContig; ++nUnpadded ) {


      int nNumberOfSingleSubcloneBasesFixedAtBase = 0;
      if (  (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 0 )
         nNumberOfSingleSubcloneBasesFixedAtBase = 1;
      else if (  (*pNumberOfSubclonesCoveringEachBase_)[ nUnpadded ] == 1 ) {
         if ( (*pWhichSingleSubcloneIsCoveringEachBase_)[ nUnpadded ]->soGetTemplateName() 
              != soSubcloneName )
            nNumberOfSingleSubcloneBasesFixedAtBase = 1;
      }
              
      

      if ( nUnpadded <= n10thBaseUnpadded ) {

         ++nPointerToLastAddedBase;
         
         nWindowOfSingleSubcloneBasesFixed[ nPointerToLastAddedBase ] =
            nNumberOfSingleSubcloneBasesFixedAtBase;

         nSingleSubcloneBasesFixedInCurrentWindow +=
            nNumberOfSingleSubcloneBasesFixedAtBase;

         if ( nUnpadded == n10thBaseUnpadded ) {
            nSingleSubcloneBasesFixedInWindowOfMostFixed = 
               nSingleSubcloneBasesFixedInCurrentWindow;

            nWindowOfMostSingleSubcloneBasesFixedUnpaddedRight =
               nUnpadded;
         }
      }
      else {

         ++nPointerToLastAddedBase;
         if ( nPointerToLastAddedBase > 9 )
            nPointerToLastAddedBase = 0;

         nSingleSubcloneBasesFixedInCurrentWindow 
            -= nWindowOfSingleSubcloneBasesFixed[ nPointerToLastAddedBase ];
         
         nSingleSubcloneBasesFixedInCurrentWindow
            += nNumberOfSingleSubcloneBasesFixedAtBase;

         nWindowOfSingleSubcloneBasesFixed[ nPointerToLastAddedBase ] 
            = nNumberOfSingleSubcloneBasesFixedAtBase;

         
         if ( nSingleSubcloneBasesFixedInCurrentWindow >
              nSingleSubcloneBasesFixedInWindowOfMostFixed ) {
            nSingleSubcloneBasesFixedInWindowOfMostFixed =
               nSingleSubcloneBasesFixedInCurrentWindow;
            nWindowOfMostSingleSubcloneBasesFixedUnpaddedRight = nUnpadded;
         }
      }
   }


   pExp->nWindowOfMostSingleSubcloneBasesFixedUnpaddedRight_ = 
      nWindowOfMostSingleSubcloneBasesFixedUnpaddedRight;

   pExp->nSingleSubcloneBasesFixedInWindow_ =
      nSingleSubcloneBasesFixedInWindowOfMostFixed;


}