src/tree/fc_status.cpp

// $Id: fc_status.cpp,v 1.14 2011/03/08 08:16:42 bobgian Exp $

/*
 * Copyright 2010 Mary Kuhner, Jon Yamato, Joseph Felsenstein, and Bob Giansiracusa
 *
 * This software is distributed free of charge for non-commercial use
 * and is copyrighted.  Of course, we do not guarantee that the software
 * works, and are not responsible for any damage you may cause or have.
 */

//------------------------------------------------------------------------------------

#include <set>
#include <map>

// Defines various symbols used to control debugging of experimental code blocks.
#include "local_build.h"

#include "fc_status.h"
#include "range.h"

//------------------------------------------------------------------------------------

using namespace std;

//------------------------------------------------------------------------------------
// Construct object representing live-site-branch-counts for arbitrary set of sites.
// Site indices are arbitrary; site indices (keys in map) are constructed as entries are made via update functions.
// Initial object is empty.

FC_Status::FC_Status()
    : m_site_index_lower_limit(MAXLONG), // LOWER end of lowest range inserted so far; the BEGINNING of the FC_Grid.
      m_site_index_upper_limit(0L)       // UPPER end of highest range inserted so far; the END of the FC_Grid.
{
    // Member m_fc_grid is an empty map, representing a live-site-branch-count of zero for all sites.
    // Member m_coalesced_sites contains an empty rangeset at object construction,
    // representing an empty set of final coalesced sites.
}

//------------------------------------------------------------------------------------
// Return set of ranges representing all sites that have coalesced, where "coalesced" means "have branch count one",
// whether the count was higher before and has decremented to one in "true coalescence" or the count just started off
// at zero and has now incremented to one.
//
// Objects (rangeset) returned are by-value copies and thus remain constant despite changes in internal state
// of FC_Status object (as horserace proceeds to other time intervals).  Client owns copied return-value objects.

rangeset FC_Status::Coalesced_Sites()
{
    // This function is called often and simply returns a value pre-computed and cached by Adjust_FC_Counts().
    return m_coalesced_sites;
}

//------------------------------------------------------------------------------------
// Verification of final-coalescence algorithm: test for equality of two FC_Status objects computed different ways.
// Simply checks that both internal subobjects are equal (equality defined by the respective container classes).

bool operator==(const FC_Status & fc_counter_lhs, const FC_Status & fc_counter_rhs)
{
    return ((fc_counter_lhs.m_fc_grid == fc_counter_rhs.m_fc_grid)
            &&
            (fc_counter_lhs.m_coalesced_sites == fc_counter_rhs.m_coalesced_sites));
}

//----------------------------------------
// Both operators are provided for completeness.

inline bool operator!=(const FC_Status & fc_counter_lhs, const FC_Status & fc_counter_rhs)
{
    return !(fc_counter_lhs == fc_counter_rhs);
}

//------------------------------------------------------------------------------------
// Increment_FC_Counts() and Decrement_FC_Counts() are defined inline in class declaration.
// They both call this function, which does all the work.

void FC_Status::Adjust_FC_Counts(const rangeset & sites_to_adjust, long int adjustment)
{
    // Even though "adjustment" can only be +1 or -1, all the operands it is added to are long ints,
    // so it will be coerced to long anyway.  Might as well start out that way.

    // Used mostly for debugging and error-checking.  Once fully debugged, may remove this.
    int iteration(0);

    // "Current range" is the range that iterator "range_iter" is pointing at for current iteration.
    for(rangeset::const_iterator range_iter(sites_to_adjust.begin()) ;
        range_iter != sites_to_adjust.end() ;
        ++range_iter)
    {
        // This loop inserts new nodes if there is none pre-existing at the site index marking the LOWER or UPPER
        // edge of a range (one of the ranges in the "sites_to_adjust" set).  A new node always indicates a change
        // in the associated branch count.  However, it is possible that one endpoint (perhaps both) of an input
        // range might coincide with pre-existing node in the FC_Grid (if the endpoint index happens to be a value
        // at which the count stored internally changes).  In such a "coincident edge endpoint" case, the
        // pre-existing node can simply be updated to reflect the new count for the region it represents (upward).

        const long int range_lower_index(range_iter->first);
        const long int range_upper_index(range_iter->second);

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 10)
        cerr << "Adjust_FC_Counts[1]: Start new iteration.  Iteration: " << ++iteration << endl
             << "Current range:            ( " << range_lower_index << ' ' << range_upper_index << " )" << endl;
        PrintFCStatus(false, sites_to_adjust, adjustment);
        cerr << endl;
#endif

        if(range_lower_index < range_upper_index) // Test for empty current range.  Shouldn't happen  ...
        {
            // Iterator "lower_edge" starts pointing to first node with key equal to or greater than current range
            // LOWER index or equal to m_fc_grid.end() if no nodes have yet been inserted into the FC_Grid above the
            // current range (because this is the first being inserted or just the highest so far being inserted).
            FC_Grid::iterator lower_edge(m_fc_grid.lower_bound(range_lower_index));

            // No nodes have yet been inserted into FC_Grid with keys HIGHER than (or EQUAL to) the LOWER end
            // of the current range (this range is either the first one inserted or the highest one inserted so far).
            // Insert one with branch count zero (will be adjusted later).
            if(lower_edge == m_fc_grid.end())
            {
                // "Hint" form of insert(), using END iterator.  Now "lower_edge" points to the node just inserted.
                lower_edge = m_fc_grid.insert(lower_edge, pair<const long int, long int>(range_lower_index, 0L));

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 7)
                cerr << "Adjust_FC_Counts[2]: Inserted new node, LOWER edge.  Iteration: "
                     << iteration << endl
                     << "Current range:            ( "
                     << range_lower_index << ' ' << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif

                if(range_lower_index < m_site_index_lower_limit)
                {
                    m_site_index_lower_limit = range_lower_index;

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 8)
                    cerr << "Adjust_FC_Counts[3]: Updated FC_Grid's LOWER range edge limit.  Iteration: "
                         << iteration << endl
                         << "Current range:            ( "
                         << range_lower_index << ' ' << range_upper_index << " )" << endl;
                    PrintFCStatus(false, sites_to_adjust, adjustment);
                    cerr << endl;
#endif
                }
            }

            // Find node marking LOWER end of current range.  If it exists, denote it with iterator "lower_edge"
            // and leave its pre-adjustment count in place (to be updated in traversal of current range later).
            else if(lower_edge->first == range_lower_index)
            {
                // Iterator "lower_edge" points to pre-existing node marking current range LOWER end.
                // This marks the BEGINNING of the upcoming traversal (the first node that WILL be adjusted).

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 5)
                cerr << "Adjust_FC_Counts[4]: Will adjust node, LOWER edge.  Iteration: "
                     << iteration << endl
                     << "Current range:            ( "
                     << range_lower_index << ' ' << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif
            }

            // If no existing node was found with key searched for above (ie, LOWER end of current range is between
            // two existing nodes in FC_Grid), find first node going lower from the one pointed at by iterator
            // "lower_edge)" (which is the first existing node with key going UP from current range LOWER end),
            // extract its count, and insert new node marking actual LOWER end of current range (using just-extracted
            // pre-adjustment branch count).  If NO node exists with a key below the value being searched for, DON'T
            // decrement the iterator and DO use zero as the branch count.
            else
            {
                long int branch_count = 0L;     // Default branch count if no node exists "below" insert location.

                // Don't attempt to point iterator "below" the LOWEST node in the entire FC_Grid.
                if(range_lower_index > m_site_index_lower_limit)
                {
                    // If a node "below" LOWER end of current range exists, read branch count and reset iterator.
                    branch_count = (--lower_edge)->second;
                    ++lower_edge;               // Reset iterator to value outside this conditional.
                }

                // "Hint" form of insert() used here, with "hint" being original position returned by lower_bound().
                // Now "lower_edge" points to the node just inserted.
                lower_edge = m_fc_grid.insert(lower_edge, pair<const long, long>(range_lower_index, branch_count));

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 6)
                cerr << "Adjust_FC_Counts[5]: Inserted new node, LOWER edge.  Iteration: " << iteration << endl
                     << "Current range:            ( " << range_lower_index << ' '
                     << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif

                if(range_lower_index < m_site_index_lower_limit)
                {
                    m_site_index_lower_limit = range_lower_index;

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 8)
                    cerr << "Adjust_FC_Counts[6]: Updated FC_Grid's LOWER range edge limit.  Iteration: "
                         << iteration << endl
                         << "Current range:            ( "
                         << range_lower_index << ' ' << range_upper_index << " )" << endl;
                    PrintFCStatus(false, sites_to_adjust, adjustment);
                    cerr << endl;
#endif
                }
            }
            // Now iterator "lower_edge" points to node (not yet adjusted) at included LOWER end of current range.

            // Iterator "upper_edge" starts pointing to first node with key equal to or greater than current range
            // UPPER index or equal to m_fc_grid.end() if no nodes have yet been inserted into the FC_Grid above the
            // current range (because this is first range inserted or just highest range so far inserted).
            FC_Grid::iterator upper_edge(m_fc_grid.lower_bound(range_upper_index));

            // No nodes have been inserted into FC_Grid with keys HIGHER than (or EQUAL to) the UPPER end of current
            // range (this range is either first inserted or highest inserted so far).  Insert one with branch count
            // zero (to be adjusted later).
            if(upper_edge == m_fc_grid.end())
            {
                // "Hint" form of insert(), with "hint" being END iterator.  Now "upper_edge" points to node inserted.
                upper_edge = m_fc_grid.insert(upper_edge, pair<const long int, long int>(range_upper_index, 0L));

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 7)
                cerr << "Adjust_FC_Counts[7]: Inserted new node, UPPER edge.  Iteration: " << iteration << endl
                     << "Current range:            ( " << range_lower_index
                     << ' ' << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif

                if(range_upper_index > m_site_index_upper_limit)
                {
                    m_site_index_upper_limit = range_upper_index;

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 8)
                    cerr << "Adjust_FC_Counts[8]: Updated FC_Grid's UPPER range edge limit.  Iteration: "
                         << iteration << endl << "Current range:            ( "
                         << range_lower_index << ' ' << range_upper_index << " )" << endl;
                    PrintFCStatus(false, sites_to_adjust, adjustment);
                    cerr << endl;
#endif
                }
            }

            // Find node marking UPPER end of current range.  If it exists, denote it with iterator "upper_edge"
            // and leave its pre-adjustment count in place (to be updated in traversal of current range later).
            else if(upper_edge->first == range_upper_index)
            {
                // Iterator "upper_edge" points to pre-existing node marking current range UPPER end.
                // This marks the END of the upcoming traversal (the first node that will NOT be adjusted).

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 5)
                cerr << "Adjust_FC_Counts[9]: Will adjust node, UPPER edge.  Iteration: "
                     << iteration << endl << "Current range:            ( "
                     << range_lower_index << ' ' << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif
            }

            // If no existing node was found with key searched for above (ie, UPPER end of current range is between
            // two existing nodes in FC_Grid), find first node going lower from the one pointed at by iterator
            // "upper_edge)" (which is the first existing node with key going DOWN from current range UPPER end),
            // extract its count, and insert new node marking actual UPPER end of current range (using just-extracted
            // pre-adjustment branch count).
            //
            // If NO node exists with a key below the value being searched for,
            // DON'T decrement the iterator and DO use zero as the branch count.
            else
            {
                long int branch_count = 0L;     // Branch count if no node exists below current insert location.

                // Don't attempt to point iterator "below" the LOWEST node in the entire FC_Grid.
                if(range_upper_index > m_site_index_lower_limit)
                {
                    // If a node "below" UPPER end of current range exists, read branch count and reset iterator.
                    branch_count = (--upper_edge)->second;
                    ++upper_edge;               // Reset iterator to value outside this conditional.
                }

                // "Hint" form of insert(), with "hint" being original position returned by lower_bound() above.
                // Now "upper_edge" points to the node just inserted.
                upper_edge = m_fc_grid.insert(upper_edge, pair<const long, long>(range_upper_index, branch_count));

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 6)
                cerr << "Adjust_FC_Counts[10]: Inserted new node, UPPER edge.  Iteration: "
                     << iteration << endl << "Current range:            ( " << range_lower_index << ' '
                     << range_upper_index << " )" << endl;
                PrintFCStatus(false, sites_to_adjust, adjustment);
                cerr << endl;
#endif

                if(range_upper_index > m_site_index_upper_limit)
                {
                    m_site_index_upper_limit = range_upper_index;

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 8)
                    cerr << "Adjust_FC_Counts[11]: Updated FC_Grid's UPPER range edge limit.  Iteration: "
                         << iteration << endl << "Current range:            ( "
                         << range_lower_index << ' ' << range_upper_index << " )" << endl;
                    PrintFCStatus(false, sites_to_adjust, adjustment);
                    cerr << endl;
#endif
                }
            }
            // Now iterator "upper_edge" points to node (not yet adjusted) at included UPPER end of current range
            // (whether pre-existing or just inserted).
            //
            // Now we can traverse all nodes within the current range, starting with the node marking the LOWER edge
            // (inclusive - will be adjusted) and continuing to (but excluded - will not be adjusted) the node
            // marking UPPER edge of current range.
            for(FC_Grid::iterator node_iter(lower_edge); node_iter != upper_edge; ++node_iter)
            {
                node_iter->second += adjustment;
#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 10)
                cerr << "Adjust_FC_Counts[12]: During traversal, per node, after adjustment: Node key: "
                     << node_iter->first << ", Node value: " << node_iter->second << endl << endl;
#endif
            }

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 9)
            cerr << "Adjust_FC_Counts[13]: After traversal and node adjustments.  Iteration: " << iteration << endl
                 << "Current range:            ( " << range_lower_index << ' ' << range_upper_index << " )" << endl;
            PrintFCStatus(false, sites_to_adjust, adjustment);
            cerr << endl;
#endif
        }
        else
        {
            // ERROR: Empty or malformed range presented.  Here now for debugging purposes.
            cerr << "Adjust_FC_Counts[14]: Empty or malformed range presented.  Iteration: " << iteration << endl;
#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 1)
            PrintFCStatus(true, sites_to_adjust, adjustment);
#endif
            cerr << endl << "Emergency exit." << endl;
            exit(1);
        }
    }

    // Now re-compute and cache the "coalesced sites", which is simple a rangeset of all sites whose branch count
    // is now one. Note that a site may have a count of one either via coalescence (decrement of count to one from a
    // higher value) or simply because the count is now one (initialized to zero by constructor, incremented
    // once, and never changed again).  Thus we must recompute the cache after a range's incremented branch count
    // (first increment, in case no more occur) as well as after a decremented branch count (a "true" coalescence).
    //
    // Hopefully soon this section will contain a beautiful functional "STL algorithm-style" method for updating a
    // data-structure.  For now, we simply start from scratch each time and build a set holding the appropriate data.
    // For now, since we are scanning the entire FC_Grid here in order to construct the cached rangeset each time,
    // we also do a bunch of consistency tests which can be eliminated once this code passes its startup sanity tests.

    // Clear it and start from scratch each time.  Later we will update rather than rebuild this object.
    m_coalesced_sites.clear();

    FC_Grid::iterator fc_grid_limit(m_fc_grid.end());
    long int previous_site_index = 0L;
    long int previous_branch_count = 0L;
    long int current_site_index = 0L;
    long int current_branch_count = 0L;

    iteration = 0;                              // Reset iteration variable to count iters through a different loop.

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 9)
    cerr << "Adjust_FC_Counts[15]: Beginning traversal on computing cached coalescences.  Previous: ("
         << previous_site_index << ' ' << previous_branch_count << "),  Current: ("
         << current_site_index << ' ' << current_branch_count << ')' << endl;
    PrintFCStatus(true, sites_to_adjust, adjustment);
    cerr << endl;
#endif

    for(FC_Grid::iterator node_iter( m_fc_grid.begin()) ; node_iter != fc_grid_limit ; /* increment inside loop */ )
    {
        current_site_index = node_iter->first;
        current_branch_count = node_iter->second;

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 10)
        cerr << "Adjust_FC_Counts[16]: Traversal iteration on computing cached coalescences.  Iteration: "
             << ++iteration << "  Previous: (" << previous_site_index << ' ' << previous_branch_count
             << "),  Current: (" << current_site_index << ' ' << current_branch_count << ')' << endl;
        PrintFCStatus(true, sites_to_adjust, adjustment);
        cerr << endl;
#endif

        if(current_branch_count < 0L)
        {
            cerr << "Adjust_FC_Counts[17]: Negative branch count on computing cached coalescences.  Iteration: "
                 << iteration << "  Previous: (" << previous_site_index << ' ' << previous_branch_count
                 << "),  Current: (" << current_site_index << ' ' << current_branch_count << ')' << endl;
#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 1)
            PrintFCStatus(true, sites_to_adjust, adjustment);
#endif
            cerr << endl << "Emergency exit." << endl;
            exit(1);
        }

        if(previous_branch_count == current_branch_count)
        {
            // If two adjacent ranges meet in a node and have the same (post-adjustment) branch counts, we can merge
            // those regions by deleting the node marking their junction.  It is also OK to delete the LOWEST and/or
            // HIGHEST node(s) if their branch count values are zero (on the LOWER end)
            // or if they match that to the left (on the UPPER end).

            FC_Grid::iterator delete_me(node_iter);
            ++node_iter;                // Increment iterator before deleting node; then continue traversal.
            m_fc_grid.erase(delete_me);

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 3)
            cerr << "Adjust_FC_Counts[18]: After region merge, erased node at LOWER edge.  Iteration: "
                 << iteration << endl << "Inferior range:           (" << previous_site_index << ' '
                 << current_site_index << ")  Branch count: " << current_branch_count << endl;
            PrintFCStatus(true, sites_to_adjust, adjustment);
            cerr << endl;
#endif
        }
        else
        {
            // Regions abutting at this node have different branch counts.  If LOWER region (to left of current node)
            // has count of 1, then it represents a final coalescence (or a count climbing up from zero).  Insert that
            // region into the cache.  Otherwise, keep looking.
            if(previous_branch_count == 1L)
            {
                // "Hinted" insertion location is always just preceding the END (UPPER end of entire FC_Grid) iterator
                // position, because we are always doing a "push_back" insertion.
                m_coalesced_sites.insert(m_coalesced_sites.end(),
                                         pair<const long, long>(previous_site_index, current_site_index));

#if (LAMARC_QA_FCTEST_DEBUGLEVEL >= 2)
                cerr << "Adjust_FC_Counts[19]: Inserted range pair " << previous_site_index << ':'
                     << (current_site_index - 1) << " on computing cached coalescences.  Iteration: "
                     << iteration << "  Previous: (" << previous_site_index << ' ' << previous_branch_count
                     << "),  Current: (" << current_site_index << ' ' << current_branch_count << ')' << endl;
                PrintFCStatus(true, sites_to_adjust, adjustment);
                cerr << endl;
#endif
            }

            ++node_iter; // Increment inside loop since we already incremented inside other branch of conditional.

            previous_site_index = current_site_index;
            previous_branch_count = current_branch_count;
        }
    }
}

//------------------------------------------------------------------------------------
// JDEBUG: quick and dirty GDB Debugging functions.

long int FC_Status::GridSize() const
{
    return m_fc_grid.size();
}

void FC_Status::PrintGrid() const
{
    for(FC_Grid::const_iterator node_iter = m_fc_grid.begin(); node_iter != m_fc_grid.end(); ++node_iter)
    {
        cerr << node_iter->first << ":" << node_iter->second << ",";
    }

    cout << endl;
}

//------------------------------------------------------------------------------------
// Testing and debugging; remove in production version.

#if (LAMARC_QA_FCTEST_DEBUGLEVEL > 0)

void FC_Status::PrintFCStatus(bool print_output, const rangeset & sites_to_adjust, long int adjustment)
{
    cerr << "Input rangeset:           " << ToString(sites_to_adjust) << "  <==>  ( ";

    for(rangeset::const_iterator range_iter(sites_to_adjust.begin()) ;
        range_iter != sites_to_adjust.end() ;
        ++range_iter)
    {
        cerr << '(' << range_iter->first << ' ' << range_iter->second << ") ";
    }

    cerr << ')' << endl
         << "Adjustment:               " << adjustment << endl;

    cerr << "site_index_limits:        " << m_site_index_lower_limit << " (lower), "
         << m_site_index_upper_limit << " (upper)" << endl;

    cerr << "m_fc_grid:                ( ";
    for(FC_Grid::const_iterator node_iter(m_fc_grid.begin()); node_iter != m_fc_grid.end(); ++node_iter)
    {
        cerr << '(' << node_iter->first << ' ' << node_iter->second << ") ";
    }
    cerr << ')' << endl;

    if(print_output)
    {
        cerr << "m_coalesced_sites:        ( ";
        for(rangeset::const_iterator range_iter(m_coalesced_sites.begin()) ;
            range_iter != m_coalesced_sites.end() ;
            ++range_iter)
        {
            cerr << '(' << range_iter->first << ' ' << range_iter->second << ") ";
        }
        cerr << ')' << endl;
    }
}

#endif

//____________________________________________________________________________________