/**
 ** Copyright (c) 2007-2010 Illumina, Inc.
 **
 ** This software is covered by the "Illumina Genome Analyzer Software
 ** License Agreement" and the "Illumina Source Code License Agreement",
 ** and certain third party copyright/licenses, and any user of this
 ** source file is bound by the terms therein (see accompanying files
 ** Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
 ** Illumina_Source_Code_License_Agreement.pdf and third party
 ** copyright/license notices).
 **
 ** This file is part of the Consensus Assessment of Sequence And VAriation
 ** (CASAVA) software package.
 **
 ** \file phageAlign.cpp
 **
 ** \brief Phage align.
 **
 ** Phage align.
 **
 ** \author Tony Cox
 **/

#include <iostream>
#include <vector>
#include <algorithm>
#include <fstream>
#include <cassert>
#include <cctype>
#include <sys/time.h>
#include <fcntl.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <unistd.h>
#include <fstream>
#include <iostream>
#include <sstream>
#include <deque>
#include <map>
#include <dirent.h>
#include <errno.h>
#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <string>
#include <ctime>

#include <boost/filesystem.hpp>

#include "common/MathCompatibility.hh"
#include "alignment/GlobalUtilities.hh"
#include "common/Alignment.hh"
#include "common/Exceptions.hh"
#include "alignment/MatchDescriptor.hh"

#define SWITCHED_TO_EXPORT_OUTPUT 1

// Simple program to match sequences against a small genome allowing for an
// arbitrary number of errors

using namespace std;
typedef unsigned int uint;
//const uint lineLength(50);


//const uint bufSize(128);
//char oligoBuffer[bufSize];
//const uint maxSeqSize(60);
//const uint numBases(4);
// Assume will encounter A, C, G, T or '.' in real sequence
const uint numReadableBases(numDifferentBases+1);

const int veryBadScore(-99999999);


const short scoreNotSet=-20000;

typedef short ScoreTable[maxSeqSize][numReadableBases][numDifferentBases];


// +20 because title line can overflow sequence line length

// Comment out following line to disable support for circularized genome
#define CIRCULARIZE
// Uncomment out following line to only search forward strand
//#define DONT_SEARCH_REVERSE_STRAND
// Uncomment next line to scan for indels as well as substitution errors
//#define CHECK_FOR_INDELS

char titleBuf[256];

// NB scores are stored as shorts, so there will be problems
// if max(defaultMatch, defaultMismatch)*seqSize>32767.
// defaultMatch=500 OK for read lengths<=65
// Assumed scoring scheme= 100*log10(odds of outcome);
const short defaultMatch=500;
const short defaultMismatch=-547;
// defaultMatch=500 is equivalent to correct:wrong=100000:1
// assuming the three possible wrong bases for any sequenced base
// are equiprobable, then the corresponding defaultMismatch is
// 100*log10(1/300000)=-547
const short defaultBlank=-47;
// "-47" comes from 100*log10(1/3) - log odds score if all four bases
// are equally likely to be sequenced as a blank


struct GenomeFrag
{
  GenomeFrag( const char* seq, const int pos ) :
    seq_(seq), pos_(pos) {}
  const char* seq_;
  uint pos_;
}; // ~GenomeFrag

bool lessThanGenomeFrag( const GenomeFrag& a, const GenomeFrag& b )
{
  return (strcmp(a.seq_, b.seq_)<0);
} // ~bool lessThanGenomeFrag

const uint reverseFlag(((uint)1)<<31);

// Sort GenomeFrags using a least-significant-digit first radix sort
void mySort( vector<GenomeFrag>& v, const int seqSize )
{
  // Written somewhat quickly, probably does more copying of data than
  // it needs to
  vector<GenomeFrag> vtemp(v);
  vector<GenomeFrag*> v1, v2;
  for (uint i(0);i<vtemp.size();i++)
    v1.push_back(&vtemp[i]);

  map<char,vector<GenomeFrag*> > m;

  for (int i(seqSize-1);i>=0;i--)
  {
    // bin entries v1 by seqSize-th digit
    for (uint j(0); j<vtemp.size(); j++)
    {
      m[(*v1[j]).seq_[i]].push_back(v1[j]);
    } // ~for

    v2.clear();
    // traverse m in order and place bin contents in v2
    // Upshot: v2 is sorted by seqSize-th digit
    for (map<char,vector<GenomeFrag*> >::iterator i(m.begin());i!=m.end();++i)
    {
      if (!i->second.empty())
      {
        v2.insert(v2.end(),i->second.begin(),i->second.end());
        i->second.clear();
      } // ~if
    } // ~for
    v1=v2;
  } // ~for

  // Now overwrite v with the finished product
  for (uint i(0);i<vtemp.size();i++)
    v[i]=*(v1[i]);

  return;
} // ~mySort






struct Seq
{
  //  Seq() {}
  char data[maxSeqSize];
}; // ~struct Seq

//typedef int PlusScore;
typedef unsigned short PlusScore;

const PlusScore badLess(60000);

typedef vector<PlusScore> ScoreByBase;
typedef vector<ScoreByBase> ScoreByCycle;
typedef vector<ScoreByCycle> ScoreBySeq;
struct QuerySeqs
{
  QuerySeqs( int seqSize ) :
    scores_( seqSize+1, ScoreByCycle(4) ) {}

  void addSeq( const OligoSource& oligos,
               const ScoreSource& scores,
               const uint seqSize );

  ScoreBySeq scores_;
  vector<BaseScore> maxScore_;
}; // ~struct QuerySeqs

void QuerySeqs::addSeq(const OligoSource& ,
                       const ScoreSource& scores,
                       const uint seqSize)
{


  BaseScore t1,t2;
  //  strncpy( data_, oligos.getLastOligo(), seqSize );
  //   for (uint i(0); i<seqSize; i++)
  //   {
  //    data_[i]=oligos.getLastOligo()[i];
  //    } // ~for i
  //  data_[seqSize]='\0';
  //   cout << data_ << endl;
  maxScore_.push_back(0);
  for (uint i(0); i<seqSize; i++)
  {
    t1=scores.getScore((uint)0,i);
    for (uint j(1); j<numDifferentBases; j++)
    {
      t2=scores.getScore(j,i);
      if (t2>t1) t1=t2;
    } // ~for j
    maxScore_.back()+=t1;
    //    cout << i << '\t' << t1;
    for (uint j(0); j<numDifferentBases; j++)
    {
      //      scores_[j][i]=t1-scores.getScore(j,i);
      scores_[i+1][j].push_back(t1-scores.getScore(j,i));
      //     cout << '\t' << scores_[j][i];
    } // ~for j
    //    cout << endl;
  } // ~for i

} // ~QuerySeqs::addSeq

struct QueryResults
{
  QueryResults( int seqSize, int numSeqs ) :
    cycle_( seqSize+1, ScoreByBase(numSeqs)),
    bestLess_(numSeqs, badLess),
    nextBestLess_(numSeqs, badLess),
    numBest_(numSeqs,0),
    bestPos_(numSeqs,0)
  {} // ~ctor

  ScoreByCycle cycle_;
  ScoreByBase bestLess_;
  ScoreByBase nextBestLess_;
  ScoreByBase numBest_;
  vector<uint> bestPos_;
}; // ~struct QueryResults


struct QuerySeq
{
  QuerySeq( const OligoSource& oligos,
            const ScoreSource& scores,
            const uint seqSize,
            const uint seqNum);


  char data_[maxSeqSize];
}; // ~struct QuerySeq

QuerySeq::QuerySeq( const OligoSource& oligos,
                    const ScoreSource& ,
                    const uint seqSize,
                    const uint )
{
  strncpy( data_, oligos.getLastOligo(), seqSize );
  data_[seqSize]='\0';

} // ~QuerySeq::c'tor



bool lessThanQuerySeq( const QuerySeq& a, const QuerySeq& b )
{
  return (strcmp(a.data_, b.data_)<0);
} // ~bool lessThanQuerySeq

// returns first position at which a and b differ (or seqSize if identical)
int firstDiff( const QuerySeq& a, const QuerySeq& b, const int seqSize )
{
  int i(0);
  for (;i<seqSize;i++)
  {
    if (a.data_[i]!=b.data_[i]) break;
  } // ~for
  return i;
} // ~int firstDiff

// returns first position at which a and b differ (or seqSize if identical)
int firstDiff( const char* a, const char* b, const int seqSize )
{
  int i(0);
  for (;i<seqSize;i++)
  {
    if (a[i]!=b[i]) break;
  } // ~for
  return i;
} // ~int firstDiff

/*****************************************************************************/
// tempFragStore must persist outside of this function because it contains
// the actual fragment strings; frags contains only pointers to them.

void storeTags(const char* tagFileName, const uint seqSize,
               vector<GenomeFrag>& frags, vector<Seq>& tempFragStore)
{
    OligoSource* pTags(OligoSource::getOligoSource(tagFileName));
    cerr << "Opened tag file " << tagFileName << endl;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf( "#TAG_FILE %s\n", tagFileName );
#endif

    const char* tagBuf;

    while((tagBuf=pTags->getNextOligo())!=NULL)
    {
        tempFragStore.push_back( Seq() );

        if (strlen(tagBuf)<seqSize)
        {
            cerr << "ERROR: all tags in tag file " << tagFileName
                 << " need to be at least " << seqSize
                 << " bases long, but tag " << tagBuf << " (line "
                 << tempFragStore.size() << ") is only " << strlen(tagBuf)
                 << " bases long." << endl;
            exit (-1);
        } // ~if

        strncpy(tempFragStore.back().data, tagBuf, seqSize);
        tempFragStore.back().data[seqSize]='\0';
    } // ~while

    frags.reserve(tempFragStore.size());

    for (uint tagNum(0);tagNum<tempFragStore.size();tagNum++)
    {
        frags.push_back(GenomeFrag( tempFragStore[tagNum].data, tagNum));
    } // ~for

    delete pTags;
}

/*****************************************************************************/
// FIXME : Use_Bases : expandLegacyUseBasesString has similar functionality.

void expandUseBasesStr(const std::string& useBasesStr,
                       std::string& expandedUseBasesStr,
                       const uint seqSize)
{
    expandedUseBasesStr.clear();

    if (useBasesStr == "all")
    {
        expandedUseBasesStr.append(seqSize, 'Y');
    } // ~if
    else if (useBasesStr == "odd")
    {
        for (uint i(0); i<2*seqSize-1; i++) {
            expandedUseBasesStr += (((i&0x1)==0) ? 'Y' : 'N');
        }
    }
    else
    {
        // Assume custom useBases argument, error check it
        const uint useBasesLen(useBasesStr.length());
        uint numYs(0), numNs(0);

        for (std::string::const_iterator chCIter(useBasesStr.begin());
             chCIter != useBasesStr.end(); ++chCIter)
        {
            numYs += (toupper(*chCIter) == 'Y');
            numNs += (toupper(*chCIter) == 'N');
        } // ~for

        if (numYs + numNs < useBasesLen)
        {
            cerr << "Error: could not parse -useBases argument "
                 << useBasesStr
                 << endl
                 << "Must be 'all', 'odd' or some combination of y,Y,n,N"
                 << endl;
            exit (1);
        } // ~if

        if (numYs<seqSize)
        {
            cerr << "Error: expecting " << seqSize
                 << " bases, -useBases argument only specifies "
                 << numYs << endl;
            exit (1);
        } // ~if

        expandedUseBasesStr = useBasesStr;
    } // ~else if

#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf("#BASES_USED %s\n", expandedUseBasesStr.c_str());
#endif
}

/*****************************************************************************/

int phageAlign( int numArgs, char* args [] )
{
  if (numArgs<5)
  {
    cerr << "syntax: " << args[0]
	 << " [-tag] genomeFile maxNumBlanks sequenceLength "
	 << endl
	 << "[-useBases basesToUse] [-scoreFile scoreFileName] sequenceFile[s]"
	 << endl
         << "-tag: genomeFile is a file of sequence tags" << endl
	 << "genomeFile: single entry FASTA format file" << endl
	 << "maxNumBlanks: max number of non-called (.) bases to tolerate"
	 << endl
	 << "sequenceLength: number of cycles to process" << endl
         << "-useBases: specify which bases to align/ignore" << endl
         << "-useBases odd: align odd-numbered bases only (default)" << endl
         << "-useBases all: align all bases" << endl
         << "-useBases YYNNYYY...: Y=align base, N=don't" << endl
         << "-scoreFile: read per-cycle scores from scoreFileName" << endl
	 << "sequenceFiles: ASCII file[s] of sequence data"
	 << endl;
    exit (1);
  } // ~if

  int thisArg=1;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
  time_t timeNow = time(NULL);
  printf( "#RUN_TIME %s", ctime(&timeNow) );
  printf( "#SOFTWARE_VERSION %s\n", id );
#endif

  int genomeSize(0);
  vector<GenomeFrag> frags;
  vector<Seq> tempFragStore;
  vector<Seq> fragStore;
  char* buf(NULL);
  char* bufR(NULL);
  char* tagFileName(NULL);
  char* genomeFileName(NULL);

  if (strcmp( args[thisArg], "-tag")==0)
  {
    tagFileName=args[++thisArg];
  } // ~if
  else
  {
    genomeFileName=args[thisArg];
  } // ~else

  uint maxNumBlanks(atoi(args[++thisArg]));

  cerr << "Will ignore sequences containing more than " << maxNumBlanks
       << " blanks" << endl;
#ifndef SWITCHED_TO_EXPORT_OUTPUT
  printf( "#MAX_BLANKS %d\n", maxNumBlanks );
#endif

  uint seqSize =atoi(args[++thisArg]);

  cerr << "Assuming sequence length of " << seqSize << " bases" << endl;
#ifndef SWITCHED_TO_EXPORT_OUTPUT
  printf( "#SEQ_LENGTH %d\n", seqSize );
#endif

  //  assert(seqSize<maxSeqSize);
  if ((seqSize==0)||(seqSize>=maxSeqSize))
  {
    cerr << "Error: specified sequence size (" << seqSize
         << ") is outside valid range (1 to " << maxSeqSize-1 << ")" << endl;
    exit (1);
  }

  std::string chrName;

  if (tagFileName!=NULL)
  {
      // sequences are matched against a set of tags

      // tempFragStore contains the actual strings - frags just points to them.
      storeTags(tagFileName, seqSize, frags, tempFragStore);
      chrName = boost::filesystem::path(tagFileName).leaf();
  } // ~if
  else
  {
    // assume target sequence is a single entry FASTA format file
    assert(genomeFileName!=NULL);

    FILE* pGenome;

    if ((pGenome=fopen(genomeFileName, "r" ))==NULL)
    {
      cerr << "Error: could not open genome file "
           << args[1] << endl;
      exit (1);
    } // ~if

    cerr << "Opened genome file " << genomeFileName << endl;
#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf( "#GENOME_FILE %s\n", genomeFileName );
#endif
    chrName = boost::filesystem::path(genomeFileName).leaf();

    int numRead(0);

    // read comment line;
    assert(fgets(titleBuf, 256, pGenome)!=NULL);
    assert(titleBuf[0]=='>');
    titleBuf[strlen(titleBuf)-1]='\0';

    if (titleBuf[strlen(titleBuf)-1]=='\r')
    {
      cerr << "ERROR: genome file " << genomeFileName
	   << " is a DOS format text file!!\n";
      exit (-1);
    } // if


    // Create a buffer to hold the genome data
    numRead=ftell(pGenome);
    fseek(pGenome,0,SEEK_END);
    genomeSize=ftell(pGenome)-numRead+seqSize;
    cerr << "Creating buffer of size " << genomeSize << " bytes" << endl;

    buf=new char[genomeSize];
    char* p(buf);
    fseek(pGenome,numRead,SEEK_SET);
    genomeSize=0;


    while (fgets(p, 256, pGenome)!=NULL)
    {
      numRead=strlen(p)-1; // -1 ignores carriage return at end
      p+=numRead;
      genomeSize+=numRead;
    } // ~while

    fclose(pGenome);

    cerr << "Read in genome " << titleBuf << " of size " << genomeSize
         << endl;

#ifdef CIRCULARIZE
    for (uint i(0);i<seqSize-1;i++)
    {
      *p++=buf[i];
    } // ~for
    cerr << "Circular genome assumed: Any matches with position >"
         << genomeSize-seqSize+1
         << " cross cut point" << endl;
    genomeSize+=seqSize-1;
#endif

    *p++='\0';

    frags.reserve(genomeSize-seqSize+1);

    bool isValid;

    for (uint i(0);i<=genomeSize-seqSize; i++)
    {
      // don't match if fragment contains non-base characters, e.g. Ns
      isValid=true;
      for (uint j(seqSize); ((j!=0)&&isValid);)
      {
        --j;
        isValid=(realBaseCodes[(uint)*(buf+i+j)]!=nc);
      } // ~for

      if (isValid)
      {
        frags.push_back( GenomeFrag( buf+i, i));
      } // ~if

    } // ~if

#ifndef DONT_SEARCH_REVERSE_STRAND

    frags.reserve(2*frags.size());

    bufR = new char[genomeSize+1];
    bufR[genomeSize]='\0';

    // reverse complement the genome not the queries
    for (int i(0);i<genomeSize;i++)
      bufR[genomeSize-1-i]=reverseCharASCII[(uint)buf[i]];

    for (uint i(0);i<=genomeSize-seqSize; i++)
    {
      // don't match if fragment contains non-base characters, e.g. Ns
      isValid=true;
      for (uint j(seqSize); ((j!=0)&&isValid);)
      {
        --j;
        isValid=(realBaseCodes[(uint)*(bufR+i+j)]!=nc);
      } // ~for
      if (isValid)
      {
        frags.push_back( GenomeFrag( bufR+i, i|reverseFlag));
      } // ~if
    } // ~if

#endif

  } // ~else


  ++thisArg;

  std::string expandedUseBasesStr;
  std::string useBasesStr("odd"); // default

  if (strcmp( args[thisArg], "-useBases")==0)
  {
    ++thisArg; // onto the USE_BASES
    assert(thisArg<numArgs);
    useBasesStr = args[thisArg];

    ++thisArg; // onto next option
    assert(thisArg<numArgs);
  }

  expandUseBasesStr(useBasesStr, expandedUseBasesStr, seqSize);

  ScoreSource* pScoresRaw(NULL);
  ScoreSourceFilter scoreFilter(expandedUseBasesStr.c_str());
  ScoreSource* pScores(&scoreFilter);


  const char* scoreFileName(NULL);

  if (strcmp( args[thisArg], "-scoreFile")==0)
  {

    ++thisArg;
    assert(thisArg<numArgs);
    cerr << "Reading score information from file " << args[thisArg]
	 << endl;

    scoreFileName=args[thisArg];

#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf( "#SCORE_FILE %s\n", args[thisArg] );
#endif

    ++thisArg;
    assert(thisArg<numArgs);

  } // ~if
  else
  {
    cerr << "Using default scores of " << defaultMatch
	 << " for match, " << defaultMismatch << " for mismatch"
	 << endl;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf( "#SCORE_FILE none\n");
#endif
  } // ~else


  uint numIgnored(0);
  uint seqNum(0);

  OligoSource* pOligosRaw(NULL);
  OligoSourceFilter oligoFilter(expandedUseBasesStr.c_str());
  OligoSource* pOligos(&oligoFilter);

  vector<casava::common::Sequence> sequenceVec;
  vector<QuerySeq> seqs;
  QuerySeqs seqs_(seqSize);


  for (int i(thisArg); i<numArgs; i++)
  {

    pOligosRaw=OligoSource::getOligoSource(args[i]);

    cerr << "Opened sequence file " << args[i] << endl;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
    printf ("#SEQ_FILE %s\n", args[i]);
#endif

    if (scoreFileName!=NULL)
    {
      oligoFilter.link(pOligosRaw);
      pScoresRaw=new ScoreSourceCycle( *pOligos, scoreFileName );
      pScores=pScoresRaw;
      //      pScoresRaw=new ScoreSourceCycle( *pOligosRaw, scoreFileName );

#ifndef SWITCHED_TO_EXPORT_OUTPUT
      printf ("#SCORE_FILE %s\n", scoreFileName);
#endif
    } // ~if
    else
    {
      pScoresRaw=dynamic_cast<ScoreSource*>(pOligosRaw);
      if (pScoresRaw!=NULL)
      {
#ifndef SWITCHED_TO_EXPORT_OUTPUT
        printf ("#SCORE_FILE %s\n", args[i]);
#endif

        scoreFilter.link(pScoresRaw);
        oligoFilter.link(pOligosRaw);
        pScores=&scoreFilter;
        pOligos=&oligoFilter;
      } // ~if
      else
      {
        //        pScoresRaw=new ScoreSourceBasic
        pScoresRaw=new ScoreSourceBasic
          ( *pOligos, defaultMatch, defaultBlank, defaultMismatch );
        pScores=pScoresRaw;
        oligoFilter.link(pOligosRaw);
        pOligos=&oligoFilter;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
        printf ("#SCORE_FILE none\n");
#endif
      } // ~else
    } // ~else

    // pScoresRaw and pOligosRaw never accessed directly from now on
    // ...until they are no longer needed and are deleted

    bool isValid(true);
    casava::common::Sequence sequence(pOligos->getNextSequence(isValid));

    while (isValid)
    {
      if (sequence.getData().size() < seqSize)
      {
        cerr << "ERROR: all query sequences in file " << args[i]
             << " need to be at least " << seqSize
             << " bases long, but sequence " << sequence.getData()
             << " (seq num " << seqNum + 1 << ") is only "
             << sequence.getData().size() << " bases long." << endl;
        exit (-1);
      } // ~if

      std::string sequenceData = sequence.getData();
      std::replace(sequenceData.begin(), sequenceData.end(), '.', 'N');
      sequence.setData(sequenceData);
      sequenceVec.push_back(sequence);
      seqs_.addSeq( *pOligos, *pScores, seqSize );
      seqs.push_back( QuerySeq( *pOligos, *pScores, seqSize, seqNum++));
      // cout << seqs[seqNum-1].num_ << " " << seqs[seqNum-1].maxScore_
      //    << seqs[seqNum-1].data_ << endl;
#ifdef DEBUG
      cout << pOligos->getLastOligo() << endl;
      for (int i(0) ; i<seqSize;i++)
      {
        cout << i << " " << pOligos->getLastOligo()[i] << endl;
        for (int j(0);j<4;j++)
        {
          cout << baseNames[j] << " " <<
            scoreTable[i]
            [baseCodes[pOligos->getLastOligo()[i]]]
            [baseCodes[baseNames[j]]]
               << " " <<pScores->getScore(baseNames[j], i) << endl;
          assert
            (scoreTable[i][baseCodes[pOligos->getLastOligo()[i]]]
             [baseCodes[baseNames[j]]]
             ==pScores->getScore(baseNames[j], i));
        } // ~for j
      } // ~for i
#endif

      sequence = pOligos->getNextSequence(isValid);
    } // ~while

    if(dynamic_cast<ScoreSource*>(pOligosRaw)==NULL) delete pScoresRaw;
    delete pOligosRaw;

    //    fclose(pSeq);
  } // ~for i

  const uint numFrags(frags.size());

  const uint numSeqs(seqs.size());
  assert(numSeqs==seqNum);
  cerr << "Sequences of length " << seqSize << endl;
  cerr << "Total sequences: " << numSeqs+numIgnored << endl;
  cerr << "Will ignore " << numIgnored << endl;
  cerr << "Will align " << numSeqs << endl;



  // standard sort is horrendously slow for large genome
  // (4Mb = ~8 million frags) - replace with radix sort
  //  sort( frags.begin(), frags.end(), lessThanGenomeFrag );
  mySort( frags, seqSize );


  //  vector<Seq> fragStore;
  fragStore.reserve(numFrags);
  fragStore.push_back( Seq() );
  strncpy(fragStore[0].data, frags[0].seq_, seqSize);
  fragStore[0].data[seqSize]='\0';
  frags[0].seq_=fragStore[0].data;

  for (uint i(1);i<numFrags;i++)
  {
    if (strncmp(frags[i].seq_, fragStore.back().data, seqSize)!=0)
    {
      fragStore.push_back( Seq() );
      strncpy(fragStore.back().data, frags[i].seq_, seqSize);
      fragStore.back().data[seqSize]='\0';
    } // ~if
    frags[i].seq_=fragStore.back().data;
  } // ~if

  // Original genome sequences no longer required
  delete [] buf;
  delete [] bufR;

  vector<uint> commonPrefix(numFrags);
  for (uint i(0);i<numFrags-1;i++)
  {
    commonPrefix[i]=firstDiff(frags[i].seq_,frags[i+1].seq_,seqSize);
  } // ~for
  commonPrefix[numFrags-1]=0;





  QueryResults results( seqSize, numSeqs );

  // now do the matching
  PlusScore thisScore;
  uint lastPrefix(0);
  uint *baseIndex = new uint[seqSize+1];

  for (uint i(0);i<numFrags; i++)
  {

    for (uint k(1);k<=seqSize;k++)
    {
        // Can use realBaseCodes here as sequences with ambiguity codes
        // were excluded from addition to frags above.
        baseIndex[k]=(uint)realBaseCodes[(uint)frags[i].seq_[k-1]];
        assert(baseIndex[k] != nc);
    }
    if ((i%10000)==0) cerr << i << endl;

    for (uint k(lastPrefix+1);k<=commonPrefix[i];k++)
    {
      //      cout << k <<endl;
      results.cycle_[k]=results.cycle_[k-1];
      for (uint j(0);j<numSeqs;j++)
      {
	//	cout << j << " j " << endl;
	results.cycle_[k][j]+=seqs_.scores_[k][baseIndex[k]][j];
	//  thisScore+=seqs_.scores_[k][baseCodes[frags[i].seq_[k-1]]][j];
      }

    } // ~for k

    const int lastCompleteCycle(max(commonPrefix[i],lastPrefix));
    //    cout << i << " " << frags[i].seq_ << " " << lastPrefix << " " << commonPrefix[i] << endl << i << " " << frags[i+1].seq_ << endl;
    lastPrefix=commonPrefix[i];
    //    cout << i << "done common pref" << endl;

    for (uint j(0);j<numSeqs;j++)
    {
      thisScore=results.cycle_[lastCompleteCycle][j];

      for (uint k(lastCompleteCycle+1);
           (thisScore<=results.nextBestLess_[j])&&(k<=seqSize);k++)
      {
        thisScore+=seqs_.scores_[k][baseIndex[k]][j];
      } // ~for

      if (thisScore<results.nextBestLess_[j])
      {
        if (thisScore<results.bestLess_[j])
        {
          results.nextBestLess_[j]=results.bestLess_[j];
          results.bestLess_[j]=thisScore;
          results.bestPos_[j]=i;
          results.numBest_[j]=1;
        } // ~if
        else if (thisScore==results.bestLess_[j])
        {
          ++results.numBest_[j];
        } // ~else if
        else
        {
          results.nextBestLess_[j]=thisScore;
        } // ~else
      } // ~if
    } // ~for j
  } // ~for i


  // Now print out results

#ifndef SWITCHED_TO_EXPORT_OUTPUT
  std::string refStr;
  refStr.resize(seqSize);
#endif

  uint thisPos(0);

  for (uint j(0);j<numSeqs;j++)
  {
      const uint posVal(frags[results.bestPos_[j]].pos_);
      const casava::common::Match::Strand strandDir((posVal & reverseFlag)
                                                ? casava::common::Match::Reverse
                                                : casava::common::Match::Forward);
      const bool align_is_unique(results.numBest_[j] == 1);

      if (align_is_unique)
      {
          thisPos = (posVal & (~reverseFlag));

          if (strandDir == casava::common::Match::Forward)
          {
              ++thisPos;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
              refStr.assign(frags[results.bestPos_[j]].seq_, seqSize);
#endif
          } // ~if
          else
          {
              thisPos = genomeSize - seqSize - thisPos + 1;

#ifndef SWITCHED_TO_EXPORT_OUTPUT
              for (uint i(0);i<seqSize;i++)
                  refStr[i]
                      =reverseCharASCII
                      [(uint)frags[results.bestPos_[j]].seq_[seqSize-1-i]];
#endif
          }
      }



#ifdef SWITCHED_TO_EXPORT_OUTPUT
      if (align_is_unique)
      {
          const std::string contigStr; // empty
          casava::common::Match match(chrName, contigStr,
                                  thisPos, strandDir, true);

          MatchDescriptor matchDesc(frags[results.bestPos_[j]].seq_,
                                    seqs[j].data_);

          std::cout << casava::common::Alignment(sequenceVec[j],
                                             match,
                                             matchDesc.str(),
                                             (((int)seqs_.maxScore_[j])
                                              -((int)results.bestLess_[j])));
      } else {
          casava::common::Match match;
          match.setNoMatch();
          std::cout << casava::common::Alignment(sequenceVec[j], match);
      }

#else
      printf("%s %i %i", seqs[j].data_,
             ((int)seqs_.maxScore_[j])-((int)results.bestLess_[j]),
             results.numBest_[j]);

      if (align_is_unique)
      {
          printf(" %d %c %s %i\n",
                 thisPos,
                 ((strandDir == casava::common::Match::Forward) ? 'F' : 'R'),
                 refStr.c_str(),
                 //	     thisPos, thisDir, frags[seqs[j].bestPos_].seq_,
                 ((int)seqs_.maxScore_[j])-((int)results.nextBestLess_[j]));
      }
      else
      {
          printf("\n");
      }

#endif // SWITCHED_TO_EXPORT_OUTPUT



  } // ~for j


  // deallocate genome buffer
  //  delete [] buf;
  delete[] baseIndex;

  return 0;
} // ~phageAlign

/*****************************************************************************/

int main( int numArgs, char* args [] )
{
    int ret_val(0);

    try {
        ret_val = phageAlign(numArgs, args);
    }
    catch (const casava::common::ExceptionData& exception) {
        cerr << exception.getMessage() << std::endl;
        ret_val = 1;
    }
    catch (...) {
        cerr << "Caught unknown exception" << std::endl;
        ret_val = 1;
    }

    return ret_val;
}

/*****************************************************************************/