HashTable/HashTableGeneric.cpp

Go to the documentation of this file.

// #######################################################################

// SSAHA : Sequence Search and Alignment by Hashing Algorithm
// Version 3.2, released 1st March 2004
// Copyright (c) Genome Research 2002

// SSAHA is free software; you can redistribute it and/or modify 
// it under the terms of version 2 of the GNU General Public Licence
// as published by the Free Software Foundation.
 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public Licence for more details.
 
// You should have received a copy of the GNU General Public Licence
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
// or see the on-line version at http://www.gnu.org/copyleft/gpl.txt

// #######################################################################

// Module Name  : HashTableGeneric
// File Name    : HashTableGeneric.cpp
// Language     : C++
// Module Author: Anthony J. Cox (ac2@sanger.ac.uk)

// Description:

// Includes:

#include "HashTableGeneric.h"
#include "HashTable.h"
#include "HashTablePacked.h"
#include "HashTableTranslated.h"
#include "SequenceReader.h"
// #include "SequenceReaderCodon.h"
#include "GlobalDefinitions.h"
#include <fstream>
#include <algorithm>
#include <map>
#include <cmath>


// ### Function Definitions ###

// Name:
// Arguments:
// TYPE  NAME  IN/OUT COMMENT
// Returns: TYPE COMMENT

enum
{
  defaultMaxNumHits=10000
};


  // Function Name: Constructor
  // Arguments:     ostream&
  // Returns: N/A
HashTableGeneric::HashTableGeneric
( ostream& monitoringStream,
  const string& name,
  Allocator<PositionInHitList>& arrayAllocator ) : 
  isInitialized_( false ), 
  monitoringStream_( monitoringStream ),
  name_( name ),
  bitsPerSymbol_( gBaseBits ), // default: may be overwritten in subclass ctor
  maxNumHits_( defaultMaxNumHits ),
  hitListFormat_( gNotSpecified ),
  pArrayAllocator_
  ( arrayAllocator.clone(&pWordPositionInHitList_,
                         name+(string)".head",
                         monitoringStream_) ),
  pSequenceSizes_(NULL),
  pNameReader_(NULL)
  {
    monitoringStream_ << "constructing HashTableGeneric\n";
    if (name_=="")
    {
      char buf[50];
      static int numTables(0);
      sprintf(buf,"SSAHAHashTable-%d-%d",getpid(),++numTables);
      name_=(string)buf;
    }
    //    pArrayAllocator_->link(&pWordPositionInHitList_,name+(string)".head");
  } // constructor 

  // Function Name:
  // Arguments:
  // TYPE  NAME  IN/OUT COMMENT
  // Returns: TYPE COMMENT
  HashTableGeneric::~HashTableGeneric()
  { 
    monitoringStream_ << "destructing HashTableGeneric ...\n";
    if ( isInitialized_ )
    {
      monitoringStream_ << "... deallocating memory\n";
      //      delete [] pHitListForAllWords_; 
      //      pArrayAllocator->deallocate();
      delete pArrayAllocator_; // also deallocs pointer array
      // deletion of HitList is handled in HashTableView
      // pHitListAllocator->deallocate();
      // delete pHitListAllocator; // NB destructor order
      //      delete [] pWordPositionInHitList_;
      delete [] pSequenceSizes_;
      //      for ( vector<string*>::iterator i( sequenceNames_.begin() );
      //    i != sequenceNames_.end() ; ++i ) delete *i;
      delete pNameReader_;
    } // ~if 
    else 
    {
      monitoringStream_ 
      << "... never initialized, no memory deallocation required\n";
    } // ~else
  } // destructor

void HashTableGeneric::createHashTable
( SequenceReader& sequenceReader, int wordLength, int maxNumHits, 
  int stepLength )
{

    monitoringStream_ << "HashTable::createHashTable" << endl;

    //    bitsPerSymbol_ = sequenceReader.getBitsPerSymbol();

    if (bitsPerSymbol_ != sequenceReader.getBitsPerSymbol() )
    {
      monitoringStream_ << "Info: reader encodes using " 
                        << sequenceReader.getBitsPerSymbol() 
                        << " bits per symbol, hash table expects "
                        << bitsPerSymbol_ << " bits per symbol.\n";
      //  throw SSAHAException("Wrong encoder for hash table");
    } // ~if


    monitoringStream_ << "Info: hash table creation will proceed using "
                      << bitsPerSymbol_ << " bits per symbol.\n";

    sourceData_ = sequenceReader.getSourceDataType();

    if ( sourceData_ == gDNAData )
    {
      monitoringStream_ 
        << "Info: hash table will be created from DNA source data.\n";
    } // ~if
    else if ( sourceData_ == gProteinData )
    {
      monitoringStream_ 
        << "Info: hash table will be created from protein source data.\n";
    } // ~else
    else
    {
      monitoringStream_ 
      << "Info: source data for hash table is of unknown type.\n";
    } // ~else


    // check wordLength is sensible

    int bitsUsedPerWord( wordLength * bitsPerSymbol_ );
    
    if ( (bitsUsedPerWord <= 0 ) || ( bitsUsedPerWord > gBitsPerWord ) ) 
       // 4 base pairs storable per byte
    {
      monitoringStream_ << "Error: bits used per hash words (" 
                        << bitsUsedPerWord << ") outside valid range ( 1 to "
                        << gBitsPerWord << ").\n";
      throw SSAHAException("Invalid hash word length.");
    } // ~if

    monitoringStream_ << "Info: " << wordLength 
                      << " symbols will be stored per hash word.\n";

    // check stepLength is sensible. If not, set it equal to the word length
    if ( ( stepLength > wordLength ) || ( stepLength <= 0 ) ) 
    { 
      monitoringStream_ << "Info: requested step length for hashing ("
                        << stepLength << ") outside valid range ( 1 to "
                        << wordLength << "), using " << wordLength
                        << " instead.\n";
        
      stepLength = wordLength; 
    }

    wordLength_ = wordLength;
    stepLength_ = stepLength;
    maxNumHits_ = maxNumHits;

    AdapterFactory factory;

    SequenceAdapter* seq
      ( factory.create( bitsPerSymbol_, wordLength_, stepLength_ ) );


    // numDifferentWords_ is the number of possible different Words for a 
    // given word length ( equal to 2^(2*wordLength_) )
    numDifferentWords_ =  ( 1 << bitsUsedPerWord );

    setupPointerArray();

    
    // Now start to extract the sequence information from the file

    //    sequenceReader.rewind();

    int numSeqs(countWordsAndGetNames( sequenceReader, seq ));

    // Create the array of sequence sizes
    monitoringStream_ << "Allocating memory for pSequenceSizes_: "
                      << numSeqs << " sequences, "
                      << numSeqs * sizeof(SequenceOffset)
                      << " bytes total ..." << endl;

    pSequenceSizes_ = new SequenceOffset [ numSeqs ];

    if (!pSequenceSizes_) 
    {
      throw SSAHAException("Memory allocation failed!");
    } // ~if

    monitoringStream_ << "... memory allocation OK" << endl;

    computePointerArray();

    setupHitList();

    // Now go through the vector of WordSequences and fill in the table

    hashAllWords( sequenceReader, seq, numSeqs );

    cleanupTempData();

    delete seq;
    //   delete [] pHitsFoundSoFar_;

    // Next flag is only set once everything is OK
    isInitialized_ = true;

} // ~HashTable::createHashTable( int wordLength, int maxNumHits )


int HashTableGeneric::countWordsAndGetNames
( SequenceReader& sequenceReader, SequenceAdapter* seq )
{

    sequenceReader.rewind();

    //    vector<WordSequence> wordSeqs; 
    // const WordSequence dummy; // needed to get this to compile on LINUX ...
    //   wordSeqs.push_back(dummy);
    WordSequence thisSeq;
    int numSeqs(0);

    //    sequenceNames_.push_back(new string);

    NameReaderLocal* pReaderLocal = new NameReaderLocal;
    pNameReader_ = pReaderLocal;

    while( sequenceReader.getNextSequence
           ( thisSeq, wordLength_) != -1 ) 
    {
      numSeqs++;

      sequenceReader.getLastSequenceName( pReaderLocal->lastName() );
      //      sequenceReader.getLastSequenceName( *(sequenceNames_.back()) );

      seq->link( thisSeq );
      countWords(*seq);

      //      sequenceNames_.push_back( new string );
      //      wordSeqs.push_back(dummy);

    } // ~while

    //    delete sequenceNames_.back();
    //  sequenceNames_.pop_back();
    //    wordSeqs.pop_back();
    //    cout << "XXX: " << name_ << numSeqs << endl;
    return numSeqs;
} // ~HashTableGeneric::countWordsAndGetNames

void HashTableGeneric::hashAllWords
( SequenceReader& sequenceReader, SequenceAdapter* seq, int numSeqs )
{

  int numWords(0);
  sequenceReader.rewind();
  WordSequence thisSeq;
    
  // NB sequences are numbered 1...n not 0...n-1
  for ( unsigned int i(1); i <= numSeqs ; i++ )
  {
    if( sequenceReader.getNextSequence( thisSeq, wordLength_) == -1 ) 
    {
      throw SSAHAException
        ("Sequence source data changed during hash table creation!");
    }

    // fill in the word count
    numWords = (int) thisSeq.size();
    pSequenceSizes_[i-1] = ( numWords > 0 ) 
      ? ( (numWords-1) * wordLength_ ) + thisSeq.getNumBasesInLast()
      : 0;

    // do the hashing
    seq->link(thisSeq);
    hashWords( *seq,i );


  } // ~for thisSeq


} // ~HashTableGeneric::hashAllWords 





void HashTableGeneric::setupPointerArray( void )
{

    // Attempt to allocate memory for the list of table positions
    monitoringStream_ << "Allocating memory for pWordPositionInHitList_: "
                      << numDifferentWords_ << " words, "
                      << numDifferentWords_ * sizeof(PositionInHitList)
                      << " bytes total ..." << endl;

   
    pArrayAllocator_->allocateAndZero
      (numDifferentWords_);

    monitoringStream_ << "... done" << endl;

  } // ~HashTableGeneric::setupPointerArray( void )

void HashTableGeneric::computePointerArray( void )
{
    unsigned long numClipped( 0 );

    // Accumulate to get position of hit list in final array
    for ( unsigned int i(1) ; i < numDifferentWords_ ; i ++ )
    {
        pWordPositionInHitList_[i] += pWordPositionInHitList_[i-1];
    } // ~for

    wordsInDatabase_ = pWordPositionInHitList_[numDifferentWords_-1];

    monitoringStream_ << numClipped
    << " words will be discarded ("
    <<   (float)(100*numClipped)
        /(float)(wordsInDatabase_+numClipped)  
    << " \% of total)." << endl;  


} // ~HashTableGeneric::computePointerArray( void )


void HashTableGeneric::setupHitList( void )
{

    // Attempt to allocate memory for the table itself
    monitoringStream_ << "Allocating memory for pHitListForAllWords_: "
                      << wordsInDatabase_ << " words, "
                      << wordsInDatabase_ * getHitTypeSize()
                      << " bytes total ..." << endl;

    allocateHitList(wordsInDatabase_);

    monitoringStream_ << "... memory allocation OK" << endl;

    // Attempt to allocate temporary memory to store hits so far for each word
    monitoringStream_ << "Allocating memory for pHitsFoundSoFar_: "
                      << numDifferentWords_ << " words, "
                      << numDifferentWords_ * sizeof(unsigned short)
                      << " bytes total ..." << endl;



    //    Allocator<PositionInHitList>* tempAlloc = 
    tempAlloc_ = 
      pArrayAllocator_->clone(&pHitsFoundSoFar_, 
                                 name_+(string)".temp",
                                 monitoringStream_);
    // ... ensures uses same allocation method as for ptr array
    tempAlloc_->allocateAndZero( numDifferentWords_  );

    monitoringStream_ << "... done" << endl;

} // ~HashTableGeneric::setupHitList( void )

void HashTableGeneric::cleanupTempData( void )
{
    delete tempAlloc_; // automatically deallocs pHitsFoundSoFar
} // ~HashTableGeneric::cleanupTempData( void )




  // Accessor Functions
  // (NB all accessor functions should be 'const')

// Function Name: getSequenceName
// Arguments: string& (out), SequenceNumber& (in)
// Places the name of the seqNum'th sequence in the database in seqName
void HashTableGeneric::getSequenceName
( string& seqName, SequenceNumber seqNum ) const
{
  pNameReader_->getSequenceName( seqName, seqNum );
#ifdef XXX
  // NB sequenceNames_ are numbered 0 ... n-1
  // but we are numbering sequences 1 ... n
  if ((seqNum <= 0 )||( seqNum > sequenceNames_.size() ))
  {
    seqName = "No sequence name stored for sequence.";
    monitoringStream_ 
      << "Warning - no sequence name for requested sequence number ("
      << seqNum << ").\n";
    
    return; 
  } // ~if

  seqName = *(sequenceNames_[seqNum - 1]);
  return;
#endif
} // ~HashTable::getSequenceName


// Function Name: getSequenceName
// Arguments: SequenceNumber& (in)
// Returns: const char*
// Returns C style string of the seqNum'th seq in the database
// NB should be used with caution as the string only lasts as long
// as the hash table, you have to copy it if you want to keep it
// after the hash table's destructor is called for any reason.
const char* HashTableGeneric::getSequenceName
( SequenceNumber seqNum ) const
{
  return pNameReader_->getSequenceName( seqNum );
#ifdef XXX
  // NB sequenceNames_ are numbered 0 ... n-1
  // but we are numbering sequences 1 ... n
  if ((seqNum <= 0 )||( seqNum > sequenceNames_.size() ))
  {
    monitoringStream_ 
      << "Warning - no sequence name for requested sequence number ("
      << seqNum << ").\n";
    return NULL; 
  } // ~if

  return sequenceNames_[seqNum - 1]->c_str();
#endif
} // ~HashTable::getSequenceName

SequenceNumber HashTableGeneric::getNumSequences() const 
{ 
  return (pNameReader_ == NULL) ? 0 : pNameReader_->size();
} // ~SequenceNumber HashTableGeneric::getNumSequences() const 





  // Function Name: getSequenceSize
  // Arguments: string& (out), SequenceNumber& (in)
  // Returns the size of the seqNum'th sequence in the database.
  SequenceOffset HashTableGeneric::getSequenceSize( SequenceNumber seqNum ) const
  {
    // NB sequenceNames_ are numbered 0 ... n-1
    // but we are numbering sequences 1 ... n
    if ((seqNum <= 0 )||( seqNum > pNameReader_->size() ))
    {
      monitoringStream_ 
      << "Warning - no sequence size for requested sequence number ("
      << seqNum << ").\n";

      return 0; 
    } // ~if
    else return pSequenceSizes_[seqNum - 1];

  } // ~HashTable::getSequenceName


  HashTableGeneric* HashTableFactory::loadHashTable( const string& name, SourceReaderIndex* pReader )
  {

    monitoringStream_ << "loading hash table from file ...\n";

    ifstream nameFile( (name+(string)".name").c_str() );
    
    if ( nameFile.fail() )
    {
      monitoringStream_ << "Error: failed to open " 
                        << name << ".name, aborting load." << endl;
      throw SSAHAException("Could not open .name file ");
    } // ~if

    char firstChar = nameFile.peek();
    string temp;
    temp += firstChar;
    
    nameFile.close();

    HitListFormatType hitListFormat((HitListFormatType)atoi(temp.c_str()));

    HashTableGeneric* pTable(NULL);

    if (hitListFormat==gStandard)
    {
      monitoringStream_ << "... creating standard format hash table\n";
      pTable = new HashTable(monitoringStream_,name);
    } // ~if
    else if (hitListFormat==g32BitPacked)
    {
      monitoringStream_ << "... creating 32 bit packed format hash table\n";
      pTable = new HashTablePacked(monitoringStream_,name);
    } // ~else if
    else if (hitListFormat==g32BitPackedProtein)
    {
      monitoringStream_ 
        << "... creating 32 bit packed format protein hash table\n";
      pTable = new HashTablePackedProtein(monitoringStream_,name);
    } // ~else if
    else if (hitListFormat==gTranslated)
    {
      monitoringStream_ << "... creating translated format hash table\n";
      pTable = new HashTableTranslated(monitoringStream_,name);
    }
    else
    {
      monitoringStream_ 
        << "Error: type of hash table (" << hitListFormat 
        << " not recognised.\n";
      throw SSAHAException("Could not recognise hash table type");
    } // ~else

    assert (pTable!=NULL);
    
    loadHashTable(*pTable, pReader);

    assert (pTable->isInitialized());

    return pTable;

  }

  // Function Name: loadHashTable
  // Arguments: const string& (in)
  // Reads a pre-computed hash table into memory from a file
  void HashTableGeneric::loadHashTable
  ( SourceReaderIndex* pSourceReader )
  {

    monitoringStream_ << "HashTable::loadHashTable" << endl;

    if ( isInitialized_ == true ) 
    {
      monitoringStream_ << "Error: hash table cannot be loaded because class"
                        << " has already been initialized." << endl;
      return;
    } // ~if

    // HashTableComponents have a dummy name reader set up already
    // so don't need the next bit, hence the next if condition
    if (pNameReader_==NULL)
    {
      if (pSourceReader==NULL) 
      {
        monitoringStream_ << "Info: sequence names will be held in local memory" 
                          << endl;
        pNameReader_ = new NameReaderLocal;
      } // ~if
      else
      {
        monitoringStream_ << "Info: sequence names will be read from file" 
                          << endl;
        pNameReader_ = new NameReaderIndex(*pSourceReader);
      } // ~else
    } // ~if
    // Now load all the sequence names

    loadSequenceNames();

    unsigned long numSeqs = pNameReader_->size();

    //    numDifferentWords_ = 1 << ( bitsPerSymbol_ * wordLength_ );


    int bitsUsedPerWord( wordLength_ * bitsPerSymbol_ );
    numDifferentWords_ = ( 1 << bitsUsedPerWord );
    unsigned long headSize = numDifferentWords_ * sizeof(PositionInHitList);

    // Attempt to allocate memory for the list of table positions
    monitoringStream_ << "Allocating memory for pWordPositionInHitList_: "
                      << numDifferentWords_ << " words, "
                      << headSize
                      << " bytes total ..." << endl;

    
      //        pArrayAllocator_->allocate(headSize);
    pArrayAllocator_->load(numDifferentWords_);
      //        loadFromFile
      //    ( 
      //      name_+(string)".head",
      //     (char*)pWordPositionInHitList_, 
      //      headSize
      //    );


    wordsInDatabase_ = pWordPositionInHitList_[numDifferentWords_-1];
    //    cout << pWordPositionInHitList_[numDifferentWords-1] << endl;

    unsigned long bodySize 
      = wordsInDatabase_ * getHitTypeSize();

    monitoringStream_ << "Allocating memory for pHitListForAllWords_: "
                      << wordsInDatabase_ << " words, "
                      << bodySize
                      << " bytes total ..." << endl;

    loadHitList(wordsInDatabase_);

    // Now load all the sequence sizes

    // Don't need a .size file for HashTableTranslated ....
    if (numSeqs != 0 )
    {    
      monitoringStream_ << "Allocating memory for pSequenceSizes_: "
                        << numSeqs << " sequences, "
                        << numSeqs*sizeof(SequenceOffset)
                        << " bytes total ...\n";

      pSequenceSizes_ = new SequenceOffset [ numSeqs  ];
      
      if (!pSequenceSizes_) 
      {
        throw SSAHAException("Memory allocation failed!");
      } // ~if

      loadFromFile(name_+(string)".size", (char*)pSequenceSizes_, 
                   numSeqs * sizeof( SequenceOffset ),
                   monitoringStream_ );
    } // ~if

    isInitialized_ = true; // set once all loading has been done
    monitoringStream_ << "Hash table loading complete.\n";

  } // ~loadHashTable

void HashTableGeneric::loadSequenceNames( void )
{

    ifstreamSSAHA nameFile( ( name_ + ".name" ).c_str() );
    
    if ( nameFile.fail() )
    {
      monitoringStream_ << "Error: failed to open " 
                        << name_ << ".name, aborting load." << endl;
      throw SSAHAException("Could not open .name file");
    } // ~if

    string hashParams;

    getline( nameFile, hashParams );

    std::string::size_type pos(0), prev_pos(0);

    pos=hashParams.find_first_of(' ',pos);
    if (pos == std::string::npos)
    {
      monitoringStream_ << "Error: could not parse " 
                        << name_ << ".name, aborting load." << endl;
      throw SSAHAException("Could not parse .name file");
    } // ~if

    if ( hitListFormat_ 
         != (HitListFormatType)atoi( hashParams.substr
                                     (prev_pos, pos-prev_pos).c_str() ) )
    {
      monitoringStream_ 
      << "Error: saved data is not compatible with this hash table" << endl;
      throw SSAHAException("Tried to load incompatible data into hash table");
    } // ~if

    prev_pos = ++pos;

    pos=hashParams.find_first_of(' ',pos);
    wordLength_ 
      = atoi(hashParams.substr(prev_pos, pos-prev_pos).c_str());
    prev_pos = ++pos;

    pos=hashParams.find_first_of(' ',pos);
    stepLength_ 
      = atoi(hashParams.substr(prev_pos, pos-prev_pos).c_str());
    prev_pos = ++pos;

    pos=hashParams.find_first_of(' ',pos);
    bitsPerSymbol_ 
      = atoi(hashParams.substr(prev_pos, pos-prev_pos).c_str());
    prev_pos = ++pos;

    pos=hashParams.find_first_of(' ',pos);
    sourceData_ = (SourceDataType) 
      atoi(hashParams.substr(prev_pos, pos-prev_pos).c_str());
    prev_pos = ++pos;

    //    sequenceNames_.push_back( new string );
 


    //    while( getline( nameFile,*(sequenceNames_.back() ) ) )
    //   {
    //    sequenceNames_.push_back( new string );
    //   } // ~while

    //    sequenceNames_.pop_back();

    pNameReader_->loadSequenceNames( nameFile );

    monitoringStream_ << "Loaded file " << name_ << ".name (" 
                       << pNameReader_->size() 
                      << " sequence names).\n";

    monitoringStream_ << "Closing .name file\n";
    nameFile.close();

    monitoringStream_ << "Loaded in word length of " 
                      << wordLength_ << " bases.\n";

    monitoringStream_ << "Loaded in step length of " 
                      << stepLength_ << " bases.\n";

    monitoringStream_ << "Hash table format uses " << bitsPerSymbol_
                      << " bits per symbol.\n";

    monitoringStream_ << "Hash table was created from ";

    switch (sourceData_)
    {
    case gDNAData:
      monitoringStream_ << "DNA"; break;
    case gProteinData:
      monitoringStream_ << "protein"; break;
    default:
      monitoringStream_ << "unknown"; break;
    }

    monitoringStream_ << " data.\n";




} // ~HashTableGeneric::loadSequenceNames( void )


  // Function Name: saveHashTable
  // Arguments: const string& (in)
  // Saves a hash table to a file (for subsequent retrieval by loadHashTable)
  void HashTableGeneric::saveHashTable( void )
  {

    monitoringStream_ << "HashTableFactory::saveHashTable" << endl;

    if ( isInitialized_ == false ) 
    {
      monitoringStream_ << "Error: hash table cannot be saved because it was"
                        << " was not successfully created." << endl;
      throw SSAHAException("Tried to save uninitialized hash table");
    } // ~if

    // Save sequence names

    saveSequenceNames();


    //        int bitsUsedPerWord( wordLength_ * bitsPerSymbol_ );
    //  unsigned long numDifferentWords ( 1 << bitsUsedPerWord );

    pArrayAllocator_->save();
    //  pHitListAllocator_->save();

    saveHitList();

    //    if (sequenceNames_.size() != 0)
    if (pNameReader_->size() != 0)
    saveToFile( name_+(string)".size", 
                (char*)pSequenceSizes_, 
                pNameReader_->size() * sizeof( SequenceOffset ),
                monitoringStream_ ); 



  } // ~saveHashTable


void HashTableGeneric::saveSequenceNames( void )
{

  ofstreamSSAHA nameFile( ( name_ + ".name" ).c_str() );

  nameFile << hitListFormat_ << " " << wordLength_ << " " 
           << stepLength_ << " " << bitsPerSymbol_ 
           << " " << sourceData_ << " " << endl;
  

  pNameReader_->saveSequenceNames( nameFile );
    
    //    for ( vector<string*>::iterator i = sequenceNames_.begin() ; 
    //         i != sequenceNames_.end() ; i ++ )
    //  {
    //   nameFile << **i << endl;
    //   }

    if ( nameFile.fail() )
    {
      monitoringStream_ << "Error: failed to write " 
                        << name_ << ".name" << endl;
      throw SSAHAException("Problem saving .name file");
    } // ~if

    nameFile.close();

    monitoringStream_ << "Saved file " << name_ << ".name\n";


} // ~HashTableGeneric::saveSequenceNames( void )


void HashTableGeneric::printHashStats( void )
{
  cout << "Analyzing data in hash table...\n";
  
  static const double thresholds[] = 
  { 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 
    0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.03, 0.02, 0.01, 
    0.001, 0.0001, 0.00001,
    0 }; 


  if (!isInitialized()) 
  {
    cout << "Hash table has not been initialized.\n";
    return;
  } // ~if

  const char* pSymbols;

  cout << "Word length: " << wordLength_ << " bases.\n";
  cout << "Bits per symbol: " << bitsPerSymbol_ << endl;

  unsigned long numDifferentWords;

  //  if (    (pHashTable->getHitListFormat()==gTranslated)
  //      || (pHashTable->getHitListFormat()==g32BitPackedProtein) )
  if (bitsPerSymbol_==gResidueBits)
  {
    numDifferentWords= (unsigned long)pow((double)gNumCodonEncodings,
                                          (int)wordLength_);
    pSymbols = gResidueNames;
  } // ~if
  else
  {
    assert(bitsPerSymbol_==gBaseBits);
    numDifferentWords= 1<<(2*wordLength_);
    pSymbols = gBaseNames;
  } // ~else


  //  int bitsUsedPerWord( wordLength_ * bitsPerSymbol_ );
  //  unsigned long numDifferentWords ( 1 << bitsUsedPerWord );
  cout << "There are " << numDifferentWords 
       << " different possible words of this length.\n";
  cout << "There are " << wordsInDatabase_ 
       << " words in this hash table.\n"; 


  int numCommon(240);

  TopList mostCommon(numCommon, wordLength_, bitsPerSymbol_);

  map<int, int> histogram;
  vector<Word> neverOccurred;


  PositionInHitList last(0); 
  int numOccs;

  for ( unsigned long i(0); i < numDifferentWords ; ++i )
  {
    numOccs = pWordPositionInHitList_[i]-last;
    histogram[numOccs]++;
    if (numOccs==0) neverOccurred.push_back((Word)i);

    mostCommon.push_back( pair<int, Word>(numOccs, (Word)i ) ); 
    last = pWordPositionInHitList_[i];
  }

  cout << neverOccurred.size() << " words never appeared.\n";
  if (neverOccurred.size() < 200)
  {
    cout << "They are:\n";
    for (vector<Word>::iterator i(neverOccurred.begin());
         i!=neverOccurred.end();i++)
      //      cout << printWord( *i, wordLength_ ) << endl;
      cout 
        << PrintFromWord( *i, wordLength_, bitsPerSymbol_, pSymbols ) 
        << endl;
  } else cout << "Not showing them - more than 200!\n";

  cout << "The " << numCommon << " most common words are:\n";

  cout << mostCommon;

  unsigned long total(0), totalWordsUsed(0);
  int thresholdNum(0);
  unsigned int n(0);

  for ( map<int, int>::iterator i( histogram.begin()); 
        i != histogram.end() ; ++i )
  {
    total += (i->first*i->second);
    totalWordsUsed += i->second;
    if (n++<2000)
    {
      cout << i->second << " words occurred " << i->first 
           << " times in the hash table (" 
           << 100.0*total/(double)wordsInDatabase_
           << "/"
           << 100.0*totalWordsUsed/(double)numDifferentWords
           << " cumulative percent).\n";
    }




    while (   ( (double)(wordsInDatabase_-total)/(double)wordsInDatabase_ )
            < thresholds[thresholdNum] )
    {
      cout << "Threshold " << i->first << " retains "
           << 100.0*(1.0-thresholds[thresholdNum++])
           << "\% of words.\n";
    } // ~while
  } // ~for

  cout << "Total: " << total << "\n";

  cout << "... analysis complete.\n";

} // ~void HashTableGeneric::printHashStats( void )


WordSequenceShifted::WordSequenceShifted
( const WordSequence& thisSeq, const HashTableGeneric& hashTable ) :
  //hashTable_( hashTable ),
size_( ( thisSeq.size() - 2 ) * hashTable.getWordLength() 
       + thisSeq.getNumBasesInLast() + 1 )
// NB size_ is not the number of bases in the sequence, it is the number of 
// hash words that may be obtained from the sequence
{
  seqs_.push_back( thisSeq );

  int wordLength( hashTable.getWordLength() );

  for( int i(1) ; i < wordLength ; i++ ) 
  { 
    seqs_.push_back( seqs_.back() );
    shiftSequence( seqs_.back(), hashTable.getBitsPerSymbol(), wordLength );  
  } // ~for

} // ~constructor


SequenceAdapterWithOverlap::SequenceAdapterWithOverlap
( int bitsPerSymbol, int wordLength, int stepLength ) :
bitsPerSymbol_( bitsPerSymbol ), 
wordLength_( wordLength ), 
stepLength_( stepLength ),
SequenceAdapter() 
{

  maskLeft_  = new Word[ wordLength_ ];
  maskRight_ = new Word[ wordLength_ ];
  //  oops! TC 29.10.2001
  //  const Word allOnes = ( 1 << ( 2 * wordLength_ ) ) - 1;
  const Word allOnes = ( 1 << ( bitsPerSymbol_ * wordLength_ ) ) - 1;

  for ( int i(0) ; i < wordLength_ ; i++ )
  {
    maskLeft_[i] = ( 1 << ( bitsPerSymbol_ * (wordLength_ - i) ) ) - 1;
    maskRight_[i] = allOnes - maskLeft_[i];
    // next two lines ensure flag bit is masked out

    // BUG FIX FROM CHUCK DILLON - AWS 21/06/05
    //maskLeft_[i] ^= gCursedWord;
    //maskRight_[i] ^= gCursedWord;
    maskLeft_[i] &= allOnes;
    maskRight_[i] &= allOnes;
    // ~ BUG FIXED

  } // ~ for it
} // constructor

SequenceAdapterWithOverlap::~SequenceAdapterWithOverlap()
{
  delete [] maskLeft_;
  delete [] maskRight_;
} // destructor

void SequenceAdapterWithOverlap::link( const WordSequence& seq )
{
  SequenceAdapter::link( seq );
  // Added check for seq.size <=1 (used to cause size_ to be 
  // (unsigned long)-1, ie huge, hence seg fault) TC 29.10.2001
   size_
     = (seq.size()<=1)
     ? 0
     : ( ( ( (int) seq.size() - 2) * wordLength_ ) + seq.getNumBasesInLast() ) 
       / stepLength_;
}
Word SequenceAdapterWithOverlap::operator[]( size_type j )
{
  j *= stepLength_;
  const int wordPos( j / wordLength_ ), wordShift( j % wordLength_ );
  Word thisWord (0); 
  thisWord 
    |= ( (*pSeq_)[wordPos]   & maskLeft_[wordShift]  ) 
      << ( bitsPerSymbol_ * wordShift );
  thisWord
    |= ( (*pSeq_)[wordPos+1] & maskRight_[wordShift] ) 
      >> ( bitsPerSymbol_ * ( wordLength_ - wordShift ) );
    thisWord 
    |= (   gCursedWord 
           * (    ( ((*pSeq_)[wordPos]  & gCursedWord) != (Word)0 ) 
               || ( ((*pSeq_)[wordPos+1]& gCursedWord) != (Word)0 ) ) );
  return thisWord;
}

// Function definitions for NameReader and subclasses

const char* NameReaderLocal::getSequenceName
( SequenceNumber seqNum ) const
{
  //  assert(seqNum<size());
  if (seqNum <= size() )
  {
    return (*this)[seqNum-1].c_str();
  } // ~if
  return NULL;

} // ~const char* NameReaderLocal::getSequenceName

void NameReaderLocal::getSequenceName
( string& seqName, SequenceNumber seqNum ) const
{
  //  assert(seqNum<size());
  if (seqNum <= size() )
  {
    seqName = (*this)[seqNum-1];
    return;
  }
  seqName = "No sequence name stored for sequence.";
} // ~NameReaderLocal::getSequenceName



const char* NameReaderIndex::getSequenceName
( SequenceNumber seqNum ) const
{
  return reader_.extractName(seqNum);
} // ~NameReaderIndex::getSequenceName


void NameReaderIndex::getSequenceName
( string& seqName, SequenceNumber seqNum ) const
{
  seqName=(string)reader_.extractName(seqNum); 
} // ~NameReaderIndex::getSequenceName

void NameReaderLocal::saveSequenceNames( ostream& nameFile )
{
  for ( vector<string>::iterator i = begin() ; 
        i != end() ; i ++ )
  {
    nameFile << *i << endl;
  }
} // ~NameReaderLocal::saveSequenceNames( ostream& nameFile )

void NameReaderLocal::loadSequenceNames( istream& nameFile )
{
  while( getline( nameFile, lastName() ) );
  pop_back();
    //    while( getline( nameFile,*(sequenceNames_.back() ) ) )
    //   {
    //    sequenceNames_.push_back( new string );
    //   } // ~while

    //    sequenceNames_.pop_back();

  //    pNameReader->loadSequenceNames( nameFile );

  //    monitoringStream_ << "Loaded in "  
  //          << size() 
  //          << " sequence names).\n";


} // ~NameReaderLocal::loadSequenceNames( istream& inFile )

SequenceNumber NameReaderIndex::size( void ) const
{
  return reader_.size();
} // ~NameReaderIndex::size( void ) const


// End of file HashTableGeneric.cpp

---