src/report/runreport.cpp

// $Id: runreport.cpp,v 1.59 2011/03/07 06:08:51 bobgian Exp $

/*
  Copyright 2002  Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein

  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 <cassert>

#include "runreport.h"
#include "chainout.h"
#include "chainpack.h"
#include "force.h"
#include "forcesummary.h"
#include "stringx.h"
#include "types.h"
#include "registry.h"  // for GetDataPack() in MakeReport()

#ifdef DMALLOC_FUNC_CHECK
#include "/usr/local/include/dmalloc.h"
#endif

using namespace std;

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

string UNKNOWN_STRING = "<unknown>";

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

/****************************************************************
   This function predicts the end of the chain process, given
   the number of steps so far and the total required number of
   steps.  It is meant to be called by the chain manager.
*****************************************************************/
RunReport::RunReport(const ForceSummary& fs, verbosity_type progress)
    : m_level(progress),
      m_filelevel(registry.GetUserParameters().GetVerbosity()),
      m_forcesum(fs)
{
    m_profiletotal = static_cast<time_t>(0.0);
};

void RunReport::PrognoseRegion(const ChainPack& chpack, long int region,
                               long int steps, long int total)
{

    if (m_level == CONCISE || m_level == NONE) return;
    if (steps == total) return; //Let's not predict the current time.

    string timestr;

#ifdef NOTIME_FUNC
    // system does not provide a clock so no prognosis can be made
    timestr = "unknown";
#else
    time_t tnow = GetTime();
    time_t tstart = chpack.GetStartTime(region);   // start of the whole thing
    if (steps == 0) timestr = "unknown";
    else
    {
        double proportion = static_cast<double>(total) / static_cast<double>(steps);
        time_t predict = static_cast<time_t>(((tnow-tstart) * proportion) + tstart);
        timestr = PrintTime(predict,string("%c"));
    }
#endif

    string msg = "Predicted end of chains for this region:  ";
    msg += timestr + "\n";
    ReportNormal(msg);

} // PrognoseRegion

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

// This functions are used to predict the end of the profile
// process and print a report on progress so far.  When first
// beginning a profile, call RecordProfileStart() to set the time;
// then call PrognoseProfiles to report on progress.

void RunReport::RecordProfileStart()
{
    // Don't even try if there is no system clock.
#ifndef NOTIME_FUNC
    m_profilestart = GetTime();
#endif
    ReportUrgent("Beginning profiling, please be patient");
} // RecordProfileStart

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

void RunReport::PrognoseProfiles(long int thisprof, long int totprofs)
{
    if (m_level == NONE || m_level == CONCISE) return;
    long int profile = thisprof + 1;   // elsewhere they are counted starting from 0
    if (profile==totprofs)
    {
        m_profiletotal += (GetTime() - m_profilestart);
        string msg = "Finished profile "+ ToString(profile)
            + " of " + ToString(totprofs) + ".\n";
        ReportNormal(msg);
        return;
    }
    string timestr;

#ifdef NOTIME_FUNC
    timestr = "unknown";
#else
    time_t tnow = GetTime();
    double proportion = static_cast<double>(totprofs) / static_cast<double>(profile);
    time_t predict = static_cast<time_t>(((tnow - m_profilestart) * proportion) + m_profilestart);
    timestr = PrintTime(predict, string("%c"));
#endif
    string msg = "Finished profile " + ToString(profile) + " of " + ToString(totprofs)
        + ".  Predicted end of this set of profiles:  " + timestr + "\n";
    ReportNormal(msg);

} // PrognoseProfiles

void RunReport::PrognoseAll(const ChainPack& chpack, long int thisregion, long int totalregions)
{
    if (m_level == CONCISE || m_level == NONE) return;
    thisregion++;
    if ((thisregion == totalregions) && (totalregions == 1)) return;
    //The program is already done, since there's no 'overall' region to profile.

    string timestr;

#ifdef NOTIME_FUNC
    // system does not provide a clock so no prognosis can be made
    timestr = "unknown";
#else
    time_t tnow = GetTime();
    time_t tstart = chpack.GetStartTime();   // start of the whole thing
    double proportion = static_cast<double>(totalregions) / static_cast<double>(thisregion);
    time_t predict = static_cast<time_t>(((tnow-tstart) * proportion) + tstart);
    if (totalregions > 1)
    {
        //We will probably spend as much time profiling over all regions as
        // we spend profiling all the regions individually.
        predict += static_cast<time_t>(m_profiletotal * proportion);
    }
    timestr = PrintTime(predict,string("%c"));
#endif

    string msg = "Predicted end of this LAMARC run: ";
    msg += timestr;
    ReportNormal(msg);

} // PrognoseAll

/**************************************************************
   This function prints the "moving bar" which shows progress
   of the run.  It is meant to be called by the chain.
   You must call SetBarParameters before calling this.  Also,
   for PrintBar to work properly it must be called for the
   first time with steps = 0 and for the last time with steps = m_totalsteps.
**************************************************************/

void RunReport::PrintBar(long int steps)
{

#ifdef NOPROGRESSBAR
    if (steps == 0)
    {
        string msg = m_chaintype + " Chain " + ToString(m_chainno+1) + " of "
            + ToString(registry.GetChainParameters().GetNChains(m_chaintype_num)) + ":";
        ReportNormal(msg, false);
    }
    return;
#else
    if (m_level == CONCISE || m_level == NONE) return;

    long int i;
    static long int previous = 0;
    long int width = 12;

    if (steps == 0)
    { // Initial display of bar
        long int chaintot = registry.GetChainParameters().GetNChains(m_chaintype_num);
        string msg = m_chaintype + " Chain " + ToString(m_chainno+1) + " of "
            + ToString(chaintot) + ":";
        SaveOutput(msg, false);
        cout << endl << PrintTime(GetTime()) << "  ";
        cout << MakeJustified(m_chaintype + string(" chain "),-14)
             << Pretty(m_chainno + 1L, 3)
             << " of " << Pretty(chaintot, 3) << ":  ";
        cout << "[|                   ] ";
        cout << Pretty(steps + 1L, width);
        cout.flush();
        previous = 0;
        return;
    }

    long int counterDisplay = steps - m_burnin;

    if(counterDisplay % m_counter_display_increment == 0)
    {
        long int percent = static_cast<long int>(100.0 * static_cast<double>(steps)
                                                 / static_cast<double>(m_totalsteps));

        if (percent > previous && ((percent/5) * 5) - percent < EPSILON)
        {
            for (i = 0; i < width+22; i++) cout << "\b";   // back up
            for (i = 0; i < percent/5; i++)
            {
                cout << "=";
            }
            if (percent != 100) cout << "|";
            for (i = 0; i < 19-(percent/5); i++) cout << " ";
            cout << "] ";
        }
        else
        {
            for (i = 0; i < width; ++i) cout << "\b";
        }
        cout << Pretty(counterDisplay, width);
        if (steps == m_totalsteps) cout << " steps " << endl;
        cout.flush();
        previous = percent;
    }
#endif
} // PrintBar

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

void RunReport::AddForceToTable(const Force* force, StringVec1d& table,
                                const ChainOut& chout) const
{
    unsigned long int i;

    StringVec1d temptable = force->MakeChainParamReport(chout);

    // add the header
    table[0] += MakeCentered(force->GetShortparamname(), temptable[0].size());

    // Code added to deal with cases where the tables are not the same
    // size, as with divergence

    if (temptable.size() + 1 > table.size())
    {
        // put in blank entries
        assert(table.size() > 0);  // no entries at all??
        string blankentry(' ',table[1].size());
        for (i = table.size(); i < temptable.size() + 1; ++i)
        {
            table.push_back(blankentry);
        }
    }
    // add the parameter values
    for (i = 1; i < table.size() && i < temptable.size() + 1; ++i)
    {
        table[i] += "  " + temptable[i-1];
    }

} // AddForceToTable

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

void RunReport::MakeReport(const ChainOut& chout) const
{
    // Prepare header lines
    unsigned long int i;

    // Prepare tables of chain MLEs
    m_scalars.erase();
    m_tables.clear();

    const vector<Force *>& forces = m_forcesum.GetAllForces();

    // initialize the cross partition tables with naming info
    // the partition tables will be handled in the loop below
    const DataPack& datapack = registry.GetDataPack();
    StringVec1d xpartnames = datapack.GetAllCrossPartitionNames();

    StringVec1d tables2d;
    StringVec1d::iterator name;
    if (xpartnames.size() > 1)
    {
        // all cross-partitions will use the Class name "Class".
        tables2d.push_back(MakeJustified("Class",-22));
        for (name = xpartnames.begin(); name != xpartnames.end(); ++name)
            tables2d.push_back(MakeJustified(*name, -20));
    }

    StringVec2d tables3d;

    for (i = 0; i < forces.size(); ++i) // fill up all tables
    {
        if (forces[i]->ReportDimensionality() == 1) // scalar
        {
            m_scalars += MakeJustified(forces[i]->GetShortparamname(), -10) + "  ";
            m_scalars += forces[i]->MakeChainParamReport(chout)[0] + "  ";
        }
        else
        {                               // tabular
            // crosswise tables
            if (forces[i]->ReportDimensionality() == 2)
            {
                AddForceToTable(forces[i],tables2d,chout);
            }
            else
            {
                // partition tables
                assert (forces[i]->ReportDimensionality() == 3);
                StringVec1d table3d;
                StringVec1d partnames = datapack.
                    GetAllPartitionNames(forces[i]->GetTag());
                table3d.push_back(MakeJustified(forces[i]->GetClassName(),-23));
                for (name = partnames.begin(); name != partnames.end(); ++name)
                    table3d.push_back(MakeJustified(*name, -22));
                AddForceToTable(forces[i],table3d,chout);
                tables3d.push_back(table3d);
            }
        }
    }

    // put all the tables together into the big table
    m_tables.push_back(string("")); // add an empty line between m_scalars
    // and tables
    m_tables.assign(tables2d.begin(), tables2d.end());
    m_tables.push_back(string("")); // add an empty line between tables
    StringVec2d::iterator tab3d;
    for(tab3d = tables3d.begin(); tab3d != tables3d.end(); ++tab3d)
    {
        m_tables.insert(m_tables.end(), tab3d->begin(), tab3d->end());
        m_tables.push_back(string("")); // add an empty line between tables
    }

    // if heating, make swapping report
    long int numtemps = chout.GetNumtemps();
    if (numtemps > 1)
    {
        DoubleVec1d swaprates = chout.GetSwaprates();
        DoubleVec1d temperatures = chout.GetTemperatures();
        assert(swaprates.size() == static_cast<unsigned long int>(numtemps));
        assert(temperatures.size() == static_cast<unsigned long int>(numtemps));
        string tableline = "Temperature: ";
        long int temp;
        for (temp = 0; temp < numtemps; ++temp)
        {
            tableline += Pretty(temperatures[temp],7) + "  ";
        }
        m_tables.push_back(tableline);
        tableline = "  Swap rate: ";
        for (temp = 0; temp < numtemps-1; ++temp)
        {
            if (swaprates[temp] >= 0.0)
            {
                tableline += Pretty(100.0 * swaprates[temp], 7) + "  ";
            }
            else
            {
                tableline += "-------  ";
            }
        }
        tableline += "-------";
        m_tables.push_back(tableline);
    }
} // MakeReport

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

vector<string> RunReport::FormatReport(const ChainOut& chout, bool current,
                                       long int linelength) const
{
    bool bayesanalysis = registry.GetChainParameters().IsBayesian();
    string skiptimestamp;
    if (current) skiptimestamp = "  ";
    else skiptimestamp = "";

    vector<string> report;
    vector<string> temptable;

    // header
    //If this is a summary over regions or replicates, we do not want the
    // header information at all.  Everything is set 'FLAGLONG' or 'FLAGDOUBLE'
    // so we check that.  It's sort of a hack, but enh.
    if (!(chout.GetAccrate() == FLAGDOUBLE &&
          chout.GetLlikedata() == FLAGDOUBLE &&
          chout.GetNumBadTrees() == FLAGDOUBLE &&
          chout.GetTinyPopTrees() == FLAGDOUBLE &&
          chout.GetZeroDLTrees() == FLAGDOUBLE &&
          chout.GetStretchedTrees() == FLAGDOUBLE ))
    {
        string tempstring = "";
        if (current) tempstring = PrintTime(GetTime()) + "  ";
        tempstring += "Accepted ";
        if (chout.GetAccrate() == FLAGDOUBLE)
        {
            tempstring += UNKNOWN_STRING;
        }
        else
        {
            tempstring += MakeJustified(ToString(100.0 * chout.GetAccrate()), 5) + "%";
        }
        if (bayesanalysis)
            tempstring += " | Point Likelihood " + Pretty(chout.GetLlikemle(), 10);
        else
            tempstring += " | Posterior lnL " + Pretty(chout.GetLlikemle(), 10);
        tempstring += " | Data lnL ";
        if (chout.GetLlikedata() == NEG_MAX)
        {
            tempstring += UNKNOWN_STRING;
        }
        else
        {
            tempstring += Pretty(chout.GetLlikedata(), 10);
        }
        report.push_back(tempstring);
        long int badtrees = chout.GetNumBadTrees();
        long int tinytrees = chout.GetTinyPopTrees();
        long int zerodltrees = chout.GetZeroDLTrees();
        long int stretchedtrees = chout.GetStretchedTrees();
        if (badtrees == 0 && tinytrees == 0 &&
            (zerodltrees == 0 || zerodltrees == FLAGLONG) &&
            (stretchedtrees == 0 || stretchedtrees == FLAGLONG))
        {
            tempstring = "No trees discarded due to limit violations.";
            report.push_back(tempstring);
        }
        else
        {
            tempstring = "Trees discarded due to too many events: ";
            if (badtrees == FLAGLONG)
            {
                tempstring += UNKNOWN_STRING;
            }
            else
            {
                tempstring += Pretty(badtrees);
            }
            report.push_back(tempstring);

            tempstring = "Trees discarded due to too small population sizes: ";
            if (tinytrees == FLAGLONG)
            {
                tempstring += UNKNOWN_STRING;
            }
            else
            {
                tempstring += Pretty(tinytrees);
            }
            report.push_back(tempstring);

            tempstring = "Trees discarded due to an infinitesimal data likelihood: ";
            if (zerodltrees == FLAGLONG)
            {
                tempstring += UNKNOWN_STRING;
            }
            else
            {
                tempstring += Pretty(zerodltrees);
            }
            report.push_back(tempstring);

            tempstring = "Trees discarded due to extremely long branch lengths: ";
            if (stretchedtrees == FLAGLONG)
            {
                tempstring += UNKNOWN_STRING;
            }
            else
            {
                tempstring += Pretty(stretchedtrees);
            }
            report.push_back(tempstring);
        }
        // If multiple Arrangers, print an Arranger report
        ratemap arrates = chout.GetAllAccrates();
        ratemap::const_iterator rate;
        if (arrates.size() > 1)
        {
            for (rate = arrates.begin(); rate != arrates.end(); ++rate)
            {
                report.push_back(MakeJustified(rate->first + " accepted ", -28) + " " +
                                 MakeJustified(ToString(rate->second.first), 8) + "/" +
                                 ToString(rate->second.second) + " proposals");
            }
        }
    }

    LongVec1d bayesunique = chout.GetBayesUnique();
    if (bayesunique.size()>0)
    {
        //We have bayesian stuff.
        report.push_back("");
        report.push_back("Number of unique sampled values for each parameter:");
        const ParamVector paramvec(true);
        assert(bayesunique.size() == paramvec.size());
        bool toofew = false;
        for (unsigned long int pnum = 0; pnum < paramvec.size(); pnum++)
        {
            ParamStatus mystatus = paramvec[pnum].GetStatus();
            if (mystatus.Varies())
            {
                string numbayes = ToString(bayesunique[pnum]);
                if (bayesunique[pnum] < 50)
                {
                    toofew = true;
                    numbayes = numbayes + "*";
                }
                report.push_back(MakeJustified(numbayes, std::max(static_cast<int>(numbayes.size()), 6))
                                 + ": " + MakeJustified(paramvec[pnum].GetName(),-70));
            }
        }
        if (toofew)
        {
            report.push_back("* Note!  This parameter has too few sampled values to provide a reasonable estimate.  "
                             "Consider collecting more samples or narrowing your prior for this value.");
        }
        report.push_back("");
    }

    // Actually print the estimates
    if (m_scalars.size())
    {
        report.push_back(skiptimestamp + m_scalars);
    }
    temptable.insert(temptable.end(),m_tables.begin(),m_tables.end());

    temptable = Linewrap(temptable, linelength);

    transform(temptable.begin(),temptable.end(),temptable.begin(),
              bind1st(plus<string>(),skiptimestamp));
    report.insert(report.end(),temptable.begin(),temptable.end());

    return(report);
} // FormatReport

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

void RunReport::DisplayReport(const ChainOut& chout)
{
    MakeReport(chout);
    vector<string> rpt = FormatReport(chout,true,80);
    for (unsigned long int it = 0; it < rpt.size(); ++it)
    {
        ReportNormal(rpt[it], false, 80);
    }
} // DisplayReport

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

void RunReport::SetBarParameters(long int totsteps, long int burn, long int chno, long int chtype)
{
    m_totalsteps = totsteps;
    m_chainno = chno;
    m_burnin = burn;
    m_counter_display_increment = 10;
    if (m_totalsteps % 10L != 0)
    {
        m_counter_display_increment = registry.GetChainParameters().GetInterval(chtype);
    }
    m_burnpercent = 100.0 * static_cast<double>(burn)/static_cast<double>(totsteps);
    m_chaintype_num = chtype;
    switch (chtype)
    {
        case 0:
            m_chaintype = "Initial";
            break;
        case 1:
            m_chaintype = "Final";
            break;
        default:
            assert(false);
            m_chaintype = "Unknown";
            break;
    }
} // SetBarParameters

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

void RunReport::ReportUrgent(const string& msg, bool printtime, long int linelength)
{
    if (m_level != NONE)
    {
        PrettyPrint(msg, printtime, linelength);
    }
    if (m_filelevel != NONE)
    {
        SaveOutput(msg, printtime);
    }
}

void RunReport::ReportNormal(const string& msg, bool printtime, long int linelength)
{
    if (m_level == NORMAL || m_level==VERBOSE)
    {
        PrettyPrint(msg, printtime, linelength);
    }
    if (m_filelevel == NORMAL || m_filelevel==VERBOSE)
    {
        SaveOutput(msg, printtime);
    }
}

void RunReport::ReportNormal(const StringVec1d& msgs, bool printtime, long int linelength)
{
    if (msgs.size() > 0)
    {
        ReportNormal(msgs[0], printtime, linelength);
    }
    for (unsigned long int line = 1; line<msgs.size(); line++)
    {
        ReportNormal("          " + msgs[line], false, linelength);
    }
}

void RunReport::ReportChat(const string& msg, bool printtime, long int linelength)
{
    if (m_level==VERBOSE)
    {
        PrettyPrint(msg, printtime, linelength);
    }
    if (m_filelevel==VERBOSE)
    {
        SaveOutput(msg, printtime);
    }
}

void RunReport::ReportDebug(const string& msg, bool printtime, long int linelength)
{
#ifdef NDEBUG
    return;
#else
    string debugmsg = "Debug message: " + msg;
    PrettyPrint(debugmsg, printtime, linelength);
    SaveOutput(debugmsg, printtime);
#endif
}

void RunReport::ReportOnce(const string& msg, OnceKey key, bool doNormal)
{
    if (m_keyfrequencies.find(key) == m_keyfrequencies.end())
    {
        //This is the first time we've seen this message.  Add it, and report
        // it to the user, either with ReportNormal or ReportDebug.
        m_keyfrequencies.insert(make_pair(key, 1L));
        if (doNormal)
        {
            ReportNormal(msg);
            m_keyindices.insert(make_pair(static_cast<long int>(m_messages.size())-1, key));
        }
        else
        {
            ReportDebug(msg);
#ifndef NDEBUG
            m_keyindices.insert(make_pair(static_cast<long int>(m_messages.size())-1, key));
#endif
        }
        //Whether we actually reported it or not, we stick it in the index.  This
        // forces each OnceKey to come from exactly one call (or that all such
        // calls be the same version of doNormal); at present, there
        // are only four such calls, and they're all unique.
    }
    else
    {
        //We've seen this message before. Increment its count, but don't report.
        KeyToLongMapiter oldkey = m_keyfrequencies.find(key);
        long int count = oldkey->second;
        count++;
        m_keyfrequencies.erase(oldkey);
        m_keyfrequencies.insert(make_pair(key, count));
    }
}

//Private--use one of the above routines instead.
void RunReport::PrettyPrint(const string& msg, bool printtime, long int linelength)
{
    vector<string> msgline;
    msgline.push_back(msg);
    long int full_line = linelength - 10;

    string timestr, skiptime;
#ifdef NOTIME_FUNC
    // system does not provide a clock
    timestr = "";
    skiptime = "";
    full_line = linelength;
#else
    if (printtime)
    {
        time_t tnow = GetTime();
        timestr = PrintTime(tnow) + "  ";
        skiptime = "          ";
    }
    else
    {
        timestr = "";
        skiptime = "";
        full_line = linelength;
    }
#endif

    vector<string> msglines = Linewrap(msgline, full_line);

    vector<string>::iterator oneline = msglines.begin();
    cout << timestr << (*oneline) << endl;
    for(oneline++; oneline != msglines.end(); ++oneline)
        cout << skiptime << (*oneline) << endl;
}

void RunReport::SaveOutput(const string& msg, bool printtime)
{
    string timestr;
#ifdef NOTIME_FUNC
    // system does not provide a clock
    timestr = "";
#else
    if (printtime)
    {
        time_t tnow = GetTime();
        timestr = PrintTime(tnow) + "  ";
    }
    else
    {
        timestr = "";
    }
#endif
    m_messages.push_back(timestr + msg);
}

//____________________________________________________________________________________