libMems/Aligner.cpp

Go to the documentation of this file.

/*******************************************************************************
 * $Id: Aligner.cpp,v 1.47 2004/04/19 23:10:30 darling Exp $
 * This file is copyright 2002-2007 Aaron Darling and authors listed in the AUTHORS file.
 * Please see the file called COPYING for licensing, copying, and modification
 * Please see the file called COPYING for licensing details.
 * **************
 ******************************************************************************/

#include "libMems/Aligner.h"
#include "libMems/Islands.h"
#include "libMems/DNAFileSML.h"
#include "libMems/MuscleInterface.h"    // it's the default gapped aligner
#include "libGenome/gnRAWSource.h"
#include "libMems/DistanceMatrix.h"
#include "libMems/Files.h"

#include <map>
#include <fstream>      // for debugging
#include <sstream>
#include <stack>
#include <algorithm>
#include <limits>

using namespace std;
using namespace genome;
namespace mems {


boolean validateLCB( MatchList& lcb );
void validateRangeIntersections( vector< MatchList >& lcb_list  );
bool debug_shite = false;

boolean validateLCB( MatchList& lcb ){
        vector< Match* >::iterator lcb_iter = lcb.begin();
        if( lcb.size() == 0 )
                return true;
        uint seq_count = (*lcb_iter)->SeqCount();
        uint seqI = 0;
        boolean complain = false;
        for(; seqI < seq_count; seqI++ ){
                lcb_iter = lcb.begin();
                int64 prev_coord = 0;
                for(; lcb_iter != lcb.end(); lcb_iter++ ){
                        if( (*lcb_iter)->Start( seqI ) == NO_MATCH )
                                continue;
                        else if( prev_coord != 0 && (*lcb_iter)->Start( seqI ) < prev_coord ){
                                complain = true;
                        }
                        prev_coord = (*lcb_iter)->Start( seqI );
                }
        }
        return !complain;
}

void EliminateOverlaps( MatchList& ml ){
        if( ml.size() < 2 )
                return;
        vector< Match* > result_matches;
        uint seq_count = ml[0]->SeqCount();
        for( uint seqI = 0; seqI < seq_count; seqI++ ){
                SingleStartComparator<AbstractMatch> msc( seqI );
                sort( ml.begin(), ml.end(), msc );
                int64 matchI = 0;
                int64 nextI = 0;
                int64 deleted_count = 0;
                vector< Match* > new_matches;

                // scan forward to first defined match
                for(; matchI != ml.size(); matchI++ )
                        if( ml[ matchI ]->Start( seqI ) != NO_MATCH )
                                break;

                for(; matchI < ml.size(); matchI++ ){
                        if( ml[ matchI ] == NULL )
                                continue;
                        
                        for( nextI = matchI + 1; nextI < ml.size(); nextI++ ){
                                if( ml[ nextI ] == NULL )
                                        continue;

                                boolean deleted_matchI = false;
                                // check for overlaps
                                int64 startI = ml[ matchI ]->Start( seqI );
                                int64 lenI = ml[ matchI ]->Length();
                                int64 startJ = ml[ nextI ]->Start( seqI );
//                              int64 diff =  absolut( startJ ) - absolut( startI ) - lenI;
                                int64 diff =  absolut( startJ ) - absolut( startI ) - lenI;

                                if( diff < 0 ){
                                        diff = -diff;
                                        Match* new_match;
                                        // delete bases from the smaller match
//                                      if( ml[ nextI ]->Length() * ml[ nextI ]->Multiplicity() >= 
//                                              lenI * ml[ matchI ]->Multiplicity() ){
                                        if( ( ml[ nextI ]->Multiplicity() > ml[ matchI ]->Multiplicity() ) ||
                                                ( ml[ nextI ]->Multiplicity() == ml[ matchI ]->Multiplicity() && ml[ nextI ]->Length() > ml[ matchI ]->Length() ) ){
                                                // mem_iter is smaller
                                                new_match = ml[matchI]->Copy();
                                                // erase base pairs from new_match
                                                if( diff >= lenI ){
//                                                      cerr << "Deleting " << **mem_iter << " at the hands of\n" << **next_iter << endl;
                                                        ml[ matchI ]->Free();
                                                        ml[ matchI ] = NULL;
                                                        matchI--;
                                                        deleted_matchI = true;
                                                        deleted_count++;
                                                }else{
                                                        if( startI > 0 ){
                                                                ml[ matchI ]->CropEnd( diff );
                                                                new_match->CropStart( new_match->Length() - diff );
                                                        }else{
                                                                ml[ matchI ]->CropStart( diff );
                                                                new_match->CropEnd( new_match->Length() - diff );
                                                        }
                                                }
                                        }else{
                                                // match_iter is smaller
                                                new_match = ml[nextI]->Copy();
                                                // erase base pairs from new_match
                                                if( diff >= ml[ nextI ]->Length() ){
//                                                      cerr << "Deleting " << **next_iter << " at the hands of\n" << **mem_iter << endl;
                                                        ml[ nextI ]->Free();
                                                        ml[ nextI ] = NULL;
                                                        deleted_count++;
                                                }else{
                                                        if( startJ > 0 ){
                                                                ml[ nextI ]->CropStart( diff );
                                                                new_match->CropEnd( new_match->Length() - diff );
                                                        }else{
                                                                ml[ nextI ]->CropEnd( diff );
                                                                new_match->CropStart( new_match->Length() - diff );
                                                        }
                                                }

                                        }
                                        new_match->SetStart( seqI, 0 );
                                        if( new_match->Multiplicity() > 1 && new_match->Length() > 0 )
                                                new_matches.push_back( new_match );
                                        else
                                        {
                                                new_match->Free();
                                                new_match = NULL;
                                        }
                                        if( deleted_matchI )
                                                break;
                                }else
                                        break;  // there are no more overlaps
                        }
//                      if( nextI > 1 )
//                              cerr << "There were " << nextI << " overlaps\n";
//                      if( nextI > config_value_2 )
//                              __asm(nop);
                }

                if( deleted_count > 0 ){
                        result_matches.reserve( ml.size() - deleted_count );
                        for( int64 copyI = 0; copyI < ml.size(); copyI++ ){
                                if( ml[ copyI ] != NULL )
                                        result_matches.push_back( ml[ copyI ] );
                        }
                        ml.clear();
                        ml.insert( ml.end(), result_matches.begin(), result_matches.end() );
                }
                ml.insert( ml.end(), new_matches.begin(), new_matches.end() );
                result_matches.clear();
                new_matches.clear();
        }
                
}


const gnSeqI default_min_r_gap_size = 200;
Aligner::Aligner( uint seq_count ) :
debug(false),
seq_count(seq_count),
min_recursive_gap_length(default_min_r_gap_size),
collinear_genomes(false),
gal(&(MuscleInterface::getMuscleInterface())),
permutation_weight(-1),
cur_min_coverage(-1),
max_extension_iters(4)
{}

Aligner::Aligner( const Aligner& al ) :
//gap_mh( al.gap_mh ),
nway_mh( al.nway_mh ),
seq_count( al.seq_count ),
debug( al.debug),
LCB_minimum_density( al.LCB_minimum_density),
LCB_minimum_range( al.LCB_minimum_range ),
cur_min_coverage( al.cur_min_coverage),
min_recursive_gap_length( al.min_recursive_gap_length ),
collinear_genomes( al.collinear_genomes ),
gal( al.gal ),
permutation_weight( al.permutation_weight ),
permutation_filename( al.permutation_filename ),
max_extension_iters( al.max_extension_iters )
{}

Aligner& Aligner::operator=( const Aligner& al )
{
        gap_mh = al.gap_mh;
        nway_mh = al.nway_mh;
        seq_count = al.seq_count;
        debug = al.debug;
        
        LCB_minimum_density = al.LCB_minimum_density;
        LCB_minimum_range = al.LCB_minimum_range;
        
        cur_min_coverage = al.cur_min_coverage;
        min_recursive_gap_length = al.min_recursive_gap_length;
        collinear_genomes = al.collinear_genomes;

        gal = al.gal;

        permutation_weight = al.permutation_weight;
        permutation_filename = al.permutation_filename;

        max_extension_iters = al.max_extension_iters;

        return *this;
}

void Aligner::SetMinRecursionGapLength( gnSeqI min_r_gap ) {
        min_recursive_gap_length = min_r_gap;
}

void Aligner::SetGappedAligner( GappedAligner& gal ){
        this->gal = &(gal);
}

void Aligner::SetMaxGappedAlignmentLength( gnSeqI len ){
        gal->SetMaxAlignmentLength( len );
}


/* returns true if all labels between start_label and end_label are contained in the no_match_labels set */
void scanLabels( set< uint >& no_match_labels, uint& start_label, boolean forward ){
        uint labelI;
        // scan no_match_labels for consecutive labels starting at start_label until one is missing
        if( forward ){
                for( labelI = start_label + 1; ; labelI++){
                        set< uint >::iterator  label_iter = no_match_labels.find( labelI );
                        if( label_iter == no_match_labels.end() ){
                                start_label = labelI - 1;
                                break;
                        }
                }
        }else{
                for( labelI = start_label; labelI > 0; labelI--){
                        set< uint >::iterator  label_iter = no_match_labels.find( labelI - 1 );
                        if( label_iter == no_match_labels.end() ){
                                start_label = labelI;
                                break;
                        }
                }
        }
}

boolean checkCollinearity( Match* m1, Match* m2 ){
        for( uint seqI = 0; seqI < m1->SeqCount(); seqI++ ){
                if( m1->Start( seqI ) == NO_MATCH ||
                        m2->Start( seqI ) == NO_MATCH )
                        continue;
                if((( m1->Start( seqI ) > 0 &&
                        m2->Start( seqI ) > 0 ) ||
                        (m1->Start( seqI ) < 0 &&
                        m2->Start( seqI ) < 0 )) &&
                        m1->Start( seqI ) <= m2->Start( seqI ) )
                        continue;
                return false;
        }
        return true;
}

void scanFit( list< LabeledMem >& pair_list, list< LabeledMem >::iterator& list_iter, Match* new_match, uint sort_seq ){

        list< LabeledMem >::iterator cur_iter = list_iter;
        list< LabeledMem >::iterator last_iter = list_iter;
//      int64 initial_start = absolut( list_iter->mem->Start( sort_seq ) );
        int64 initial_start = absolut( list_iter->mem->Start( sort_seq ) );

        uint match_count = 0;
        for(; last_iter != pair_list.end(); ++last_iter ){
                if( last_iter->mem->Start( sort_seq ) == NO_MATCH ){
                        ++match_count;
                        continue;
                }
//              if( absolut( last_iter->mem->Start( sort_seq ) ) < initial_start ||
//                      absolut( last_iter->mem->Start( sort_seq ) ) > new_match->Start( sort_seq ) )
                if( absolut( last_iter->mem->Start( sort_seq ) ) < initial_start ||
                        absolut( last_iter->mem->Start( sort_seq ) ) > new_match->Start( sort_seq ) )
                        break;
                ++match_count;
        }
        vector< vector< int > > score_vector;
        score_vector.reserve( new_match->SeqCount() - sort_seq - 1 );
        for( uint seqI = sort_seq + 1; seqI < new_match->SeqCount(); ++seqI ){
                vector< int > sv;
                score_vector.push_back( sv );
                score_vector[ score_vector.size() - 1 ].reserve( match_count );
        }
        uint matchI = 0;
        for(; cur_iter != last_iter; ++cur_iter ){
                
                for( uint seqI = sort_seq + 1; seqI < new_match->SeqCount(); ++seqI ){
                        int64 p_start = cur_iter->mem->Start( seqI );
                        int64 m_start = new_match->Start( seqI );
                        p_start = p_start < 0 ? -p_start : p_start;
                        m_start = m_start < 0 ? -m_start : m_start;
                        if( m_start == NO_MATCH ){
                                score_vector[ seqI - sort_seq - 1 ].push_back( 0 );
                        }else if( p_start == NO_MATCH ){
                                score_vector[ seqI - sort_seq - 1 ].push_back( 0 );
                        }else if( p_start < m_start ){
                                score_vector[ seqI - sort_seq - 1 ].push_back( 1 );
                        }else
                                score_vector[ seqI - sort_seq - 1 ].push_back( -1 );
                }
        }
        vector< int > scores;
        scores.reserve( match_count );
        for( matchI = match_count; matchI > 0; matchI-- )
                scores.push_back( 0 );
        for( uint seqI = 0; seqI < new_match->SeqCount() - sort_seq - 1; ++seqI ){
                boolean redefined = false;
                for( matchI = match_count; matchI > 0; matchI-- ){
                        if( !redefined ){
                                if( score_vector[ seqI ][ matchI - 1 ] >= 0 ){
                                        if( score_vector[ seqI ][ matchI - 1 ] == 1 )
                                                redefined = true;
                                        ++scores[ matchI - 1 ];
                                }
                        }else{
                                if( score_vector[ seqI ][ matchI - 1 ] == -1 )
                                        redefined = false;
                        }
                }
        }
        // find the first highest scoring match
        cur_iter = list_iter;
        int max_score = 0;
        for( matchI = 0; matchI < match_count; ++matchI ){
                if( scores[ matchI ] > max_score ){
                        max_score = scores[ matchI ];
                        list_iter = cur_iter;
                }
                ++cur_iter;
        }
}

void AaronsLCB( MatchList& mlist, set<uint>& breakpoints ){
        breakpoints.clear(); // make sure this is empty
        if( mlist.size() == 0 )
                return;
        // can only look for breakpoints if there is more than one match!!
        if( mlist.size() == 1 ){
                breakpoints.insert( 0 );
                return;
        }
        uint seq_count = mlist[0]->SeqCount();

        SingleStartComparator<AbstractMatch> msc( 0 );
        sort( mlist.begin(), mlist.end(), msc );
        vector<Match*>::iterator mem_iter = mlist.begin();
        list<LabeledMem> pair_list;
        
        map<uint, Match*> debug_label_map;
        boolean debugging = false;
        
        
        list< PlacementMatch > placement_list;
        
        for(; mem_iter != mlist.end(); ++mem_iter ){
                if( (*mem_iter)->Start( 0 ) != NO_MATCH ){              
                        // add this one to the list.
                        LabeledMem lm;
                        lm.mem = *mem_iter;
                        lm.label = 0;
                        pair_list.push_back( lm );
                }else{
                        PlacementMatch pm;
                        pm.mem = *mem_iter;
                        pm.iter = pair_list.end();
                        placement_list.push_back( pm );
                }
        }
        LabeledMemComparator lmc( 0 );
        pair_list.sort( lmc );
        list< LabeledMem >::iterator pair_iter = pair_list.begin();
        for(; pair_iter != pair_list.end(); ++pair_iter ){
                PlacementMatch pm;
                pm.mem = pair_iter->mem;
                pm.iter = pair_iter;
                placement_list.push_back( pm );
        }
        
        // place all the subset matches from each sequence in the correct place in pair_list.
        for( uint seqI = 1; seqI < seq_count; ++seqI ){
                PlacementMatchComparator pmc( seqI );
                placement_list.sort( pmc );
                list< PlacementMatch >::iterator placement_prev;
                list< PlacementMatch >::iterator placement_iter = placement_list.begin();
                if( placement_iter->iter == pair_list.end() &&
                        placement_iter->mem->Start( seqI ) != NO_MATCH ){
                        LabeledMem lm;
                        lm.mem = placement_iter->mem;
                        lm.label = 0;
                        pair_list.insert( pair_list.begin(), lm );
                        placement_iter->iter = pair_list.begin();
                }

                for( ++placement_iter; placement_iter != placement_list.end(); ++placement_iter ){
                        placement_prev = placement_iter;
                        placement_prev--;
                        
                        if( placement_iter->iter != pair_list.end() )
                                continue;
                        
                        if( placement_iter->mem->Start( seqI ) == NO_MATCH )
                                continue;
                        
                        list< LabeledMem >::iterator insert_iter = placement_prev->iter;
                        if( insert_iter == pair_list.end() || placement_prev->mem->Start( seqI ) == NO_MATCH )
                                insert_iter = pair_list.begin();
                        else{
                                if( insert_iter->mem->Start( seqI ) < 0 ){
                                        // invert if necessary and insert before
                                        if( placement_iter->mem->Start( seqI ) > 0 )
                                                placement_iter->mem->Invert();
                                        if( !checkCollinearity( placement_iter->mem, insert_iter->mem ) ){
                                                placement_iter->mem->Invert();
                                                scanFit( pair_list, insert_iter, placement_iter->mem, seqI );
                                                ++insert_iter;
                                        }
                                }else{
                                        // insert in the earliest place this match fits with surrounding matches
                                        scanFit( pair_list, insert_iter, placement_iter->mem, seqI );
                                        ++insert_iter;
                                }
                        }
                        
                        LabeledMem lm;
                        lm.mem = placement_iter->mem;
                        lm.label = 0;
                        pair_list.insert( insert_iter, lm );
                        placement_iter->iter = insert_iter;
                        placement_iter->iter--;
                }
        }
        boolean debug_labels = false;
        ofstream debug_label_file;
        if( debug_labels )
                debug_label_file.open( "label_debug.txt" );
        // number the LabeledMems in the pair_list
        uint cur_label = 0;
        mlist.clear();
        vector< LabeledMem > pair_vec;
        pair_vec.reserve( pair_list.size() );
        mlist.reserve( pair_list.size() );
        for( pair_iter = pair_list.begin(); pair_iter != pair_list.end(); ++pair_iter ){
                pair_iter->label = cur_label++;
                mlist.push_back( pair_iter->mem );
                pair_vec.push_back( *pair_iter );
                if( debug_labels ){
                        debug_label_map.insert( map<uint, Match*>::value_type( pair_iter->label, pair_iter->mem ) );
                        debug_label_file << pair_iter->label << '\t' << (*pair_iter->mem) << endl;
                }
        }
        if( debug_labels )
                debug_label_file.close();
        
        breakpoints.clear();
        pair_list.clear();
        vector< LabeledMem >::iterator pair_vec_iter;
        for( uint seqI = 1; seqI < seq_count; seqI++ ){
                // sort the list on the current genome
                LabeledMemComparator lmc( seqI );
                sort( pair_vec.begin(), pair_vec.end(), lmc );
                set< uint > no_match_labels;

                // debugging code
/*              stringstream debug_filename;
                debug_filename << "label_sort_" << seqI << ".txt";
                ofstream debug_file( debug_filename.str().c_str() );
                for( uint pairI = 0; pairI < pair_vec.size(); pairI++ ){
                        debug_file << pair_vec[ pairI ].label << *pair_vec[ pairI ].mem << endl;
                }
                debug_file.close();
*/              // end debugging code
                
                pair_vec_iter = pair_vec.begin();
                uint block_start = pair_vec_iter->label;
                uint break_label = 0;
                for( ++pair_vec_iter; pair_vec_iter != pair_vec.end(); ++pair_vec_iter ){
                        vector<LabeledMem>::iterator pair_prev = pair_vec_iter;
                        pair_prev--;
                        break_label = 0;
                        uint scan_label = 0;
                        if( pair_prev->mem->Start( seqI ) == NO_MATCH ){
                                no_match_labels.insert( set< uint >::value_type( pair_prev->label ) );
                                // get the correct block start
                                if( pair_vec_iter->mem->Start( seqI ) < 0 ){
                                        block_start = pair_vec_iter->label;
                                        scanLabels( no_match_labels, block_start, true );
                                }else if( pair_vec_iter->mem->Start( seqI ) > 0 ){
                                        block_start = pair_vec_iter->label;
                                        scanLabels( no_match_labels, block_start, false );
                                }
                                
                                continue;
                        }

                        if( pair_prev->mem->Start( seqI ) < 0 ){
                                // this block would break at its start
                                break_label = block_start;
                        }else{
                                // this block would break at its end
                                break_label = pair_prev->label;
                                scanLabels( no_match_labels, break_label, true );
                        }
                        if( pair_vec_iter->mem->Start( seqI ) < 0 ){
                                // scan forward to the beginning of new block
                                scan_label = pair_vec_iter->label;
                                scanLabels( no_match_labels, scan_label, true );
                        }else{
                                // scan back to the beginning of new block
                                scan_label = pair_vec_iter->label;
                                scanLabels( no_match_labels, scan_label, false );
                        }

                        if( pair_vec_iter->mem->Start( seqI ) < 0 &&
                                pair_prev->mem->Start( seqI ) < 0 ){
                                if( scan_label + 1 == pair_prev->label )
                                        continue;
                                if( debugging ){
                                        map< uint, Match* >::const_iterator debug_iter = debug_label_map.find( pair_vec_iter->label );
                                        while( debug_iter->first <= pair_prev->label ){
                                                cout << debug_iter->first << '\t' << *(debug_iter->second) << endl;
                                                ++debug_iter;
                                        }
                                }
                        }else
                        if( pair_vec_iter->mem->Start( seqI ) > 0 &&
                                pair_prev->mem->Start( seqI ) > 0 ){
                                
                                if( scan_label - 1 == pair_prev->label )
                                        continue;
                                if( debugging ){
                                        map< uint, Match* >::const_iterator debug_iter = debug_label_map.find( pair_prev->label );
                                        while( debug_iter->first <= pair_vec_iter->label ){
                                                cout << debug_iter->first << '\t' << *(debug_iter->second) << endl;
                                                ++debug_iter;
                                        }
                                }
                        }
                        // check if the missing matches are in the set of non-matches

                        // since it didn't meet any of the above
                        // criteria it's a breakpoint.  insert the label of the end of the current block
                        // note that if it's a reverse complement block, the end label is really the start label
                        breakpoints.insert( break_label );
                        block_start = scan_label;
                }

                // insert the correct block ending
                if( pair_vec_iter != pair_vec.begin() ){
                        pair_vec_iter--;
                        
                        if( pair_vec_iter->mem->Start( seqI ) < 0 ){
                                break_label = block_start;
                        }else{
                                break_label = pair_vec_iter->label;
                                scanLabels( no_match_labels, break_label, true );
                        }
                        breakpoints.insert( break_label );
                }
        }
}

void Aligner::SetPermutationOutput( std::string& permutation_filename, int64 permutation_weight )
{
        this->permutation_filename = permutation_filename;
        this->permutation_weight = permutation_weight;
}


void GetLCBCoverage( MatchList& lcb, uint64& coverage ){
        vector< Match* >::iterator match_iter = lcb.begin();
        coverage = 0;
        bool debug = true;
        for( ; match_iter != lcb.end(); ++match_iter ){
                coverage += (*match_iter)->Length() * (*match_iter)->Multiplicity();

                // if we have sequence information then
                // subtract the coverage for any position that contains an N
                if( lcb.seq_table.size() > 0 )
                {
                        for( uint seqI = 0; seqI < (*match_iter)->SeqCount(); ++seqI )
                        {
                                gnSeqI lend = absolut((*match_iter)->Start(seqI));
                                gnSeqI length = (*match_iter)->Length();
                                if( lend == 0 )
                                        continue;
                                string match_seq = lcb.seq_table[seqI]->ToString(length, lend);
                                for( size_t s = 0; s < match_seq.size(); ++s )
                                        if( match_seq[s] == 'n' || match_seq[s] == 'N' )
                                                if( (*match_iter)->Start(seqI) > 0 )
                                                        coverage--;
                        }
                }
        }
}


void computeLCBAdjacencies_v2( vector<MatchList>& lcb_list, vector< int64 >& weights, vector< LCB >& adjacencies ){
        IntervalList iv_list;
        for( uint lcbI = 0; lcbI < lcb_list.size(); ++lcbI ){
                vector<AbstractMatch*> asdf;
                asdf.push_back( lcb_list[ lcbI ].front() );
                if( lcb_list[lcbI].size() > 1 )
                        asdf.push_back( lcb_list[ lcbI ].back() );
                Interval iv( asdf.begin(), asdf.end() );
                iv_list.push_back( iv );
        }
        computeLCBAdjacencies_v2( iv_list, weights, adjacencies );
}

const uint NO_ADJACENCY = (std::numeric_limits<uint>::max)();

void computeLCBAdjacencies_v2( IntervalList& iv_list, vector< int64 >& weights, vector< LCB >& adjacencies ){
        adjacencies.clear(); // start with no LCB adjacencies
        if( iv_list.size() == 0 )
                return; // there aren't any LCBs so there aren't any adjacencies!

        uint seq_count = iv_list[0].SeqCount();
        uint seqI;
        uint lcbI;
        adjacencies.resize(iv_list.size());
        for( lcbI = 0; lcbI < iv_list.size(); ++lcbI ){
                LCB& lcb = adjacencies[lcbI];
                lcb.left_end.resize(seq_count);
                lcb.right_end.resize(seq_count);
                lcb.left_adjacency.resize(seq_count);
                lcb.right_adjacency.resize(seq_count);
                for( seqI = 0; seqI < seq_count; seqI++ ){
                        // support "ragged edges" on the ends of LCBs
                        int64 leftI = iv_list[lcbI].LeftEnd(seqI);
                        int64 rightI = NO_MATCH;
                        if( leftI != NO_MATCH )
                        {
                                leftI = iv_list[lcbI].Orientation(seqI) == AbstractMatch::forward ? leftI : -leftI;
                                rightI = iv_list[lcbI].RightEnd(seqI)+1;
                                rightI = iv_list[lcbI].Orientation(seqI) == AbstractMatch::forward ? rightI : -rightI;
                        }

                        lcb.left_end[seqI] = leftI;
                        lcb.right_end[seqI] = rightI;
                        lcb.left_adjacency[seqI] = NO_ADJACENCY;
                        lcb.right_adjacency[seqI] = NO_ADJACENCY;
                }
                lcb.lcb_id = lcbI;
                lcb.weight = weights[ lcbI ];
                lcb.to_be_deleted = false;
        }

        for( seqI = 0; seqI < seq_count; seqI++ ){
                LCBLeftComparator llc( seqI );
                sort( adjacencies.begin(), adjacencies.end(), llc );
                for( lcbI = 1; lcbI + 1 < iv_list.size(); lcbI++ ){
                        adjacencies[ lcbI ].left_adjacency[ seqI ] = adjacencies[ lcbI - 1 ].lcb_id;
                        adjacencies[ lcbI ].right_adjacency[ seqI ] = adjacencies[ lcbI + 1 ].lcb_id;
                }
                if( lcbI == iv_list.size() )
                        lcbI--; // need to decrement when there is only a single LCB

                // set first and last lcb adjacencies to -1
                adjacencies[ 0 ].left_adjacency[ seqI ] = NO_ADJACENCY;
                adjacencies[ lcbI ].right_adjacency[ seqI ] = NO_ADJACENCY;
                if( lcbI > 0 ){
                        adjacencies[ 0 ].right_adjacency[ seqI ] = adjacencies[ 1 ].lcb_id;
                        adjacencies[ lcbI ].left_adjacency[ seqI ] = adjacencies[ lcbI - 1 ].lcb_id;
                }
        }
        LCBIDComparator lic;
        sort( adjacencies.begin(), adjacencies.end(), lic );
        
}


void scanLeft( int& left_recurseI, vector< LCB >& adjacencies, int min_weight, int seqI ){
        while( left_recurseI != -1 && adjacencies[ left_recurseI ].weight < min_weight )
                left_recurseI = adjacencies[ left_recurseI ].left_adjacency[ seqI ];
}
void scanRight( int& right_recurseI, vector< LCB >& adjacencies, int min_weight, int seqI ){
        while( right_recurseI != -1 && adjacencies[ right_recurseI ].weight < min_weight )
                right_recurseI = adjacencies[ right_recurseI ].right_adjacency[ seqI ];
}



void CreateGapSearchList( vector< LCB >& adjacencies, const vector< gnSequence* >& seq_table, vector< vector< int64 > >& iv_regions, boolean entire_genome ) 
{
        iv_regions.clear();
        if( adjacencies.size() == 0 )
                return;         // there aren't any intervening LCB regions!
        if( adjacencies.size() == 1 && !entire_genome )
                return;         // there aren't any interveniing LCB regions in the local area
        boolean debug_lcb_extension = false;    
        const uint seq_count = seq_table.size();

        uint seqI = 0;
        int lcbI = 0;
        iv_regions = vector< vector< int64 > >( seq_count );

        // extract a gnSequence containing only the intervening regions
        for( seqI = 0; seqI < seq_count; seqI++ ){

                // find the first LCB in this sequence
                for( lcbI = 0; lcbI < adjacencies.size(); lcbI++ ){
                        if( adjacencies[ lcbI ].left_adjacency[ seqI ] == -1 )
                                break;
                }
                // start concatenating the intervening regions
                // scan right
                int right_recurseI = lcbI;
                lcbI = -1;
                if( !entire_genome && right_recurseI != -1 ){
                        lcbI = right_recurseI;
                        right_recurseI = adjacencies[ lcbI ].right_adjacency[ seqI ];
                }
                gnSeqI seq_len = 0;
                while( (lcbI != -1 || right_recurseI != -1 ) && right_recurseI < (int)adjacencies.size() ){
                        int64 l_end = lcbI == -1 ? 1 : adjacencies[ lcbI ].right_end[ seqI ];
                        int64 r_end = right_recurseI == -1 ? seq_table[ seqI ]->length() : adjacencies[ right_recurseI ].left_end[ seqI ];

                        // break out if outside the last LCB and not searching the entire genome
                        if( !entire_genome && right_recurseI == -1 )
                                break;

                        l_end = absolut( l_end );
                        r_end = absolut( r_end );
                        
                        if( l_end > r_end && !( r_end + 1 == l_end && right_recurseI == -1 ) ){
                                std::cerr << "Overlapping LCBs.  lcbI " << lcbI << " right_recurseI " << right_recurseI << endl;
                                std::cerr << "lend: " << l_end << " rend: " << r_end << endl;
                                l_end = r_end;
                                
                        }
                        
                        lcbI = right_recurseI;
                        if( right_recurseI != -1 )
                                right_recurseI = adjacencies[ right_recurseI ].right_adjacency[ seqI ];
                        if( r_end + 1 == l_end && right_recurseI == -1 )
                                continue;       // we're at the right end and there's nothing to add
                        seq_len += r_end - l_end;
                        iv_regions[ seqI ].push_back( l_end );
                        iv_regions[ seqI ].push_back( r_end );
                }
                if( debug_lcb_extension )
                        std::cerr << "seqI " << seqI << " seq_len: " << seq_len << endl;
        }

}

void SearchLCBGaps( MatchList& new_matches, const std::vector< std::vector< int64 > >& iv_regions, MaskedMemHash& nway_mh ) {
        if( iv_regions.size() == 0 )
                return;         // there aren't any intervening LCB regions!
        size_t sI = 0;
        for( ; sI < iv_regions.size(); sI++ )
                if( iv_regions[sI].size() > 0 )
                        break;
        if( sI == iv_regions.size() )
                return;         // there aren't any intervening LCB regions!

        boolean debug_lcb_extension = false;    
        const uint seq_count = new_matches.seq_table.size();
        uint seqI = 0;
        int lcbI = 0;
        MatchList gap_list;
        gap_list.seq_table = vector< gnSequence* >( seq_count );        
        gap_list.sml_table = vector< SortedMerList* >( seq_count );

        // extract a gnSequence containing only the intervening regions
        for( seqI = 0; seqI < seq_count; seqI++ ){
                gap_list.seq_table[ seqI ] = new gnSequence();
                gap_list.sml_table[ seqI ] = new DNAMemorySML();
                gnSeqI seq_len = 0;
                for( size_t ivI = 0; ivI < iv_regions[seqI].size(); ivI += 2 )
                {
                        int64 l_end = iv_regions[seqI][ivI];
                        int64 r_end = iv_regions[seqI][ivI+1];
                        try{
                        if( debug_lcb_extension )
                                cerr << "Adding " << seqI << "\t" << l_end << "\t" << r_end << "\t(" << r_end - l_end << " bp)" << endl;
                        gap_list.seq_table[ seqI ]->append( new_matches.seq_table[ seqI ]->ToString(r_end - l_end, l_end ) );
//                      gap_list.seq_table[ seqI ]->append( new_matches.seq_table[ seqI ]->subseq( l_end, r_end - l_end ) );
                        }catch(...){
                                cout << "";
                        }
                        seq_len += r_end - l_end;
                }
                if( debug_lcb_extension )
                        cerr << "seqI " << seqI << " seq_len: " << seq_len << endl;
        }
        //
        // search for MUMs in the intervening sequence regions
        //

        // calculate potential mer sizes for searches
        gnSeqI total_iv_length = 0;
        for( seqI = 0; seqI < seq_count; seqI++ ){
                total_iv_length += gap_list.seq_table[ seqI ]->length();
/*              cerr << "seqI: " << seqI << " length: " << gap_list.seq_table[ seqI ]->length();
                cerr << "\n";
*/
        }
        total_iv_length /= seq_count;

        uint search_mer_size = getDefaultSeedWeight( total_iv_length );
        if( search_mer_size < MIN_DNA_SEED_WEIGHT )
                return;         // The seed size is too small to be significant
        uint64 default_seed = getSeed( search_mer_size );
        
        //      Create sorted mer lists for the intervening gap region
        vector< boost::filesystem::path > delete_files;
        boolean create_succeeded = true;
        for( seqI = 0; seqI < seq_count; seqI++ ){
                gap_list.sml_table[ seqI ]->Clear();
                try{
                        if( debug_lcb_extension )
                                cerr << "Creating memory SML for seqI " << seqI << endl;
                        gap_list.sml_table[ seqI ]->Create( *(gap_list.seq_table[ seqI ]), default_seed );
                }catch(...){
                        create_succeeded = false;
                        break;
                }
        }
        if( !create_succeeded ){        
                // free memory consumed by any SMLs     
                for( seqI = 0; seqI < seq_count; seqI++ ){
                        gap_list.sml_table[ seqI ]->Clear();
                        delete gap_list.sml_table[ seqI ];
                }

                for( seqI = 0; seqI < seq_count; seqI++ ){
                        cerr << "Creating dmSML for seqI " << seqI << endl;
                        // presumably we ran out of memory and couldn't use a MemorySML.        
                        // try using a FileSML with external sort
                        string concat_file = CreateTempFileName("seqconcat");

                        concat_file += ".raw";  // need .raw extension to tell stupid libGenome it's a raw file
                        gnRAWSource::Write( *(gap_list.seq_table[ seqI ]), concat_file.c_str() );
                        delete_files.push_back( concat_file );
                        delete gap_list.seq_table[ seqI ];      // make sure memory gets freed!
                        cerr << "Wrote raw sequence for seqI " << seqI << endl;
                        gap_list.seq_table[ seqI ] = new gnSequence();
                        gap_list.seq_table[ seqI ]->LoadSource( concat_file.c_str() );
                        cerr << "Loaded sequence " << seqI << gap_list.seq_table[ seqI ]->length() << "b.p.\n";
                        string sml_file = CreateTempFileName("dmsml");  
                        DNAFileSML* sml = new DNAFileSML( sml_file.c_str() );   
                        gap_list.sml_table[ seqI ] = sml;       
                        sml->dmCreate( *(gap_list.seq_table[ seqI ]), default_seed );   
                        delete_files.push_back( sml_file );
                        delete_files.push_back( sml_file + ".coords" );
                }
        }

        //      Find all exact matches in the gap region
        nway_mh.Clear();
        nway_mh.FindMatches( gap_list );
        gap_list.MultiplicityFilter( seq_count );
//      nway_mh.GetMatchList( gap_list );

        // free memory used by SMLs!
        for( seqI = 0; seqI < seq_count; seqI++ ){
                gap_list.sml_table[ seqI ]->Clear();
                delete gap_list.sml_table[ seqI ];
        }
        
        if( debug_lcb_extension ){
                ofstream debug_extension_out( "new_extension_matches.txt" );
                WriteList( gap_list, debug_extension_out );
                debug_extension_out.close();
        }

        //      
        // If an N mask was used, transpose MUMs back into the previous 
        // sequence coordinates 
        //      
        if( !create_succeeded ){
                for( seqI = 0; seqI < seq_count; seqI++ )       
                        transposeMatches( gap_list, seqI, ((FileSML*)gap_list.sml_table[ seqI ])->getUsedCoordinates() );
        }
        //
        // Transpose MUMs back into their original sequence coordinates
        //
        for( seqI = 0; seqI < seq_count; seqI++ )
                transposeMatches( gap_list, seqI, iv_regions[ seqI ] );

        EliminateOverlaps( gap_list );
        gap_list.MultiplicityFilter( seq_count );
        // filter out matches that are too short
        gap_list.LengthFilter( MIN_ANCHOR_LENGTH );

        // free memory used by sequences!
        for( seqI = 0; seqI < seq_count; seqI++ )
                delete gap_list.seq_table[ seqI ];

        for( int delI = 0; delI < delete_files.size(); delI++ ) 
                boost::filesystem::remove( delete_files[delI] );

        new_matches.insert( new_matches.end(), gap_list.begin(), gap_list.end() );
}



class MatchLeftEndComparator {
public:
        MatchLeftEndComparator( unsigned seq = 0 ){
                m_seq = seq;
        }
        MatchLeftEndComparator( MatchLeftEndComparator& msc ){
                m_seq = msc.m_seq;
        }
        // TODO??  make this do a wraparound comparison if all is equal?
        boolean operator()(const AbstractMatch* a, const AbstractMatch* b) const{
                int32 start_diff = max( a->FirstStart(), m_seq ) - max( b->FirstStart(), m_seq );
                if(start_diff == 0){
                        uint32 m_count = a->SeqCount();
                        m_count = m_count <= b->SeqCount() ? m_count : b->SeqCount();
                        for(uint32 seqI = m_seq; seqI < m_count; seqI++){
                                int64 a_start = absolut( a->Start( seqI ) ), b_start = absolut( b->Start( seqI ) );
                                int64 diff = a_start - b_start;
                                if(a_start == (int64)NO_MATCH || b_start == (int64)NO_MATCH)
                                        continue;
                                else if(diff == 0)
                                        continue;
                                else
                                        return diff < 0;
                        }
                }
                return start_diff < 0;
        }
private:
        unsigned m_seq;
};

void transposeMatches( MatchList& mlist, uint seqI, const vector< int64 >& seq_regions ){
        if( seq_regions.size() < 2 )    
                return; // no work to be done here...

        uint matchI = 0;
        MatchLeftEndComparator msc( seqI );
        sort( mlist.begin(), mlist.end(), msc );
        uint regionI = 0;
        gnSeqI region_sum = seq_regions[ 1 ] - seq_regions[ 0 ];
        gnSeqI region_start_sum = 0;
        MatchList new_matches;

        for( ; matchI < mlist.size(); matchI++ ){
                // find the translated start coordinate for this match
                int64 trans_start = mlist[ matchI ]->Start( seqI );
                int64 iv_orig_start = trans_start;
                if( trans_start == 0 )
                        continue;
                while( region_sum < absolut( trans_start ) && regionI + 2 < seq_regions.size() ){
                        regionI += 2;
                        region_start_sum = region_sum;
                        region_sum += seq_regions[ regionI + 1 ] - seq_regions[ regionI ];
                }

                if( trans_start < 0 )
                        trans_start = -seq_regions[ regionI ] - ( -trans_start - region_start_sum ) + 1;
                else if( trans_start > 0 )
                        trans_start = seq_regions[ regionI ] + ( trans_start - region_start_sum ) - 1;

                int64 trans_end = mlist[ matchI ]->Start( seqI );
                trans_end += trans_end > 0 ? mlist[ matchI ]->Length() - 1: -mlist[ matchI ]->Length() + 1;
                
                mlist[ matchI ]->SetStart( seqI, trans_start );
                
                // this bad boy may need to be split
                gnSeqI end_region_sum = region_sum;
                gnSeqI end_prev_sum = region_start_sum;
                uint end_regionI = regionI;
                Match* cur_match = mlist[ matchI ];
                while( end_region_sum < absolut( trans_end ) && end_regionI + 2 < seq_regions.size() ){
                        end_regionI += 2;

                        Match* left_match = new Match( *cur_match );
                        // clip off the part going to the other match
                        if( left_match->Start( seqI ) < 0 ){
                                cur_match->CropStart( absolut( iv_orig_start ) + left_match->Length() - end_region_sum - 1);
                                left_match->CropEnd( cur_match->Length() );
                        }else{
                                cur_match->CropEnd( absolut( iv_orig_start ) + left_match->Length() - end_region_sum - 1);
                                left_match->CropStart( cur_match->Length() );
                        }

                        iv_orig_start += iv_orig_start > 0 ? cur_match->Length(): -cur_match->Length();

                        if( trans_start < 0 )
                                trans_start = -seq_regions[ end_regionI ] - ( -iv_orig_start - end_region_sum ) + 1;
                        else if( trans_start > 0 )
                                trans_start = seq_regions[ end_regionI ] + ( iv_orig_start - end_region_sum ) - 1;
                        
                        left_match->SetStart( seqI, trans_start );

                        cur_match = left_match;
                        new_matches.push_back( left_match );

                        end_prev_sum = end_region_sum;
                        end_region_sum += seq_regions[ end_regionI + 1 ] - seq_regions[ end_regionI ];

                }
//              if( end_region_sum == absolut( trans_end ) )
//                      cerr << "Beware of a possible bug in transposeMatches()\n";
        }
        
        // voila... coordinates are translated
        mlist.insert( mlist.end(), new_matches.begin(), new_matches.end() );
}

void ComputeLCBs( MatchList& meml, set<uint>& breakpoints, vector<MatchList>& lcb_list, vector<int64>& weights ){

        // there must be at least one end of a block defined
        if( breakpoints.size() < 1 )
                return;
                
        lcb_list.clear();
        weights.clear();
        
        // organize the LCBs into different MatchList instances

        set<uint>::iterator break_iter = breakpoints.begin();
        uint prev_break = 0;    // prev_break is the first match in the current block
        MatchList lcb = meml;
        for( ; break_iter != breakpoints.end(); break_iter++ ){
                lcb.clear();
                lcb.insert( lcb.begin(), meml.begin() + prev_break, meml.begin() + *break_iter + 1 );
                prev_break = *break_iter + 1;
                
                // code to filter LCBs based on their coverage
                uint64 coverage;
                GetLCBCoverage( lcb, coverage );
                weights.push_back( coverage );

                // add the new MatchList to the set if it made the cut
                lcb_list.push_back( lcb );
        }
}

void Aligner::Recursion( MatchList& r_list, Match* r_begin, Match* r_end, boolean nway_only ){
        try{
        gnSeqI gap_size = 0;
        uint seqI = 0;
//      gnSeqI min_gap_size = 0;
        boolean create_ok = true;
        // create gnSequences for each intervening region
        // create a MatchList for the intervening region
        MatchList gap_list;
        
        gap_list.seq_table.reserve( seq_count );
        gap_list.sml_table.reserve( seq_count );
        vector< int64 > starts;
        uint below_cutoff_count = 0;
// 
//      Get the sequence in the intervening gaps between these two matches
//
        for( seqI = 0; seqI < seq_count; seqI++ ){
                int64 gap_end = 0;
                int64 gap_start = 0;
                getInterveningCoordinates( r_list.seq_table, r_begin, r_end, seqI, gap_start, gap_end );
                if( (r_end && r_end->Start( seqI ) == NO_MATCH) ||
                        (r_begin && r_begin->Start( seqI ) == NO_MATCH )){
                        below_cutoff_count++;
                        cerr << "It's screwed up\n";
                        gap_list.seq_table.push_back( new gnSequence() );
                        gap_list.sml_table.push_back( new DNAMemorySML() );
                        continue;
                }
                if( gap_end < 0 && gap_start > 0 ){
                        create_ok = false;
                        cerr << "It's screwed up 2\n";
                        break; // bail out on directional inconsistency
                }else if( gap_end < 0 && gap_start > 0 ){
                        cerr << "It's screwed up 3\n";
                        create_ok = false;
                        break;  // bail out on directional inconsistency
                }
                int64 diff = gap_end - gap_start;
                diff = 0 < diff ? diff : 0;
                gap_size = diff < gap_size ? gap_size : diff;

                if( gap_start == 0 )
                        cerr << "scheiss\n";

                if( debug )
                        cout << r_list.seq_table[ seqI ]->length() << endl;

                if( diff < min_recursive_gap_length )
                        below_cutoff_count++;
                starts.push_back( gap_start );
                gnSequence* new_seq = new gnSequence( r_list.seq_table[ seqI ]->subseq( gap_start, diff ) );
                gap_list.seq_table.push_back( new_seq );
                gap_list.sml_table.push_back( new DNAMemorySML() );
        }
        
        // only perform recursive anchoring if the gapped regions are long enough
        // otherwise just let ClustalW do the work
        if( below_cutoff_count + 1 < seq_count ){
                if( nway_only )
                        nway_mh.Clear();
                else
                        gap_mh.get().Clear();

                multimap< uint, uint > mer_sizes;
                // calculate potential mer sizes for searches
                for( seqI = 0; seqI < seq_count; seqI++ ){
                        uint search_mer_size = getDefaultSeedWeight( gap_list.seq_table[ seqI ]->length() );
                        mer_sizes.insert( multimap< uint, uint >::value_type( search_mer_size, seqI ) );
                }
                multimap< uint, uint >::iterator mer_iter = mer_sizes.end();
                mer_iter--;
                vector< uint > search_seqs;
                while( mer_iter != mer_sizes.end() ){
                        uint prev_mer = mer_iter->first;
                        uint new_seqs = 0;
                        while( true ){
                                if( mer_iter->first < MIN_DNA_SEED_WEIGHT )
                                        break;
                                if( mer_iter->first == prev_mer || search_seqs.size() < 2 ){
                                        search_seqs.push_back( mer_iter->second );
                                        new_seqs++;
                                        if( mer_iter == mer_sizes.begin() ){
                                                mer_iter = mer_sizes.end();     // signify that the scan is complete
                                                break;
                                        }
                                        prev_mer = mer_iter->first;
                                        mer_iter--;
                                }else
                                        break;
                        }

                        if( search_seqs.size() < 2 )
                                break;
                        // look for MUMs
                        
                        //
                        //      Create sorted mer lists for the intervening gap region
                        //

                        uint64 default_seed = getSeed( prev_mer );
                        if( prev_mer < MIN_DNA_SEED_WEIGHT )
                                break;
                        for( uint seqI = 0; seqI < gap_list.seq_table.size(); seqI++ ){
                                gap_list.sml_table[ seqI ]->Clear();
                                gap_list.sml_table[ seqI ]->Create( *(gap_list.seq_table[ seqI ]), default_seed );
                        }
                        //
                        //      Find all exact matches in the gap region
                        //
                        MatchList cur_mems = gap_list;
                        cur_mems.clear();
                        if( nway_only ){
                                // no sense in searching for matches in subsets!!
                                if( search_seqs.size() < seq_count )
                                        continue;
                                nway_mh.ClearSequences();
                                nway_mh.FindMatches( cur_mems );
                        }else{
                                gap_mh.get().ClearSequences();
                                gap_mh.get().FindMatches( cur_mems );
                        }
                        for( size_t mI = 0; mI < cur_mems.size(); ++mI )
                                cur_mems[mI]->Free();
                        cur_mems.clear();
                }
                if( nway_only )
                        nway_mh.GetMatchList( gap_list );
                else
                        gap_mh.get().GetMatchList( gap_list );
                

                // delete overlaps/inclusions           
                EliminateOverlaps( gap_list );
                // mult. filter after EliminateOverlaps because e.o. may generate some subset matches
                if( nway_only )
                        gap_list.MultiplicityFilter( seq_count );
                
                // for anchor accuracy, throw out any anchors that are shorter than the minimum
                // anchor length after EliminateOverlaps()
                gap_list.LengthFilter( MIN_ANCHOR_LENGTH );

        //      if( min_gap_size < search_mer_size )
        //              create_ok = false;
                if( gap_list.size() > 0 && create_ok ){

        /*              if( debug ){
                                cout << "Starting mem: " << *r_begin << endl;
                                cout << "Next mem: " << *r_end << endl;
                                list<Match*>::iterator gappy_iter = gap_list.begin();
                                while( gappy_iter != gap_list.end() ){
                                        cout << **gappy_iter;
                                        cout << endl;
                                        gappy_iter++;
                                }
                        }
        */

                        // move all the matches that were found
                        vector< Match* >::iterator mum_iter = gap_list.begin();
                        for( ; mum_iter != gap_list.end(); ){
                                boolean add_ok = true;
                                for( uint seqI = 0; seqI < (*mum_iter)->SeqCount(); ++seqI ){
                                        int64 gap_start;
                                        if( (*mum_iter)->Start( seqI ) == NO_MATCH )
                                                continue;
                                        else if( (*mum_iter)->Start( seqI ) < 0 ){
                                                gap_start = r_begin != NULL ? -r_begin->End( seqI ) : 0;
                                                if( gap_start > 0 )
        //                                              gap_start = -r_end->Start( seqI ) + r_end->Length() - 1;
                                                        gap_start = r_end != NULL ? r_end->Start( seqI ) - r_end->Length() + 1 : 0;
                                                else if( r_begin )
                                                        add_ok = false;
                                                (*mum_iter)->SetStart( seqI, (*mum_iter)->Start( seqI ) + gap_start );
                                        }else{
                                                // insert them all before mem_iter
                                                gap_start = r_begin != NULL ? r_begin->End( seqI ) : 0;
                                                if( gap_start < 0 ){
                                                        gap_start = r_end != NULL ? r_end->Start( seqI ) - r_end->Length() + 1 : 0;
                                                        add_ok = false;
                                                }
                                                (*mum_iter)->SetStart( seqI, (*mum_iter)->Start( seqI ) + gap_start );
                                        }
                                }
                                if( add_ok )
                                        r_list.push_back( *mum_iter );
                                else{
                                        (*mum_iter)->Free();
                                        (*mum_iter) = NULL;
                                }
                                ++mum_iter;
                        }
        //              for( ; mum_iter != gap_list.end(); )
        //                      match_allocator.Free( *mum_iter );
                }
        }
        // delete sequences and smls
        for( uint seqI = 0; seqI < gap_list.seq_table.size(); ++seqI )
                delete gap_list.seq_table[ seqI ];
        for( uint seqI = 0; seqI < gap_list.sml_table.size(); ++seqI )
                delete gap_list.sml_table[ seqI ];
                
        gap_list.seq_table.clear();
        gap_list.sml_table.clear();
        
        }catch( gnException& gne ){
                cerr << gne << endl;
        }catch( exception& e ){
                cerr << e.what() << endl;
        }catch(...){
                cerr << "When I say 'ohhh' you say 'shit'!\n";
        }
}

// compute the gapped alignments between anchors in an LCB
void AlignLCBInParallel( bool collinear_genomes, mems::GappedAligner* gal, MatchList& mlist, Interval& iv, AlnProgressTracker& apt )
{
        // check whether this function can do anything useful...
        if( !collinear_genomes && mlist.size() < 2 ){
                iv.SetMatches( mlist );
                return;
        }
        size_t galI = 0;
        vector<GappedAlignment*> gapped_alns(mlist.size()+1, NULL);
        vector<int> success(gapped_alns.size(), 0);
        gnSeqI progress_base = apt.cur_leftend;
//#pragma omp parallel for
        for( int mI = 0; mI < mlist.size()-1; mI++ )
        {
                // align the region between mI and mI+1
                GappedAlignment ga(mlist.seq_table.size(),0);
                gapped_alns[mI] = ga.Copy();

                bool align_success = gal->Align( *(gapped_alns[mI]), mlist[mI], mlist[mI+1], mlist.seq_table );
                if(align_success)
                        success[mI] = 1;
                if(mI % 50 == 0 && mI > 0)
                {
                        // update and print progress
                        int done = 0;
                        for( int i = 0; i < gapped_alns.size(); i++ )
                                if(gapped_alns[i] != NULL)
                                        done++;
//#pragma omp critical
{
                        double cur_progress = ((double)(progress_base+done) / (double)apt.total_len)*100.0;
                        printProgress((uint)apt.prev_progress, (uint)cur_progress, cout);
                        apt.prev_progress = cur_progress;
}
                }
        }
        apt.cur_leftend += mlist.size()-1;

        // merge the alignments and anchors back together
        vector<AbstractMatch*> merged(mlist.size()*2 + 1);
        size_t mlistI = 0;
        size_t gappedI = 0;
        bool turn = true;
        size_t mJ = 0;

        // check if genomes are collinear and get the start and end alignments if necessary
        if(collinear_genomes)
        {
                GappedAlignment ga_tmp(mlist.seq_table.size(),0);
                GappedAlignment* ga = ga_tmp.Copy();
                bool align_success = gal->Align( *ga, NULL, mlist[0], mlist.seq_table );
                if(align_success)
                        merged[mJ++] = ga;
                gapped_alns[mlist.size()] = ga_tmp.Copy();
                align_success = gal->Align( *(gapped_alns[mlist.size()]), mlist.back(), NULL, mlist.seq_table );
                if(align_success)
                        success[mlist.size()] = 1;
        }
        for( ; mJ < merged.size() && mlistI < mlist.size();  )
        {
                if(turn)
                        merged[mJ++] = mlist[mlistI++];
                else if(success[gappedI])
                        merged[mJ++] = gapped_alns[gappedI++];
                else
                        gappedI++;
                turn = !turn;
        }
        // add the last alignment
        if( success[mlist.size()]==1 )
                merged[mJ++] = gapped_alns.back();
        merged.resize(mJ);

        iv.SetMatches(merged);
}

// compute the gapped alignments between anchors in an LCB
void Aligner::AlignLCB( MatchList& mlist, Interval& iv ){
        // check whether this function can do anything useful...
        if( !collinear_genomes && mlist.size() < 2 ){
                iv.SetMatches( mlist );
                return;
        }

        vector< AbstractMatch* > iv_matches;
        boolean debug_recurse = false;
        int64 config_value = 138500;
        int print_interval = 50;
        try{
        list< Match* > match_list;
        match_list.insert( match_list.end(), mlist.begin(), mlist.end() );
        mlist.clear();
        MatchList r_list = mlist;

        list< Match* >::iterator recurse_iter = match_list.begin();
        list< Match* >::iterator recurse_prev = match_list.begin();
        // scan ahead to the first n-way matches
        while( recurse_prev != match_list.end() && (*recurse_prev)->Multiplicity() != seq_count )
                ++recurse_prev;

        recurse_iter = recurse_prev;
        if( !collinear_genomes ){
                if( recurse_iter != match_list.end() )
                        ++recurse_iter;
                while( recurse_iter != match_list.end() && (*recurse_iter)->Multiplicity() != seq_count )
                        ++recurse_iter;
        }else
                cout << "Assuming collinear genomes...\n";
        
        uint memI = 0;
        uint matchI = 0;
        while( true ){
                if( memI >= print_interval && memI % print_interval == 0 || debug)
                        cout << "Number: " << memI << " match " << **recurse_prev << endl;
                ++memI;
                if( debug_recurse ){
                        cout << "Recursing on " << endl;
                        if( recurse_prev != match_list.end() )
                                cout << **recurse_prev << " and " << endl;
                        if( recurse_iter != match_list.end() )
                                cout << **recurse_iter << endl;
                }
                
                if( recurse_prev != match_list.end() && (*recurse_prev)->Start( 0 ) == config_value )
                        cout << "";
                
                // recurse on a pair of matches! 
                // this function should locate all matches between the two iterators
                // and add them to r_list               
                r_list.clear();
                GappedAlignment* cr = NULL;
                boolean align_success = false;
                
                Match* r_lend = NULL;
                Match* r_rend = NULL;
                if( recurse_iter != recurse_prev )
                        r_lend = *recurse_prev;
                if( recurse_iter != match_list.end() )
                        r_rend = *recurse_iter;

                // attempt a clustalW alignment
                cr = new GappedAlignment();
                align_success = gal->Align( *cr, r_lend, r_rend, r_list.seq_table );

                // add the gapped alignment to the Interval
                if( r_lend != NULL )
                        iv_matches.push_back( r_lend );
                if( align_success )
                        iv_matches.push_back( cr );

                // scan ahead to the next pair of n-way matches
                recurse_prev = recurse_iter;
                if( recurse_iter != match_list.end() )
                        ++recurse_iter;
                while( recurse_iter != match_list.end() && (*recurse_iter)->Multiplicity() != seq_count )
                        ++recurse_iter;

                if( ( recurse_iter == match_list.end() && !collinear_genomes ) ||
                                ( recurse_prev == match_list.end() && collinear_genomes ) )
                                break;
        }
        // get the last little bit at the end of the LCB.
        list< Match* >::iterator iter = recurse_prev;
        for( ; iter != recurse_iter; ++iter )
                iv_matches.push_back(*iter);

        mlist.insert( mlist.end(), match_list.begin(), match_list.end() );
        iv.SetMatches(iv_matches); 

        }catch( gnException& gne ){
                cerr << gne << endl;
        }catch(exception& e){
                cerr << e.what();
        }catch(...){
                cerr << "matrix exception?\n";
        }
}

// just search each intervening region once for matches, no gapped alignment...
void Aligner::SearchWithinLCB( MatchList& mlist, std::vector< search_cache_t >& new_cache, bool leftmost, bool rightmost){
        // check whether this function can do anything useful...
        if( !(leftmost || rightmost) && mlist.size() < 2 )
                return;

        boolean debug_recurse = false;
        int64 config_value = 138500;
        int print_interval = 50;

        try{
        list< Match* > match_list;
        match_list.insert( match_list.end(), mlist.begin(), mlist.end() );
        mlist.clear();
        MatchList r_list = mlist;

        list< Match* >::iterator recurse_iter = match_list.begin();
        list< Match* >::iterator recurse_prev = match_list.begin();
        if( !leftmost && recurse_iter != match_list.end() )
                ++recurse_iter;
        
        uint memI = 0;
        uint matchI = 0;
        while( recurse_prev != match_list.end() ){
                if( memI >= print_interval && memI % print_interval == 0 || debug)
                        cout << "Number: " << memI << " match " << **recurse_prev << endl;
                ++memI;
                if( debug_recurse ){
                        cout << "Recursing on " << endl;
                        if( recurse_prev != match_list.end() )
                                cout << **recurse_prev << " and " << endl;
                        if( recurse_iter != match_list.end() )
                                cout << **recurse_iter << endl;
                }
                
                
                // recurse on a pair of matches! 
                // this function should locate all matches between the two iterators
                // and add them to r_list               
                r_list.clear();
                Match* r_left = NULL;
                Match* r_right = NULL;
                if( recurse_iter == match_list.begin() && leftmost ){
                        r_left = NULL;
                        r_right = *recurse_iter;
                }else if( recurse_iter == match_list.end() && rightmost ){
                        r_left = *recurse_prev;
                        r_right = NULL;
                }else{
                        r_left = *recurse_prev;
                        r_right = *recurse_iter;
                }
                // check the cache to see whether this search has already been done!

                search_cache_t cacheval = make_pair( r_left, r_right );
                if( cacheval.first != NULL )
                        cacheval.first = cacheval.first->Copy();
                if( cacheval.second != NULL )
                        cacheval.second = cacheval.second->Copy();
                std::vector< search_cache_t >::iterator cache_entry = std::upper_bound( search_cache.begin(), search_cache.end(), cacheval, cache_comparator );
                if( cache_entry == search_cache.end() || 
                        (cache_comparator( cacheval, *cache_entry ) || cache_comparator( *cache_entry, cacheval )) )
                {
                        // search this region
                        Recursion( r_list, r_left, r_right, true );
                }
                new_cache.push_back( cacheval );

                if( debug_recurse ){
                        vector< Match* >::iterator r_iter = r_list.begin();
                        cout << "Found matches " << endl;
                        for(; r_iter != r_list.end(); ++r_iter )
                                cout << **r_iter << endl;
                }

                // insert any n-way matches into the match list
                for( matchI = 0; matchI < r_list.size(); ++matchI ){
                        if( r_list[ matchI ]->Multiplicity() == seq_count ){
                                match_list.insert( recurse_iter, r_list[ matchI ] );
                        }else
                        {
                                r_list[matchI]->Free();
                                r_list[matchI] = NULL;
                        }
                }

                // move ahead to the next pair of n-way matches
                recurse_prev = recurse_iter;
                if( recurse_iter != match_list.end() )
                        ++recurse_iter;
                
                // break early if we aren't assuming genome collinearity
                if( !rightmost && recurse_iter == match_list.end() )
                        break;
                        
        }

        mlist.insert( mlist.begin(), match_list.begin(), match_list.end() );

        }catch( gnException& gne ){
                cerr << gne << endl;
        }catch(exception& e){
                cerr << e.what();
        }catch(...){
                cerr << "matrix exception?\n";
        }

        // Multiplicity Filter...
        mlist.MultiplicityFilter( seq_count );
        EliminateOverlaps( mlist );
        // E.O. can create some matches of lower multiplicity
        mlist.MultiplicityFilter( seq_count );
}

void Aligner::consistencyCheck( uint lcb_count, vector< LCB >& adjacencies, vector< MatchList >& lcb_list, vector< int64 >& weights ){
        vector< LCB > tmp_adj = adjacencies;
        vector< MatchList > tmp_lcbs = lcb_list;
        vector< int64 > tmp_weights = weights;
        filterMatches( tmp_adj, tmp_lcbs, tmp_weights );
        MatchList emmlist;
        for( uint lcbI = 0; lcbI < tmp_lcbs.size(); lcbI++ )
                emmlist.insert( emmlist.end(), tmp_lcbs[ lcbI ].begin(), tmp_lcbs[ lcbI ].end() );
        set< uint > breakpoints;
        AaronsLCB( emmlist, breakpoints );
        
        // do the correct number of LCBs exist?
        if( lcb_count != tmp_lcbs.size() ){
                cerr << "lcb_count: " << lcb_count << "\ttmp_lcbs.size(): " << tmp_lcbs.size() << endl;
        }
        if( lcb_count != breakpoints.size() ){
                cerr << "lcb_count: " << lcb_count << "\tbreakpoints.size(): " << breakpoints.size() << endl;
        }
        if( tmp_lcbs.size() != breakpoints.size() ){
                cerr << "tmp_lcbs.size(): " << tmp_lcbs.size() << "\tbreakpoints.size(): " << breakpoints.size() << endl;
        }
}


int64 greedyBreakpointElimination( gnSeqI minimum_weight, vector< LCB >& adjacencies, vector< int64 >& weights, ostream* status_out ){
        // repeatedly remove the low weight LCBs until the minimum weight criteria is satisfied
        uint lcbI = 0;
        vector< uint > low_weight;
        bool have_weight = false;
        gnSeqI min_weight = 0;
        gnSeqI prev_min_weight = 0;
        uint min_lcb = 0;
        uint lcb_count = adjacencies.size();
        boolean debug_bp_elimination = false;
        uint current_lcbI = 0;  
        if( adjacencies.size() == 0 )
                return 0;       // nothing can be done
        uint seq_count = adjacencies[0].left_end.size();
        
        while( min_weight < minimum_weight ){
                if( lcb_count == 1 )
                        break;  // if only a single LCB remains, don't remove it

                while(true){
                        have_weight = false;
                        min_weight = 0;
                        current_lcbI = 0;       // always scan the entire set

                        // start with current_lcbI since everything up to it has already been scanned
                        for( lcbI = current_lcbI; lcbI < weights.size(); lcbI++ ){
                                if( adjacencies[ lcbI ].lcb_id != lcbI ){
                                        // this lcb has been removed or merged with another lcb
                                        continue;
                                }
                                if( weights[ lcbI ] < min_weight || !have_weight ){
                                        min_weight = weights[ lcbI ];
                                        min_lcb = lcbI;
                                        have_weight = true;
                                        if( min_weight == prev_min_weight && current_lcbI > 0 )
                                                break;  // we've already found a minimum
                                                                // weight LCB, stop here to save some searching
                                }
                        }
                        lcbI = min_lcb;
                        have_weight = false;
                        // if the min weight changed then scan the entire set from the beginning
                        if( prev_min_weight != min_weight ){
                                if( status_out != NULL )
                                        *status_out << "There are " << lcb_count << " LCBs with minimum weight " << min_weight << endl;

                                current_lcbI = 0;
                                prev_min_weight = min_weight;
                                continue;
                        }

                        // save time by skipping LCBs that have already been scanned
                        current_lcbI = min_lcb;
                        break;
                }
                
//              consistencyCheck( lcb_count, adjacencies, lcb_list, weights );
                if( min_weight >= minimum_weight )
                        break;

                // actually remove the LCBs now
                // (only remove a single LCB for now -- it's easier to calculate adjacencies)

                // remove this LCB
                adjacencies[ lcbI ].lcb_id = -2;
                
                // update adjacencies
                uint seqI;
                uint left_adj;
                uint right_adj;
                for( seqI = 0; seqI < seq_count; seqI++ ){
                        left_adj = adjacencies[ lcbI ].left_adjacency[ seqI ];
                        right_adj = adjacencies[ lcbI ].right_adjacency[ seqI ];
                        if( debug_bp_elimination ){
                                if( left_adj == -2 || right_adj == -2 ){
                                        cerr << "improper linking\n";
                                }
                                // for debugging, check for consistency:
                                if( left_adj != -1 && adjacencies[ left_adj ].right_adjacency[ seqI ] != lcbI )
                                        cerr << "Mutiny on the bounty!\n";
                                // for debugging, check for consistency
                                if( right_adj == adjacencies.size() )
                                        cerr << "Horrible Error -399a\n";
                                if( right_adj != -1 && adjacencies[ right_adj ].left_adjacency[ seqI ] != lcbI )
                                        cerr << "Mutiny on the bounty!\n";
                        }
                        if( left_adj != -1 )
                                adjacencies[ left_adj ].right_adjacency[ seqI ] = right_adj;
                        if( right_adj != -1 && right_adj != adjacencies.size() )
                                adjacencies[ right_adj ].left_adjacency[ seqI ] = left_adj;
                        
                }
                // just deleted an lcb, drop the lcb count
                lcb_count--;

                // check for collapse
                for( seqI = 0; seqI < seq_count; seqI++ ){
                        left_adj = adjacencies[ lcbI ].left_adjacency[ seqI ];
                        right_adj = adjacencies[ lcbI ].right_adjacency[ seqI ];
                        if( left_adj == -1 || right_adj == -1 )
                                continue;       // can't collapse with a non-existant LCB!

                        if( debug_bp_elimination ){
                                if( right_adj == adjacencies.size() )
                                        cerr << "Horrible Error -399a\n";
                                // check whether this LCB has already been merged
                                if( left_adj != adjacencies[ left_adj ].lcb_id ||
                                        right_adj != adjacencies[ right_adj ].lcb_id ){
                                        // because adjacency pointers are always updated to point to the 
                                        // representative entry of an LCB, the lcb_id and the array index
                                        // should always be identical
                                        cerr << "improper linking\n";
                                        continue;
                                }
                                if( left_adj == -2 || right_adj == -2 ){
                                        cerr << "improper linking\n";
                                }
                        }

                        // check whether the two LCBs are adjacent in each sequence
                        boolean orientation = adjacencies[ left_adj ].left_end[ seqI ] > 0 ? true : false;
                        uint seqJ;
                        for( seqJ = 0; seqJ < seq_count; seqJ++ ){
                                boolean j_orientation = adjacencies[ left_adj ].left_end[ seqJ ] > 0;
                                if( j_orientation == orientation &&
                                        adjacencies[ left_adj ].right_adjacency[ seqJ ] != right_adj )
                                        break;
                                if( j_orientation != orientation &&
                                        adjacencies[ left_adj ].left_adjacency[ seqJ ] != right_adj )
                                        break;
                                // check that they are both in the same orientation
                                if( adjacencies[ right_adj ].left_end[ seqJ ] > 0 != j_orientation )
                                        break;
                        }

                        if( seqJ != seq_count )
                                continue;
                        

                        // these two can be collapsed
                        // do it.  do it now.
                        adjacencies[ right_adj ].lcb_id = left_adj;
                        if( adjacencies[ right_adj ].lcb_id == -1 ||
                                adjacencies[ right_adj ].lcb_id == -2 )
                                cerr << "Trouble in the eleventh circle\n";
                        weights[ left_adj ] += weights[ right_adj ];
                        // unlink right_adj from the adjacency list and
                        // update left and right ends of left_adj
                        for( seqJ = 0; seqJ < seq_count; seqJ++ ){
                                boolean j_orientation = adjacencies[ left_adj ].left_end[ seqJ ] > 0;
                                uint rr_adj = adjacencies[ right_adj ].right_adjacency[ seqJ ];
                                uint rl_adj = adjacencies[ right_adj ].left_adjacency[ seqJ ];
                                if( j_orientation == orientation ){
                                        adjacencies[ left_adj ].right_end[ seqJ ] = adjacencies[ right_adj ].right_end[ seqJ ];
                                        adjacencies[ left_adj ].right_adjacency[ seqJ ] = rr_adj;
                                        if( rr_adj == adjacencies.size() )
                                                cerr << "Horrible Error -399a\n";
                                        if( rr_adj != -1 )
                                                adjacencies[ rr_adj ].left_adjacency[ seqJ ] = left_adj;
                                }else{
                                        adjacencies[ left_adj ].left_end[ seqJ ] = adjacencies[ right_adj ].left_end[ seqJ ];
                                        adjacencies[ left_adj ].left_adjacency[ seqJ ] = rl_adj;
                                        if( rl_adj == adjacencies.size() )
                                                cerr << "Horrible Error -399a\n";
                                        if( rl_adj != -1 )
                                                adjacencies[ rl_adj ].right_adjacency[ seqJ ] = left_adj;
                                }
                                // update lcbI's adjacency links to point nowhere
                                if( adjacencies[ lcbI ].left_adjacency[ seqJ ] == right_adj )
                                        adjacencies[ lcbI ].left_adjacency[ seqJ ] = left_adj;
                                if( adjacencies[ lcbI ].right_adjacency[ seqJ ] == right_adj )
                                        adjacencies[ lcbI ].right_adjacency[ seqJ ] = left_adj;


                        }
                        // just collapsed an lcb, decrement
                        lcb_count--;
                }
        }
        return min_weight;
}

class LCBLeftEndComp
{
public:
        LCBLeftEndComp() : ssc(0) {};
        bool operator()( const MatchList& a, const MatchList& b )
        {
                return ssc(a.front(), b.front());
        }
protected:
        SingleStartComparator<AbstractMatch> ssc;
};

void filterMatches( vector< LCB >& adjacencies, vector< MatchList >& lcb_list, vector< int64 >& weights ){
        if( lcb_list.size() < 1 )
                return;
        MatchList lcb_tmp = lcb_list[ 0 ];
        lcb_tmp.clear();
        vector< MatchList > filtered_lcbs = vector< MatchList >( lcb_list.size(), lcb_tmp );
        uint lcbI;
        for( lcbI = 0; lcbI < adjacencies.size(); lcbI++ ){
                if( adjacencies[ lcbI ].lcb_id == lcbI ){
                        filtered_lcbs[ lcbI ].insert( filtered_lcbs[ lcbI ].end(), lcb_list[ lcbI ].begin(), lcb_list[ lcbI ].end() );
                        continue;
                }
                if( adjacencies[ lcbI ].lcb_id == -1 ){
                        cerr << "weird";
                        continue;       // this one was removed
                }
                if( adjacencies[ lcbI ].lcb_id == -2 )
                        continue;       // this one was removed

                // this one points elsewhere
                // search and update the union/find structure for the target
                stack< uint > visited_lcbs;
                visited_lcbs.push( lcbI );
                uint cur_lcb = adjacencies[ lcbI ].lcb_id;
                while( adjacencies[ cur_lcb ].lcb_id != cur_lcb ){
                        visited_lcbs.push( cur_lcb );
                        cur_lcb = adjacencies[ cur_lcb ].lcb_id;
                        if( cur_lcb == -1 || cur_lcb == -2 ){
//                              cerr << "improper hoodidge\n";
                                break;  // this one points to an LCB that got deleted
                        }
                }
                while( visited_lcbs.size() > 0 ){
                        adjacencies[ visited_lcbs.top() ].lcb_id = cur_lcb;
                        visited_lcbs.pop();
                }
                // add this LCB's matches to the target LCB.
                if( cur_lcb != -1 && cur_lcb != -2 )
                        filtered_lcbs[ cur_lcb ].insert( filtered_lcbs[ cur_lcb ].end(), lcb_list[ lcbI ].begin(), lcb_list[ lcbI ].end() );
        }


        lcb_list.clear();
        vector< int64 > new_weights;
        for( lcbI = 0; lcbI < filtered_lcbs.size(); lcbI++ ){
                if( filtered_lcbs[ lcbI ].size() > 0 ){
                        lcb_list.push_back( filtered_lcbs[ lcbI ] );
                        uint64 wt = 0;
                        GetLCBCoverage( filtered_lcbs[lcbI], wt );
                        new_weights.push_back( wt );
//                      if( new_weights[ new_weights.size() - 1 ] != weights[ lcbI ] ){
//                              cerr << "Error: Have you lost weight Susan? difference: " << new_weights[ new_weights.size() - 1 ] - weights[ lcbI ] << "\n";
//                      }
                }
        }

        // sort the matches inside consolidated LCBs
        MatchStartComparator<AbstractMatch> msc( 0 );
        for( lcbI = 0; lcbI < lcb_list.size(); lcbI++ ){
                sort( lcb_list[ lcbI ].begin(), lcb_list[ lcbI ].end(), msc );
        }

        // sort the LCBs themselves
        LCBLeftEndComp llec;
        std::sort( lcb_list.begin(), lcb_list.end(), llec );

        // calculate the LCB adjacencies
        weights = new_weights;
        computeLCBAdjacencies_v2( lcb_list, weights, adjacencies );

}

void Aligner::WritePermutation( vector< LCB >& adjacencies, std::string out_filename )
{
        ofstream permutation_out( out_filename.c_str() );
        if( !permutation_out.is_open() )
        {
                cerr << "Error opening " << out_filename << endl;
                return;
        }
        for( int seqI = 0; seqI < seq_count; seqI++ )
        {
                // find the left-most LCB in this genome
                int left_lcb = 0;
                for( ; left_lcb < adjacencies.size(); left_lcb++ )
                {
                        uint left_adj = adjacencies[left_lcb].left_adjacency[seqI];
                        if( left_adj == -1 )
                                break;
                }
                // write out lcb id's in order
                for( uint lcbI = left_lcb; lcbI < adjacencies.size(); )
                {
                        if( lcbI != left_lcb )
                                permutation_out << '\t';
                        if( adjacencies[lcbI].left_end[seqI] < 0 )
                                permutation_out << "-";
                        permutation_out << adjacencies[lcbI].lcb_id;
                        lcbI = adjacencies[lcbI].right_adjacency[seqI];
                }
                permutation_out << endl;
        }
}

void WritePermutationCoordinates( IntervalList& perm_iv_list, std::string out_filename )
{
        ofstream perm_out( out_filename.c_str() );
        if( !perm_out.is_open() )
        {
                cerr << "Error opening \"" << out_filename << "\"\n";
                return;
        }
        perm_out << "#";
        for( size_t seqI = 0; seqI < perm_iv_list.seq_table.size(); ++seqI )
        {
                if( seqI > 0 )
                        perm_out << '\t';
                perm_out << "seq" << seqI << "_leftend\tseq" << seqI << "_rightend";
        }
        perm_out << endl;
        for( size_t ivI = 0; ivI < perm_iv_list.size(); ++ivI )
        {
                for( size_t seqI = 0; seqI < perm_iv_list.seq_table.size(); ++seqI )
                {
                        if( seqI > 0 )
                                perm_out << '\t';
                        if( perm_iv_list[ivI].Orientation(seqI) == AbstractMatch::reverse )
                                perm_out << '-';
                        perm_out << perm_iv_list[ivI].LeftEnd(seqI) << '\t';
                        if( perm_iv_list[ivI].Orientation(seqI) == AbstractMatch::reverse )
                                perm_out << '-';
                        perm_out << perm_iv_list[ivI].RightEnd(seqI);
                }
                perm_out << endl;
        }
}

void Aligner::RecursiveAnchorSearch( MatchList& mlist, gnSeqI minimum_weight, vector< MatchList >& LCB_list, boolean entire_genome, ostream* status_out ){

//
// Step 4) Identify regions of collinearity (LCBs) among the remaining n-way multi-MUMs
//
        uint lcbI;
        set<uint> breakpoints;
        vector< int64 > weights;
        vector< LCB > adjacencies;
        MatchList new_matches;
        new_matches.seq_table = mlist.seq_table;
        new_matches.seq_filename = mlist.seq_filename;

        if( mlist.size() == 0 )
                return;

        AaronsLCB( mlist, breakpoints );
        if( status_out )
                *status_out << "The " << mlist.size() << " matches constitute " << breakpoints.size() << " breakpoints\n";
        // organize the LCBs into different MatchList instances (inside of LCB_list)
        ComputeLCBs( mlist, breakpoints, LCB_list, weights );
        uint weightI;
        for( weightI = 0; weightI < weights.size(); weightI++ )
                if( weights[weightI] < cur_min_coverage || cur_min_coverage == -1 )
                        cur_min_coverage = weights[weightI];

        computeLCBAdjacencies_v2( LCB_list, weights, adjacencies );

        int cur_extension_round = 0;
        int64 total_weight = 0;
        int64 prev_total_weight = 0;
        weightI = 0;
        vector< vector< int64 > > prev_iv_regions;
        do {

//              for( ; weightI < weights.size(); weightI++ )
//                      total_weight += weights[ weightI ];

                int64 extension_weight = total_weight;
                int64 prev_extension_weight = total_weight;

                // only search outside existing LCBs on the whole-genome scale to save time
                if( entire_genome && extend_lcbs && total_weight != 0 &&
                        cur_extension_round < this->max_extension_iters )
                {
                        cur_extension_round++;
                        if( status_out )
                                *status_out << "Performing LCB extension\n";
                        vector< vector< int64 > > cur_iv_regions;
                        CreateGapSearchList( adjacencies, new_matches.seq_table, cur_iv_regions, entire_genome );
                        // only do the search if there's something new to search
                        if( prev_iv_regions != cur_iv_regions )
                        {
                                int local_round = 0;
                                do {
                                        local_round++;
                                        // search the gaps between the LCBs to extend the ends of LCBs
                                        new_matches.clear();
                                        vector< vector< int64 > > new_iv_regions;
                                        CreateGapSearchList( adjacencies, new_matches.seq_table, new_iv_regions, entire_genome );
                                        SearchLCBGaps( new_matches, new_iv_regions, nway_mh );
                                        mlist.insert( mlist.end(), new_matches.begin(), new_matches.end() );
                                        
                                        AaronsLCB( mlist, breakpoints );
                                        ComputeLCBs( mlist, breakpoints, LCB_list, weights );
                                        cur_min_coverage = *(std::min_element(weights.begin(), weights.end()));
                                        computeLCBAdjacencies_v2( LCB_list, weights, adjacencies );

                                        // calculate the new total LCB weight
                                        prev_extension_weight = extension_weight;
                                        extension_weight = 0;
                                        for( weightI = 0; weightI < weights.size(); weightI++ )
                                                extension_weight += weights[ weightI ];
                                        if( status_out )
                                                *status_out << "Previous weight: " << prev_extension_weight << " new weight: " << extension_weight << endl;
                                        if( prev_extension_weight > extension_weight ){
                                                cerr << "Error! Previous weight: " << prev_extension_weight << " new weight: " << extension_weight << endl;
                                        }
                                }while( extension_weight > prev_extension_weight && local_round < this->max_extension_iters);
                        }
                        swap( prev_iv_regions, cur_iv_regions );
                }
                
                // now search within LCBs
                if( currently_recursing && total_weight != 0 ){
                        vector< search_cache_t > new_cache;
                        for( lcbI = 0; lcbI < LCB_list.size(); lcbI++ ){
//                              if( status_out )
//                                      *status_out << "Searching in LCB: " << lcbI << endl;
                                int prev_size = LCB_list[ lcbI ].size();
                                bool leftmost = true;
                                for( int i = 0; leftmost && i < adjacencies[lcbI].left_adjacency.size(); i++ )
                                        if(adjacencies[lcbI].left_adjacency[i] != NO_ADJACENCY)
                                                leftmost = false;
                                bool rightmost = true;
                                for( int i = 0; rightmost && i < adjacencies[lcbI].right_adjacency.size(); i++ )
                                        if(adjacencies[lcbI].right_adjacency[i] != NO_ADJACENCY)
                                                rightmost = false;
                                SearchWithinLCB( LCB_list[ lcbI ], new_cache, leftmost, rightmost );
//                              if( status_out )
//                                      *status_out << "Gained " << LCB_list[ lcbI ].size() - prev_size << " matches\n";

                        }

                        // delete the previous search cache
                        swap( search_cache, new_cache );
                        for( size_t mI = 0; mI < new_cache.size(); mI++ )
                        {
                                if( new_cache[mI].first != NULL )
                                        new_cache[mI].first->Free();
                                if( new_cache[mI].second != NULL )
                                        new_cache[mI].second->Free();
                        }
                        new_cache.clear();
                        std::sort( search_cache.begin(), search_cache.end(), cache_comparator );
                }
                
                mlist.clear();
                for( lcbI = 0; lcbI < LCB_list.size(); lcbI++ ){
                        mlist.insert( mlist.end(), LCB_list[ lcbI ].begin(), LCB_list[ lcbI ].end() );
                }

                if( currently_recursing && total_weight != 0 ){
                        // remove low weight LCBs, while searching coalesced regions
                        AaronsLCB( mlist, breakpoints );
                        ComputeLCBs( mlist, breakpoints, LCB_list, weights );
                        computeLCBAdjacencies_v2( LCB_list, weights, adjacencies );
                        cur_min_coverage = *(std::min_element(weights.begin(), weights.end()));
                }

                
                // write  alist for debugging
//              ofstream debug_match_list( "debug_match_list.txt" );
//              mlist.WriteList( debug_match_list );
//              debug_match_list.close();

//
// Step 6) Use greedy breakpoint elimination to remove low-weight LCBs
//
                int64 cur_perm_weight = permutation_weight != -1 ? permutation_weight : minimum_weight;
                do{
                        vector<double> m_weights(weights.size());
                        std::copy( weights.begin(), weights.end(), m_weights.begin());
                        SimpleBreakpointScorer sbs(adjacencies, cur_perm_weight, this->collinear_genomes);
                        if( status_out )
                                (*status_out) << "Performing greedy breakpoint elimination (this may take some time)\n";

                        greedyBreakpointElimination_v4(adjacencies, m_weights, sbs, NULL, false);
//                      cur_min_coverage = greedyBreakpointElimination( cur_perm_weight, adjacencies, weights, status_out );
//                      MatchList deleted_matches;
                        filterMatches( adjacencies, LCB_list, weights );
                        cur_min_coverage = *(std::min_element(weights.begin(), weights.end()));
                        
                        mlist.clear();
                        for( lcbI = 0; lcbI < LCB_list.size(); lcbI++ ){
                                mlist.insert( mlist.end(), LCB_list[ lcbI ].begin(), LCB_list[ lcbI ].end() );
                        }
                        if( status_out )
                                *status_out << "Greedy breakpoint elimination leaves " << mlist.size() << " matches constituting " << LCB_list.size() << " LCBs covering at least " << cur_min_coverage << "b.p.\n";
                        
                        if( permutation_weight != -1 ){
                                // construct a filename
                                stringstream cur_perm_filename;
                                cur_perm_filename << permutation_filename << "." << cur_perm_weight / seq_count;
                                // output the permutation
                                WritePermutation( adjacencies, cur_perm_filename.str() );

                                // also write out condensed interval data for the permutation
                                cur_perm_filename << ".lcbs";
                                IntervalList perm_iv_list;
                                perm_iv_list.seq_filename = mlist.seq_filename;
                                perm_iv_list.seq_table = mlist.seq_table;
                                for( int permI = 0; permI < LCB_list.size(); permI++ ){
                                        vector< AbstractMatch* > perm_vector;
                                        perm_vector.push_back( LCB_list[permI].front() );
                                        if( LCB_list[permI].size() > 1 )
                                                perm_vector.push_back( LCB_list[permI].back() );
                                        Interval perm_iv(perm_vector.begin(), perm_vector.end());
                                        perm_iv_list.push_back(perm_iv);
                                }
                                WritePermutationCoordinates( perm_iv_list, cur_perm_filename.str() );

                                // get the current min weight
                                vector< int64 >::iterator min_w = std::min_element( weights.begin(), weights.end() );
                                // increment the current weight
                                cur_perm_weight = *min_w + seq_count;
                        }
                }while( cur_perm_weight < minimum_weight );
                // only enable recursive anchor search once we achieve
                // the desired weight threshold once -- for speed's sake
                if( recursive && entire_genome ){
                        currently_recursing = true;
                }

                // calculate the new total LCB weight
                prev_total_weight = total_weight;
                total_weight = 0;
                for( weightI = 0; weightI < weights.size(); weightI++ )
                        total_weight += weights[ weightI ];
                if( status_out )
                        *status_out << "Previous weight: " << prev_total_weight << " new weight: " << total_weight << endl;
        // the weight can shrink--this isn't an error condition
//              if( prev_total_weight > total_weight ){
//                      cerr << "Error! Previous weight: " << prev_total_weight << " new weight: " << total_weight << endl;
                        // write out the lcb lists
//              }

//
// Step 7) Repeat 4, 5 and 6 until the total weight stabilizes
//
        }while( total_weight != prev_total_weight );

        // delete the search cache
        for( size_t mI = 0; mI < search_cache.size(); mI++ )
        {
                if( search_cache[mI].first != NULL )
                        search_cache[mI].first->Free();
                if( search_cache[mI].second != NULL )
                        search_cache[mI].second->Free();
        }
}

void Aligner::align( MatchList& mlist, IntervalList& interval_list, double LCB_minimum_density, double LCB_minimum_range, boolean recursive, boolean extend_lcbs, boolean gapped_alignment, string tree_filename ){
        seq_count = mlist.seq_table.size();
        this->LCB_minimum_density = LCB_minimum_density;
        this->LCB_minimum_range = LCB_minimum_range;
        this->recursive = recursive;
        this->currently_recursing = false;
        this->extend_lcbs = extend_lcbs;
        this->gapped_alignment = gapped_alignment;

        // use LCB_minimum_range == -1 to indicate that all genomes are 
        // expected to be collinear
        this->collinear_genomes = LCB_minimum_range == -1;
        if( collinear_genomes )
                cout << "\nAssuming collinear genomes...\n";

        // set the nway_mh mask
        uint64 nway_mask = 1;
        nway_mask <<= seq_count;
        nway_mask--;
        nway_mh.SetMask( nway_mask );
                
        cout << "Starting with " << mlist.size() << " MUMs\n";
        
//
// Step 1) Eliminate overlaps among the multi-MUMs
//      
        // Remove linked inclusions
        EliminateOverlaps( mlist );
        cout << "Eliminating overlaps yields " << mlist.size() << " MUMs\n";

//
// Step 2) Compute a phylogenetic guide tree using the multi-MUMs
//

        bool guide_tree_loaded = false;
        MuscleInterface& mi = MuscleInterface::getMuscleInterface();    

        if( !guide_tree_loaded && (recursive || tree_filename != "") ){
                // Make a phylogenetic tree for ClustalW
                interval_list.seq_table = mlist.seq_table;
                interval_list.seq_filename = mlist.seq_filename;
                // use the identity matrix method and convert to a distance matrix
                NumericMatrix< double > distance;
                DistanceMatrix( mlist, distance );
                if( tree_filename == "" )
                        tree_filename = CreateTempFileName("guide_tree");
                mi.CreateTree( distance, tree_filename );
        }

//
// Step 3) Remove subset multi-MUMs
//
        // Multiplicity Filter...
        mlist.MultiplicityFilter( seq_count );
        cout << "Multiplicity filter gives " << mlist.size() << " MUMs\n";

        if( mlist.size() == 0 )
                return;
        
//
// Steps 4 through 7 are contained in RecursiveAnchorSearch
//
        vector< MatchList > LCB_list;
        RecursiveAnchorSearch( mlist, (gnSeqI)LCB_minimum_range, LCB_list, true, &cout );


//
// Step 8) Perform gapped alignment on each LCB using the anchors
//
        if( gapped_alignment && recursive )
                cout << "\nMaking final gapped alignment...\n";
        interval_list.clear();
        AlnProgressTracker apt;
        apt.cur_leftend = 0;
        apt.prev_progress = 0;
        apt.total_len = 0;
        for( uint lcbI = 0; lcbI < LCB_list.size(); lcbI++ )
                apt.total_len += LCB_list[lcbI].size()-1;
        for( uint lcbI = 0; lcbI < LCB_list.size(); lcbI++ ){
                Interval new_iv;
                interval_list.push_back( new_iv );
                Interval& iv = interval_list.back();
                if( !gapped_alignment || !recursive ){
                        iv.SetMatches( LCB_list[lcbI] );
                }else{
//                      AlignLCB( LCB_list[ lcbI ], iv );
                        AlignLCBInParallel( collinear_genomes || (LCB_list.size()==1), gal, LCB_list[ lcbI ], iv, apt );
                }
        }
        
        // finally add any unaligned regions to the interval list       
        if( gapped_alignment )
                addUnalignedIntervals( interval_list );
}

}       // namespace mems

---