xref: /dports/biology/lamarc/lamarc-2.1.8/src/force/forcesummary.cpp

// $Id: forcesummary.cpp,v 1.61 2011/04/23 02:02:49 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 <cmath>

#include "local_build.h"

#include "chainout.h"
#include "datapack.h"                   // for access to Region functions in SummarizeData()
#include "force.h"                      // to create Force objects
#include "forceparam.h"
#include "forcesummary.h"
#include "plotstat.h"
#include "region.h"
#include "runreport.h"
#include "summary.h"
#include "tree.h"                       //to create Tree objects
#include "ui_vars.h"
#include "vector_constants.h"
#include "xml_strings.h"                // for ToXML()
#include "event.h"                      // to create Event objects
#include "regiongammainfo.h"
#include "timemanager.h"

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

using namespace std;

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

ForceSummary::ForceSummary(ForceVec fvec, ForceParameters fparams, DataPack& dpack)
    : m_forcevec(fvec),
      m_startParameters(fparams),
      m_llikemles(dpack.GetNRegions()), // likelihoods at the MLEs per region.
      m_overallLlikemle()
{
    // first set the partition info
    LongVec1d partitioncounts;
    long ind = 0;
    ForceVec::iterator fit;
    for(fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
    {
        if ((*fit)->SetPartIndex(ind))
        {
            ++ind;
            long partCount = (*fit)->GetNPartitions();
            partitioncounts.push_back(partCount);
        }
        //Set up the identical groups;
        ULongVec2d newgroups = (*fit)->GetIdenticalGroupedParams();
        for (unsigned long gnum=0; gnum<newgroups.size(); gnum++)
        {
            m_identicalGroups.push_back(newgroups[gnum]);
        }
        //Set up the multiplicative groups;
        newgroups = (*fit)->GetMultiplicativeGroupedParams();
        for (unsigned long gnum=0; gnum<newgroups.size(); gnum++)
        {
            m_multiplicativeGroups.push_back(newgroups[gnum]);
        }
        //Set up the multiplicative multipliers;
        DoubleVec1d newmults = (*fit)->GetParamGroupMultipliers();
        m_multgroupmultiplier.insert(m_multgroupmultiplier.end(),newmults.begin(),newmults.end());
    }
    dpack.SetFinalPartitionCounts(partitioncounts);

} // ForceSummary constructor

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

ForceSummary::~ForceSummary()
{
    ForceVec::iterator fit;
    for (fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
        delete(*fit);
    m_forcevec.clear();
} // ForceSummary destructor

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

vector<Parameter> ForceSummary::GetAllParameters()
// a private function; public users should create a ParamVector
{
    ForceVec::iterator force;
    vector<Parameter> result;

    for (force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
    {
        vector<Parameter>& tempvec = (*force)->GetParameters();
        result.insert(result.end(), tempvec.begin(), tempvec.end());
    }

    return result;

} // GetAllParameters

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

void ForceSummary::SetAllParameters(const vector<Parameter>& src)
// a private function; public users should create a ParamVector
{
    ForceVec::iterator force;
    vector<Parameter>::const_iterator startparam = src.begin();
    vector<Parameter>::const_iterator endparam;
    for (force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
    {
        endparam = startparam + (*force)->GetNParams();
        vector<Parameter> result(startparam, endparam);
        (*force)->SetParameters(result);
        startparam = endparam;
    }
    //LS NOTE:  We might not be at the end of the param vector if there's a gamma, but
    // we throw that information away.  Or keep it somewhere else.  Or who knows what.

} // SetAllParameters

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

Tree *ForceSummary::CreateProtoTree() const
{
    Tree *tree;

    if (CheckForce(force_REC))
    {
        tree = new RecTree();
    }

    else
    {
        tree = new PlainTree();
    }

    tree->SetTreeTimeManager(CreateTimeManager());

    return tree;
}

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

TreeSummary* ForceSummary::CreateProtoTreeSummary() const
{
    bool shortform = true;

    // if growth is in effect, short form summaries cannot be used
    if (HasGrowth() || HasLogisticSelection())
    {
        shortform = false;
    }

    // if recombination is in effect, a recombinant summary must be
    // created
    TreeSummary* treesum;
    if (CheckForce(force_REC))
    {
        treesum = new RecTreeSummary();
    }
    else
    {
        treesum = new TreeSummary();
    }

    IntervalData& interval = treesum->GetIntervalData();

    ForceVec::const_iterator fit = m_forcevec.begin();
    for ( ; fit != m_forcevec.end(); ++fit)
    {
        Summary* sum;
        if (force_REC == (*fit)->GetTag() && !shortform)
        {
            if (CheckForce(force_DISEASE))
                sum = (*fit)->CreateSummary(interval, false);
            else
                sum = (*fit)->CreateSummary(interval, true);
        }
        else
            sum = (*fit)->CreateSummary(interval, shortform);
        if (sum)
        {
            treesum->AddSummary((*fit)->GetTag(), sum);
        }
    }

    return treesum;

} // CreateProtoTreeSummary

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

vector<Event*> ForceSummary::CreateEventVec() const
{
    vector<Event*> events;
    ForceVec::const_iterator fit(m_forcevec.begin());
    for( ; fit != m_forcevec.end(); ++fit)
    {
        vector<Event*> ev((*fit)->MakeEvents(*this));
        events.insert(events.end(),ev.begin(),ev.end());
    }

    for(fit = m_forcevec.begin(); fit != m_forcevec.end(); ++fit)
        (*fit)->ModifyEvents(*this,events);

    return events;
}

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

long ForceSummary::GetNLocalPartitionForces() const
{
    long nforces(0L);
    ForceVec::const_iterator fit(m_forcevec.begin());
    for( ; fit != m_forcevec.end(); ++fit)
    {
        if ((*fit)->IsLocalPartitionForce()) ++nforces;
    }
    return nforces;
}

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

ForceVec::const_iterator ForceSummary::GetForceByTag(force_type tag) const
{
    // returns an iterator to the element, or end() if none is found
    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->GetTag() == tag) return (it);
    }
    it = m_forcevec.end();
    return (it);

} // GetForceByTag

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

const StickForce& ForceSummary::GetStickForce() const
{
    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->IsStickForce())
        {
            StickForce* fp(dynamic_cast<StickForce*>(*it));
            return *fp;
        }
    }

    assert(false);  // asked for a stickforce when there wasn't one to be had!!

    return *(dynamic_cast<StickForce*>(*it)); // this is probably illegal,
    // but we should never be here
    // see assert above

} // GetStickForce

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

LongVec1d ForceSummary::GetLocalPartitionIndexes() const
{
    LongVec1d lpindexes;

    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->IsLocalPartitionForce())
        {
            lpindexes.push_back(dynamic_cast<PartitionForce*>(*it)->GetPartIndex());
        }
    }

    return lpindexes;

} // GetLocalPartitionIndexes

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

LongVec1d ForceSummary::GetNonLocalPartitionIndexes() const
{
    LongVec1d lpindexes;

    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->IsPartitionForce() && (!(*it)->IsLocalPartitionForce()))
        {
            lpindexes.push_back(dynamic_cast<PartitionForce*>(*it)->GetPartIndex());
        }
    }

    return lpindexes;

} // GetNonLocalPartitionIndexes

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

force_type ForceSummary::GetNonLocalPartitionForceTag() const
{
    if (CheckForce(force_MIG)) return force_MIG;
    if (CheckForce(force_DIVMIG)) return force_DIVMIG;
    assert(false);
    return force_NONE;  // should never happen
} // GetNonLocalPartitionForceTag

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

bool ForceSummary::IsMissingForce(force_type tag) const
{
    ForceVec::const_iterator it = GetForceByTag(tag);
    return (it == m_forcevec.end());
} // IsMissingForce

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

bool ForceSummary::AnyForcesOtherThanGrowCoal() const
{
    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->GetTag() != force_COAL && (*it)->GetTag() != force_GROW)
            return true;
    }

    return false;

} // AnyForcesOtherThanGrowCoal

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

bool ForceSummary::UsingStick() const
{
    ForceVec::const_iterator it;
    for (it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
    {
        if ((*it)->IsStickForce()) return true;
    }

    return false;
}

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

bool ForceSummary::CheckForce(force_type tag) const
{
    if (GetForceByTag(tag) == m_forcevec.end()) return false;
    return true;
}

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

long ForceSummary::GetMaxEvents(force_type tag) const
{
    ForceVec::const_iterator it = GetForceByTag(tag);
    assert(it != m_forcevec.end());
    return((*it)->GetMaxEvents());
} // GetMaxEvents

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

void ForceSummary::SetRegionMLEs(const ChainOut& chout, long region)
{
    // we do this by way of the ParamVector linearized form; it's easier

    ParamVector paramvec(false);
    vector<double> estimates = chout.GetEstimates().GetGlobalParameters();

    assert(estimates.size() == paramvec.size());  // need one estimate per parameter

    ParamVector::iterator param;
    vector<double>::const_iterator est = estimates.begin();

    for (param = paramvec.begin(); param != paramvec.end(); ++param, ++est)
    {
        if (param->IsValidParameter()) param->AddMLE(*est, region);
    }

    m_llikemles[region] = chout.GetLlikemle();

} // SetRegionMLEs

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

void ForceSummary::SetOverallMLE(const ChainOut& chout)
{
    ParamVector paramvec(false);
    vector<double> estimates = chout.GetEstimates().GetGlobalParameters();

    assert(estimates.size() == paramvec.size());  // need one estimate per parameter

    ParamVector::iterator param;
    vector<double>::const_iterator est = estimates.begin();

    for (param = paramvec.begin(); param != paramvec.end(); ++param, ++est)
    {
        if (param->IsValidParameter()) param->AddOverallMLE(*est);
    }

    m_overallLlikemle = chout.GetLlikemle();

} // SetOverallMLEs

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

MethodTypeVec1d ForceSummary::GetMethods(force_type tag) const
{
    ForceVec::const_iterator it = GetForceByTag(tag);
    assert(it != m_forcevec.end());
    return (*it)->GetMethods();

} // GetMethods

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

const vector<Epoch>* ForceSummary::GetEpochs() const
{
    if (!CheckForce(force_DIVERGENCE)) return NULL;
    ForceVec::const_iterator it = GetForceByTag(force_DIVERGENCE);
    const DivForce* divforce = dynamic_cast<const DivForce*>(*it);
    return divforce->GetEpochs();
} // GetEpochs

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

StringVec1d ForceSummary::ToXML(unsigned long nspaces) const
{
    StringVec1d xmllines;

    string line = MakeIndent(MakeTag(xmlstr::XML_TAG_FORCES),nspaces);
    xmllines.push_back(line);
    nspaces += INDENT_DEPTH;
    ForceVec::const_iterator force;
    for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
    {
        StringVec1d forcexml((*force)->ToXML(nspaces));
        xmllines.insert(xmllines.end(),forcexml.begin(),forcexml.end());
    }
    if (registry.GetRegionGammaInfo())  //checks a pointer value
    {
        StringVec1d gammaforcexml((*registry.GetRegionGammaInfo()).ToXML(nspaces));
        xmllines.insert(xmllines.end(),gammaforcexml.begin(),gammaforcexml.end());
    }
    nspaces -= INDENT_DEPTH;
    line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_FORCES),nspaces);
    xmllines.push_back(line);

    return xmllines;

} // ToXML

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

bool ForceSummary::ConstrainParameterValues(ForceParameters& fp) const
{
    bool parameters_constrained = false;
    const ParamVector pvec(true);
    ForceVec::const_iterator fit;
    DoubleVec1d params = fp.GetGlobalParameters();
    bool islog = false;
    for (unsigned long pnum = 0; pnum< pvec.size(); pnum++)
    {
        if (pvec[pnum].IsVariable())
        {
            force_type ftype = pvec[pnum].WhichForce();
            fit = GetForceByTag(ftype);
            double trimvalue = (*fit)->Truncate(params[pnum]);
            if (trimvalue != params[pnum])
            {
                string msg = "The parameter \"" + pvec[pnum].GetName() +
                    "\" maximized at " + ToString(params[pnum]) +
                    " but is being constrained to " + ToString(trimvalue) + ".\n";
                RunReport& runreport = registry.GetRunReport();
                runreport.ReportNormal(msg);
                SetParamWithConstraints(pnum, trimvalue, params, islog);
                //params[param] = trimvalue;
                parameters_constrained = true;
            }
        }
    }
    fp.SetGlobalParameters(params);

    return parameters_constrained;
} // ConstrainParameterValues

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

MethodTypeVec2d ForceSummary::Get2DMethods(force_type tag) const
{
    ForceVec::const_iterator it = GetForceByTag(tag);
    assert(it != m_forcevec.end());
    MethodTypeVec2d result = SquareOffVector((*it)->GetMethods());
    return result;
} // Get2DMethods

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

const ForceVec ForceSummary::GetPartitionForces() const
{
    ForceVec partitions;
    ForceVec::const_iterator it;
    for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
        if ((*it)->IsPartitionForce()) partitions.push_back(*it);

    return partitions;
} // GetPartitionForces

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

const ForceVec ForceSummary::GetLocalPartitionForces() const
{
    ForceVec partitions;
    ForceVec::const_iterator it;
    for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
        if ((*it)->IsLocalPartitionForce()) partitions.push_back(*it);

    return partitions;
} // GetLocalPartitionForces

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

ULongVec1d ForceSummary::GetParamsIdenticalGroupedWith(unsigned long pindex) const
{
    for (unsigned long gnum=0; gnum<m_identicalGroups.size(); gnum++)
    {
        for (unsigned long gpnum=0; gpnum<m_identicalGroups[gnum].size(); gpnum++)
        {
            if (m_identicalGroups[gnum][gpnum] == pindex)
            {
                return m_identicalGroups[gnum];
            }
        }
    }
    assert(false); //It's probably not a good idea to check this for *all*
    // pindices.
    ULongVec1d emptyvec;
    emptyvec.push_back(pindex);
    return emptyvec;
}

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

ULongVec1d ForceSummary::GetParamsMultiplicativelyGroupedWith(unsigned long pindex) const
{
    for (unsigned long gnum=0; gnum<m_multiplicativeGroups.size(); gnum++)
    {
        for (unsigned long gpnum=0; gpnum<m_multiplicativeGroups[gnum].size(); gpnum++)
        {
            if (m_multiplicativeGroups[gnum][gpnum] == pindex)
            {
                return m_multiplicativeGroups[gnum];
            }
        }
    }
    assert(false); //It's probably not a good idea to check this for *all*
    // pindices.
    ULongVec1d emptyvec;
    emptyvec.push_back(pindex);
    return emptyvec;
}

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

double ForceSummary::GetGroupMultiplier(unsigned long pindex) const
{
    for (unsigned long gnum=0; gnum<m_multiplicativeGroups.size(); gnum++)
    {
        for (unsigned long gpnum=0; gpnum<m_multiplicativeGroups[gnum].size(); gpnum++)
        {
            if (m_multiplicativeGroups[gnum][gpnum] == pindex)
            {
                return m_multgroupmultiplier[gnum];
            }
        }
    }
    assert(false); //It's probably not a good idea to check this for *all*
    // pindices.
    return FLAGDOUBLE;
}

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

vector<force_type> ForceSummary::GetForceTags() const
{
    vector<force_type> retVec;
    for (unsigned long fnum = 0; fnum<m_forcevec.size(); fnum++)
    {
        retVec.push_back(m_forcevec[fnum]->GetTag());
    }
    return retVec;
}

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

vector<long> ForceSummary::GetForceSizes() const
{
    vector<long> retVec;
    for (unsigned long fnum = 0; fnum<m_forcevec.size(); fnum++)
    {
        retVec.push_back(m_forcevec[fnum]->GetNParams());
    }
    return retVec;
}

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

void ForceSummary::SetParamWithConstraints(long pindex, double newval,
                                           DoubleVec1d& pvecnums, bool islog) const
{
    //Note that this function is sometimes called with 'pvecnums' equal to
    // *log* parameters, not always 'normal' parameters.

    const ParamVector pvec(true);
    // MDEBUG a good refactor would be to unify identical-group and multiplicative-group logic here
    ParamStatus mystatus = pvec[pindex].GetStatus();
    ULongVec1d groupmembers;
    double multiplier = 1.0;
    if (mystatus.Status() == pstat_identical_head)
    {
        groupmembers = GetParamsIdenticalGroupedWith(pindex);
    }
    if (mystatus.Status() == pstat_multiplicative_head)
    {
        groupmembers = GetParamsMultiplicativelyGroupedWith(pindex);
        multiplier = GetGroupMultiplier(pindex);
        if (islog) multiplier = log(multiplier);
    }
    mystatus.SetWithConstraints(pindex, newval, pvecnums, groupmembers, multiplier);
}

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

long ForceSummary::GetNParameters(force_type tag) const
{
    ForceVec::const_iterator it = GetForceByTag(tag);
    assert(it != m_forcevec.end());
    return (*it)->GetNParams();
} // GetNParameters

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

long ForceSummary::GetAllNParameters() const
{
    long total = 0;
    ForceVec::const_iterator it;
    for(it = m_forcevec.begin(); it != m_forcevec.end(); ++it)
        total += (*it)->GetNParams();
    return total;
} // ForceSummary::GetAllNParameters

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

vector<double> ForceSummary::GetModifiers(long posinparamvec) const
{
    const ParamVector paramvec(true);  // read-only copy
    vector<double> results;

    // diagnose whether we have fixed or percentile profiles.

    proftype profiletype = paramvec[posinparamvec].GetProfileType();

    //Note:  The profiling type must be consistent within a force, but can
    // change from force to force.  Also within a force, profiling can be turned
    // on or off.

    if (profiletype == profile_NONE) return results;  // no profiles, hence no modifiers

    verbosity_type verbosity = registry.GetUserParameters().GetVerbosity();

    if (profiletype == profile_PERCENTILE) // percentile
    {
        if (verbosity == CONCISE || verbosity == NONE)
            results = vecconst::percentiles_short;
        else
            results = vecconst::percentiles;
        return results;
    }
    else
    {                                   // fixed
        if (verbosity == CONCISE || verbosity == NONE)
        {
            results = vecconst::multipliers_short; //growth or not, just do this one.
        }
        else
        {
            if(paramvec[posinparamvec].IsForce(force_GROW))
            {
                results =  vecconst::growthmultipliers;
                DoubleVec1d gfixed = vecconst::growthfixed;
                for (unsigned long i=0; i<gfixed.size(); i++)
                    results.push_back(gfixed[i]);
                //There are two components to this vector--the first the normal
                // multipliers of growth, but the second actual values for growth
                // which we want to test.  We'll have to catch this later in
                // Analyzer::CalcProfileFixed (in profile.cpp)
            }
            else
                results = vecconst::multipliers;
        }
        return results;
    }
} // GetModifiers

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

DoubleVec1d ForceSummary::GetModifiersByForce(force_type tag) const
{
    const ParamVector paramvec(true);  // read-only copy
    for (unsigned long i=0; i<paramvec.size(); ++i)
        if (paramvec[i].IsValidParameter() && paramvec[i].IsForce(tag))
            if (paramvec[i].GetProfileType() != profile_NONE)
                return GetModifiers(i);
    return GetModifiers(0);
    //This assumes that all parameters of a given force have either the
    // same profiling type or none at all.
} // GetModifiersByForce

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

proftype ForceSummary::GetOverallProfileType() const
{
    unsigned long i;
    for (i = 0; i < m_forcevec.size(); ++i)
    {
        proftype mytype = m_forcevec[i]->SummarizeProfTypes();
        if (mytype != profile_NONE) return mytype;
    }

    return profile_NONE;

} // GetOverallProfileType

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

double ForceSummary::GetLowParameter(long posinparamvec) const
{
    const ParamVector paramvec(true);  // read-only copy
    force_type tag = paramvec[posinparamvec].WhichForce();
    if (force_REGION_GAMMA == tag)
    {
        const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
        if (!pRegionGammaInfo)
        {
            string msg = "ForceSummary::GetLowParameter() attempted to get the low ";
            msg += "parameter value for force \"" + ToString(tag);
            msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
            throw implementation_error(msg);
        }
        return pRegionGammaInfo->GetLowValue();
    }
    ForceVec::const_iterator fit = GetForceByTag(tag);
    if (m_forcevec.end() == fit)
    {
        string msg = "ForceSummary::GetLowParameter() failed to get a force object ";
        msg += "corresponding to tag \"" + ToString(tag) + ".\"";
        throw implementation_error(msg);
    }
    return (*fit)->GetLowVal();
}

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

double ForceSummary::GetHighParameter(long posinparamvec) const
{
    const ParamVector paramvec(true);  // read-only copy
    force_type tag = paramvec[posinparamvec].WhichForce();
    if (force_REGION_GAMMA == tag)
    {
        const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
        if (!pRegionGammaInfo)
        {
            string msg = "ForceSummary::GetHighParameter() attempted to get the high ";
            msg += "parameter value for force \"" + ToString(tag);
            msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
            throw implementation_error(msg);
        }
        return pRegionGammaInfo->GetHighValue();
    }
    ForceVec::const_iterator fit = GetForceByTag(tag);
    if (m_forcevec.end() == fit)
    {
        string msg = "ForceSummary::GetHighParameter() failed to get a force object ";
        msg += "corresponding to tag \"" + ToString(tag) + ".\"";
        throw implementation_error(msg);
    }
    return (*fit)->GetHighVal();
}

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

double ForceSummary::GetLowMult(long posinparamvec) const
{
    const ParamVector paramvec(true);  // read-only copy
    force_type tag = paramvec[posinparamvec].WhichForce();
    if (force_REGION_GAMMA == tag)
    {
        const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
        if (!pRegionGammaInfo)
        {
            string msg = "ForceSummary::GetLowMult() attempted to get the low ";
            msg += "multiplier value for force \"" + ToString(tag);
            msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
            throw implementation_error(msg);
        }
        return pRegionGammaInfo->GetLowMultiplier();
    }
    ForceVec::const_iterator fit = GetForceByTag(tag);
    if (m_forcevec.end() == fit)
    {
        string msg = "ForceSummary::GetLowMult() failed to get a force object ";
        msg += "corresponding to tag \"" + ToString(tag) + ".\"";
        throw implementation_error(msg);
    }
    return (*fit)->GetLowMult();
}

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

double ForceSummary::GetHighMult(long posinparamvec) const
{
    const ParamVector paramvec(true);  // read-only copy
    force_type tag = paramvec[posinparamvec].WhichForce();
    if (force_REGION_GAMMA == tag)
    {
        const RegionGammaInfo *pRegionGammaInfo = registry.GetRegionGammaInfo();
        if (!pRegionGammaInfo)
        {
            string msg = "ForceSummary::GetHighMult() attempted to get the high ";
            msg += "multiplier value for force \"" + ToString(tag);
            msg += ",\" but the necessary RegionGammaInfo was not found in the registry.";
            throw implementation_error(msg);
        }
        return pRegionGammaInfo->GetHighMultiplier();
    }
    ForceVec::const_iterator fit = GetForceByTag(tag);
    if (m_forcevec.end() == fit)
    {
        string msg = "ForceSummary::GetHighMult() failed to get a force object ";
        msg += "corresponding to tag \"" + ToString(tag) + ".\"";
        throw implementation_error(msg);
    }
    return (*fit)->GetHighMult();
}

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

long ForceSummary::GetPartIndex(force_type forcename) const
{
    assert(dynamic_cast<PartitionForce*>(*GetForceByTag(forcename)));
    return dynamic_cast<PartitionForce*>(*GetForceByTag(forcename))->
        GetPartIndex();
} // ForceSummary::GetPartIndex

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

TimeManager* ForceSummary::CreateTimeManager() const
{
    if (UsingStick())
    {
        if (HasGrowth())
        {
            if (HasSelection())
            {
                string msg = "ForceSummary::CreateTimeManager(), detected the ";
                msg += "simultaneous presence of the growth and selection ";
                msg += "forces.  If LAMARC has been updated to allow both of these ";
                msg += "forces to act at the same time, then ";
                msg += "ForceSummary::CreateTimeManager() needs to be updated ";
                msg += "accordingly.";
                throw implementation_error(msg);
            }
            return new ExpGrowStickTimeManager();
        }

        if (HasSelection())
            return new SelectStickTimeManager();

    }
    else
    {
        if (HasGrowth())
        {
            if (HasSelection())
            {
                string msg = "ForceSummary::CreateTimeManager(), detected the ";
                msg += "simultaneous presence of the growth and selection ";
                msg += "forces.  If LAMARC has been updated to allow both of these ";
                msg += "forces to act at the same time, then ";
                msg += "ForceSummary::CreateTimeManager() needs to be updated ";
                msg += "accordingly.";
                throw implementation_error(msg);
            }
            return new ExpGrowTimeManager();
        }

        if (HasLogisticSelection())
            return new LogSelectTimeManager();
    }

    return new ConstantTimeManager();

} // ForceSummary::CreateTimeManager

//------------------------------------------------------------------------------------
// WARNING warning--butt ugly code!  maybe add new classes of
// forces LocalPartitionForces and SexualForces (REC && CONVERSION)

LongVec1d ForceSummary::GetSelectedSites() const
{
    bool foundsex(false);
    ForceVec::const_iterator force;
    for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
    {
        if ((*force)->IsSexualForce())
        {
            foundsex = true;
            break;
        }
    }

    LongVec1d selsites;

    if (foundsex)
    {
        for(force = m_forcevec.begin(); force != m_forcevec.end(); ++force)
        {
            if ((*force)->IsLocalPartitionForce())
                selsites.push_back(dynamic_cast<LocalPartitionForce*>(*force)->GetLocalSite());
        }
    }

    return selsites;

} // ForceSummary::GetSelectedSites

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

bool ForceSummary::IsValidForceSummary() const
{
    // Some parameters must be present!
    if (m_forcevec.size() == 0) return false;

    // unsigned long npops = nparams;

    unsigned long ncoal = 0, nmig = 0, ndisease = 0, nrec = 0, ngrow = 0, ngamma = 0;
    unsigned long nlselect = 0, nepoch = 0;

    if (!CheckForce(force_COAL)) return false;
    ncoal = m_startParameters.GetGlobalThetas().size();

    if (CheckForce(force_MIG))
    {
        nmig = m_startParameters.GetMigRates().size();
        if (nmig < 2) return false;
    }

    if (CheckForce(force_DIVMIG))
    {
        nmig = m_startParameters.GetMigRates().size();
        if (nmig < 2) return false;
    }

    if (CheckForce(force_DISEASE))
    {
        ndisease = m_startParameters.GetDiseaseRates().size();
        if (ndisease < 2) return false;
    }

    //nmig and ndisease are set to 1 because I later take the square root to
    // determine the number of partitions for each.

    if (CheckForce(force_DIVERGENCE))
    {
        nepoch = m_startParameters.GetEpochTimes().size();
        unsigned long nmigparts = static_cast<long>(sqrt(static_cast<double>(nmig)));
        if (nepoch != (nmigparts - 1)/2) return false;   // truncation okay
    }

    if (CheckForce(force_REC))
    {
        nrec = m_startParameters.GetRecRates().size();
        if (nrec != 1) return false;
        // Until we implement distinct rec rates per population, at least.
    }

    if (CheckForce(force_GROW))
    {
        ngrow = m_startParameters.GetGrowthRates().size();
        if (ngrow != ncoal) return false;
    }

    if (CheckForce(force_LOGISTICSELECTION))
    {
        nlselect = m_startParameters.GetLogisticSelectionCoefficient().size();
        if (1 != nlselect) return false;
    }

    if (CheckForce(force_EXPGROWSTICK))
    {
        ngrow = m_startParameters.GetGrowthRates().size();
        if (ngrow != ncoal) return false;
    }

    unsigned long nmigparts;
    if (nmig)
    {
        nmigparts = static_cast<long>(sqrt(static_cast<double>(nmig)));
        if (nmigparts*nmigparts != nmig) return false;
    }
    else nmigparts = 1;

    unsigned long ndiseaseparts;
    if (ndisease)
    {
        ndiseaseparts = static_cast<long>(sqrt(static_cast<double>(ndisease)));
        if (ndiseaseparts*ndiseaseparts != ndisease) return false;
    }
    else ndiseaseparts = 1;

    if (ncoal != nmigparts * ndiseaseparts) return false;

    //Now check out a parameter vector and check to be sure the number of total
    // parameters is accurate, and that each parameter's index value is right.
    const ParamVector testvec(true);

    if (testvec.size() != ncoal + nmig + ndisease + nrec + nepoch + ngrow + ngamma + nlselect)
        return false;

    for (unsigned long index=0; index < testvec.size(); index++)
    {
        if (index != testvec[index].GetParamVecIndex()) return false;
    }

    // WARNING DEBUG should check validity of the parameterlist
    // (but I don't know how) -- Mary

    return true;
} // IsValidForceSummary

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

void ForceSummary::ValidateForceParamOrBarf(const ForceParameters& fp)
{
    DoubleVec1d params;
    vector<force_type> allforces;
    allforces.push_back(force_COAL);
    allforces.push_back(force_MIG);
    allforces.push_back(force_DIVMIG);
    allforces.push_back(force_DISEASE);
    allforces.push_back(force_REC);
    allforces.push_back(force_GROW);
    allforces.push_back(force_LOGISTICSELECTION);
    allforces.push_back(force_DIVERGENCE);

    for (vector<force_type>::iterator force = allforces.begin();
         force != allforces.end(); force++)
    {
        params = fp.GetGlobalParametersByTag((*force));
        if (CheckForce((*force)))
        {
            //This force is active, and should have a non-zero-sized vector in the forceparameters object.
            if (params.size() == 0)
            {
                string msg = "Error:  The force " + ToString((*force)) +
                    " was turned on in this run of LAMARC, but\n" +
                    " was not on in the LAMARC run that created this summary file.\n\n"
                    + "Re-run LAMARC with " + ToString((*force)) +
                    " turned off to accurately continue the old run.";
                throw data_error(msg);
            }
        }
        else
        {
            //This force is inactive, and should have a zero-sized vector in the forceparameters object.
            if (params.size() > 0)
            {
                string msg = "Error:  The force " + ToString((*force)) +
                    " was turned off in this run of LAMARC, but\n" +
                    " was on in the LAMARC run that created this summary file.\n\n" +
                    "Re-run LAMARC with " + ToString((*force)) +
                    " turned on to accurately continue the old run.";
                throw data_error(msg);
            }
        }
    }
}

//____________________________________________________________________________________