xref: /dports/biology/lamarc/lamarc-2.1.8/src/ui_vars/ui_vars_forces.cpp

// $Id: ui_vars_forces.cpp,v 1.90 2012/03/09 22:55:20 jmcgill Exp $

/*
 *  Copyright 2004  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 <algorithm>
#include <deque>
#include <map>
#include <vector>

#include "calculators.h"
#include "constants.h"          // for method_type
#include "defaults.h"
#include "errhandling.h"
#include "registry.h"
#include "stringx.h"
#include "ui_vars.h"
#include "ui_constants.h"
#include "ui_vars_forces.h"
#include "ui_vars_prior.h"
#include "ui_warnings.h"
#include "vectorx.h"

using std::deque;
using std::make_pair;
using std::map;
using std::vector;
using std::list;

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

UIVarsSingleForce::UIVarsSingleForce
(
    UIVars * myUIVars,
    force_type ftype,
    long int nparams,
    double defaultVal,
    method_type methodType,
    bool canBeSetUnset,
    bool isOn,
    long int eventMax,
    UIVarsPrior defaultPrior
    )
    :
    UIVarsComponent(myUIVars),
    m_numParameters(nparams),
    m_defaultValue(defaultVal),
    m_defaultMethodType(methodType),
    m_canSetOnOff(canBeSetUnset),
    m_onOff(isOn),
    m_maxEvents(eventMax),
    m_profileType(defaults::profileType),
    m_doProfile(nparams,true),
    m_userSuppliedStartValues(nparams,FLAGDOUBLE),
    // we CAN'T have user supplied start values at
    // the time this structure is created. So we
    // set them to FLAGDOUBLE. We should never access
    // these FLAGDOUBLE values because that would
    // require m_startValueMethodTypes to contain
    // method_USER, and we assert against that below
    m_startValueMethodTypes(nparams,methodType),
    m_pstatusValues(nparams,ParamStatus(pstat_unconstrained)),
    m_defaultPrior(defaultPrior),
    m_priors(nparams, defaultPrior),
    m_useDefaultPrior(nparams, true),
    m_ftype(ftype)
{
    long i;
    for (i = 0; i < nparams; ++i)
    {
        m_pstatusValues.push_back(ParamStatus(pstat_unconstrained));
    }
    assert(methodType        != method_USER);
    assert(m_defaultMethodType != method_USER);
    // must be a type that comes with a default value
    // (e.g. method_PROGRAMDEFAULT) or requires calculation
    // (method_FST or method_WATTERSON)
}

UIVarsSingleForce::UIVarsSingleForce
(
    UIVars * myUIVars,
    const UIVarsSingleForce& singleForce)
    :
    UIVarsComponent(myUIVars),
    m_numParameters                     (singleForce.m_numParameters),
    m_defaultValue                      (singleForce.m_defaultValue),
    m_defaultMethodType                 (singleForce.m_defaultMethodType),
    m_canSetOnOff                       (singleForce.m_canSetOnOff),
    m_onOff                             (singleForce.m_onOff),
    m_maxEvents                         (singleForce.m_maxEvents),
    m_profileType                       (singleForce.m_profileType),
    m_doProfile                         (singleForce.m_doProfile),
    m_userSuppliedStartValues           (singleForce.m_userSuppliedStartValues),
    m_startValueMethodTypes             (singleForce.m_startValueMethodTypes),
    m_groups                            (singleForce.m_groups),
    m_pstatusValues                     (singleForce.m_pstatusValues),
    m_defaultPrior                      (singleForce.m_defaultPrior),
    m_priors                            (singleForce.m_priors),
    m_useDefaultPrior                   (singleForce.m_useDefaultPrior),
    m_ftype                             (singleForce.m_ftype),
    m_calculatedStartValues             (singleForce.m_calculatedStartValues)
{
}

void UIVarsSingleForce::checkIndexValue(long int index) const
{
    // if you get an assert here, the menu or xml logic
    // isn't checking something it should be
    assert(index >= 0);
    assert(index < GetNumParameters());
    // However, allowing the user to put in parameter values in the Group
    // tag means we check here for that validity, and we need to
    // throw if it fails:
    if (index < 0)
    {
        throw data_error("Invalid (negative) parameter index in checkIndexValue.");
    }
    if (index >= GetNumParameters())
    {
        throw data_error("Invalid (too large) parameter index in checkIndexValue.");
    }
}

void UIVarsSingleForce::AssertOnIllegalStartMethod(method_type method)
{
    assert(method == method_PROGRAMDEFAULT);
    // assume that there are no calculation start method
    // types (such as method_FST or method_WATTERSON) for
    // an arbitrary force. Those that do have them will
    // override this method
}

// used to set a start method without providing a value (thus not
// suitable for use in Set
void UIVarsSingleForce::SetStartMethod(method_type method, long int index)
{
    checkIndexValue(index);
    AssertOnIllegalStartMethod(method);
    m_startValueMethodTypes[index] = method;

}

void UIVarsSingleForce::SetStartMethods(method_type method)
{
    for (int index = 0; index < GetNumParameters(); index++)
    {
        // use this method since we want to pick up
        // sanity checking done in virtual method
        // SetStartMethod
        SetStartMethod(method,index);
    }
}

string
UIVarsSingleForce::GetParamName(long int pindex)
{
    return GetConstUIVars().GetParamNameWithConstraint(m_ftype, pindex);
}

bool
UIVarsSingleForce::GetLegal() const
{
    return (GetOnOff() || m_canSetOnOff);
}

bool
UIVarsSingleForce::GetDoProfile(long int index) const
{
    checkIndexValue(index);
    ParamStatus mystatus = GetParamstatus(index);
    if (mystatus.Inferred())
    {
        if (GetConstUIVars().chains.GetDoBayesianAnalysis())
        {
            return true;
        }
        return m_doProfile[index];
    }
    else
        return false;
}

proftype UIVarsSingleForce::GetProfileType() const
{
    return m_profileType;
}

proftype UIVarsSingleForce::GetProfileType(long int index) const
{
    checkIndexValue(index);
    if( GetDoProfile(index))
    {
        return m_profileType;
    }
    else
    {
        return profile_NONE;
    }
}

string
UIVarsSingleForce::GetProfileTypeSummaryDescription() const
{
    proftype ptype = GetProfileType();
    paramlistcondition condition = GetParamListCondition();
    if(condition == paramlist_NO)
    {
        return ToString(condition) + "       ";
    }
    else
    {
        return ToString(ptype) +" (" +ToString(condition) +")";
    }
}

string
UIVarsSingleForce::GetPriorTypeSummaryDescription() const
{
    bool somelinear = false;
    bool somelog = false;
    for (unsigned long pnum=0; pnum<m_priors.size(); pnum++)
    {
        ParamStatus mystatus = GetParamstatus(pnum);
        if (mystatus.Inferred())
        {
            switch(GetPrior(pnum).GetPriorType())
            {
                case LINEAR:
                    somelinear = true;
                    break;
                case LOGARITHMIC:
                    somelog = true;
                    break;
            }
        }
    }
    if (somelinear && somelog)
    {
        return "(mixed linear/log)";
    }
    else if (somelinear && (!somelog))
    {
        return "(all linear)";
    }
    else if (somelog && (!somelinear))
    {
        return "(all logarithmic)";
    }
    else
    {
        return "(none eligible)";
    }
}

paramlistcondition
UIVarsSingleForce::GetParamListCondition() const
{
    bool seenOff = false;
    bool seenOn  = false;
    for (int index = 0; index < GetNumParameters(); index++)
    {
        // premature optimization being the root of all
        // evil, I've not chosen to break out of the loop when
        // seenOn and seenOff are both true, but you can change
        // it if you want -- ewalkup 8-19-2004
        if(GetParamValid(index))
        {
            if (GetDoProfile(index))
            {
                seenOn = true;
            }
            else
            {
                seenOff = true;
            }
        }
    }
    if(seenOn)
    {
        if(seenOff)
        {
            return paramlist_MIX;
        }
        else
        {
            return paramlist_YES;
        }
    }
    else
    {
        return paramlist_NO;
    }
}

void  UIVarsSingleForce::SetDoProfile(bool doIt, long int index)
{
    checkIndexValue(index);
    long gindex = ParamInGroup(index);
    ParamStatus mystatus = GetParamstatus(index);
    if (doIt && !mystatus.Varies())
    {
        // Do nothing--we can't turn on profiling for this parameter.
        if (GetParamValid(index))
        {
            string warning = "Warning:  " + GetParamName(index) +
                " may not be profiled because it is currently set to be "
                + ToString(mystatus.Status()) + ".";
            GetConstUIVars().GetUI()->AddWarning(warning);
        }
        return;
    }
    if (GetConstUIVars().chains.GetDoBayesianAnalysis())
    {
        if (doIt==false && mystatus.Inferred())
        {
            string warning = "Warning:  Profiling is always on for all valid parameters in a Bayesian analysis, since profiling takes no extra time.";
            GetConstUIVars().GetUI()->AddWarning(warning);
            return;
        }
    }
    m_doProfile[index] = doIt;
    if (gindex != FLAGLONG)
    {
        SetDoProfilesForGroup(doIt, gindex);
    }
}

void  UIVarsSingleForce::SetDoProfile(bool doIt)
{
    if (GetConstUIVars().chains.GetDoBayesianAnalysis())
    {
        if (doIt==false)
        {
            string warning = "Warning:  Profiling is always on for all valid parameters in a Bayesian analysis, since profiling takes no extra time.";
            GetConstUIVars().GetUI()->AddWarning(warning);
            return;
        }
    }
    for (int index = 0; index < GetNumParameters(); index++)
    {
        SetDoProfile(doIt,index);
    }
}

void UIVarsSingleForce::SetProfileType(proftype p)
{
    m_profileType = p;
}

ParamStatus UIVarsSingleForce::GetParamstatus(long int pindex)  const
{
    long int gindex = ParamInGroup(pindex);
    if (gindex != FLAGLONG)
    {
        ParamStatus mystatus = m_groups[gindex].first;
        if ((mystatus.Status() == pstat_identical) && (pindex == m_groups[gindex].second[0]))
        {
            return ParamStatus(pstat_identical_head);
            //MFIX probably needs changing, need to interface here
            //with finished UI design
        }
        else
            if ((mystatus.Status() == pstat_multiplicative) && (pindex == m_groups[gindex].second[0]))
            {
                return ParamStatus(pstat_multiplicative_head);
            }
            else
            {
                return mystatus;
            }
    }
    else return m_pstatusValues[pindex];
};

void UIVarsSingleForce::SetParamstatus(const ParamStatus& mystatus, long int index)
{
    long int gindex = ParamInGroup(index);
    if (gindex == FLAGLONG)
    {
        m_pstatusValues[index] = mystatus;
    }
    else
    {
        SetGroupParamstatus(mystatus.Status(), gindex);
    }
};

double UIVarsSingleForce::GetStartValue(long int pindex) const
{
    long int gindex = ParamInGroup(pindex);
    if (gindex == FLAGLONG)
    {
        return GetUngroupedStartValue(pindex);
    }
    double sum = 0;
    for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
    {
        sum += GetUngroupedStartValue(m_groups[gindex].second[gpindex]);
    }
    return (sum/m_groups[gindex].second.size());
}

double UIVarsSingleForce::GetUngroupedStartValue(long int index) const
{
    // start value might have been set by the user, or we might
    // have to grab it from the calculated start values

    double retval = 0.0;
    if (!GetParamstatus(index).Valid()) return 0.0;

    checkIndexValue(index);
    method_type thisMethod = m_startValueMethodTypes[index];

    switch(thisMethod)
    {
        case method_PROGRAMDEFAULT:
            retval = m_defaultValue;
            break;
        case method_USER:
            retval = m_userSuppliedStartValues[index];
            break;
        default:
        {
            // see if we have calculated values for the
            // start value method type for this parameter
            map<method_type,DoubleVec1d>::const_iterator it;
            it = m_calculatedStartValues.find(thisMethod);
            if(it == m_calculatedStartValues.end())
            {
                assert(false);
                throw implementation_error("Unable to calculate start values for " +
                                           ToString(thisMethod) + ".");
            }
            // we've populated the m_calculatedStartValues for
            // this method, so grab them and return the correct
            // one for this parameter
            retval = (it->second)[index];
            break;
        }
    }
    double minval = GetMinStartValue(index);
    double maxval = GetMaxStartValue(index);
    assert (maxval > minval);
    if (retval < minval)
    {
        return minval;
    }
    if (retval > maxval)
    {
        return maxval;
    }
    return retval;
}

DoubleVec1d
UIVarsSingleForce::GetStartValues() const
{
    DoubleVec1d retValue(GetNumParameters(),FLAGDOUBLE);
    // the user may have calculated start values and
    // then hand-edited some of them, so we need to
    // populate a vector with the individual values
    for(size_t i = 0; i < retValue.size(); i++)
    {
        retValue[i] = GetStartValue(i);
    }
    return retValue;
}

double UIVarsSingleForce::GetMinStartValue(long int pindex) const
{
    ParamStatus mystatus = GetParamstatus(pindex);
    if (!mystatus.Valid()) return 0;
    if (GetConstUIVars().chains.GetDoBayesianAnalysis())
    {
        if (mystatus.Status() != pstat_constant)
        {
            return GetPrior(pindex).GetLowerBound();
        }
    }
    switch(m_ftype)
    {
        case force_COAL:
            return defaults::minTheta;
        case force_MIG:
        case force_DIVMIG:
            return defaults::minMigRate;
        case force_DISEASE:
            return defaults::minDiseaseRate;
        case force_REC:
            return defaults::minRecRate;
        case force_EXPGROWSTICK:
        case force_GROW:
            return defaults::minGrowRate;
        case force_REGION_GAMMA:
            return defaults::minGammaOverRegions;
            break;
        case force_LOGSELECTSTICK:
        case force_LOGISTICSELECTION:
            return defaults::minLSelectCoeff;
            break;
        case force_DIVERGENCE:
            return defaults::minEpoch;
            break;
        default:
            assert(false);              //Uncaught force type.
    }
    return 0.0;  //To prevent compiler warning.
}

double UIVarsSingleForce::GetMaxStartValue(long int pindex) const
{
    ParamStatus mystatus = GetParamstatus(pindex);
    if (!mystatus.Valid()) return 0;
    if (GetConstUIVars().chains.GetDoBayesianAnalysis())
    {
        if (mystatus.Status() != pstat_constant)
        {
            return GetPrior(pindex).GetUpperBound();
        }
    }
    switch(m_ftype)
    {
        case force_COAL:
            return defaults::maxTheta;
        case force_MIG:
        case force_DIVMIG:
            return defaults::maxMigRate;
        case force_DISEASE:
            return defaults::maxDiseaseRate;
        case force_REC:
            return defaults::maxRecRate;
        case force_EXPGROWSTICK:
        case force_GROW:
            return defaults::maxGrowRate;
        case force_REGION_GAMMA:
            return defaults::maxGammaOverRegions;
            break;
        case force_LOGSELECTSTICK:
        case force_LOGISTICSELECTION:
            return defaults::maxLSelectCoeff;
            break;
        case force_DIVERGENCE:
            return defaults::maxEpoch;
            break;
        default:
            assert(false);              //Uncaught force type.
    }
    return 1000.0; // To prevent compiler warning.
}

void UIVarsSingleForce::SetUserStartValue(double startVal, long int index)
{
    checkIndexValue(index);
    bool isinvalid = !GetParamstatus(index).Valid();
    if (isinvalid)
    {
        if (startVal != 0)
        {
            string warning = "Warning:  " + GetParamName(index) +
                " ignores any attempt to set the starting value, since it has been set invalid.";
            GetConstUIVars().GetUI()->AddWarning(warning);
        }
        return;
    }
    double minval = GetMinStartValue(index);
    double maxval = GetMaxStartValue(index);
    assert ((maxval > minval) || isinvalid);
    if (startVal < minval)
    {
        startVal = minval;
        string warning = "Warning:  the minimum value for " + GetParamName(index) +
            " is " + ToString(minval) + ":  setting the starting value there.";
        GetConstUIVars().GetUI()->AddWarning(warning);
    }
    if (startVal > maxval)
    {
        startVal = maxval;
        string warning = "Warning:  the maximum value for " + GetParamName(index) +
            " is " + ToString(maxval) + ":  setting the starting value there.";
        GetConstUIVars().GetUI()->AddWarning(warning);
    }
    long int gindex = ParamInGroup(index);
    if (gindex != FLAGLONG)
    {
        for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
        {
            long int pindex = m_groups[gindex].second[gpindex];
            m_userSuppliedStartValues[pindex]=startVal;
            m_startValueMethodTypes[pindex] = method_USER;
        }
    }
    else
    {
        m_userSuppliedStartValues[index]=startVal;
        m_startValueMethodTypes[index] = method_USER; // don't use SetStartMethod
        // here -- it can block use
        // of method_USER since
        // the UI shouldn't set
        // method_USER without
        // going through this method
    }
}

void UIVarsSingleForce::SetUserStartValues(double startVal)
{
    for (int pindex = 0; pindex < GetNumParameters(); pindex++)
    {
        if (!GetParamstatus(pindex).Valid())
        {
            SetUserStartValue(startVal,pindex);
        }
    }
}

method_type UIVarsSingleForce::GetStartMethod(long int index) const
{
    checkIndexValue(index);
    return m_startValueMethodTypes[index];
}

void UIVarsSingleForce::SetTrueValue(double trueVal, long int index)
{
    checkIndexValue(index);
    if (!GetParamstatus(index).Valid())
    {
        if (trueVal != 0)
        {
            string warning = "Warning:  " + GetParamName(index) +
                " ignores any attempt to set the starting value, since it has been set invalid.";
            GetConstUIVars().GetUI()->AddWarning(warning);
        }
        return;
    }
    if (m_truevalues.empty())
    {
        m_truevalues.resize(GetNumParameters());
    }
    m_truevalues[index] = trueVal;
}

double UIVarsSingleForce::GetTrueValue(long int index) const
{
    checkIndexValue(index);
    if (m_truevalues.empty()) return FLAGDOUBLE; // JDEBUG
    else return m_truevalues[index];
}

const UIVarsPrior& UIVarsSingleForce::GetPrior(long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return m_defaultPrior;
    }
    checkIndexValue(pindex);
    if (GetUseDefaultPrior(pindex))
    {
        return m_defaultPrior;
    }
    return m_priors[pindex];
};

priortype UIVarsSingleForce::GetPriorType(long int pindex) const
{
    checkIndexValue(pindex);
    return GetPrior(pindex).GetPriorType();
};

double UIVarsSingleForce::GetLowerBound(long int pindex) const
{
    checkIndexValue(pindex);
    return GetPrior(pindex).GetLowerBound();
};

double UIVarsSingleForce::GetUpperBound(long int pindex) const
{
    checkIndexValue(pindex);
    return GetPrior(pindex).GetUpperBound();
};

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
long UIVarsSingleForce::GetRelativeSampling(long int pindex) const
{
    checkIndexValue(pindex);
    return GetPrior(pindex).GetRelativeSampling();
};
#endif

void UIVarsSingleForce::SetPriorType(priortype ptype, long int pindex)
{
    checkIndexValue(pindex);
    long int gindex = ParamInGroup(pindex);
    if (gindex == FLAGLONG)
    {
        SetUngroupedPriorType(ptype, pindex);
    }
    else for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
         {
             SetUngroupedPriorType(ptype, m_groups[gindex].second[gpindex]);
         }
}

void UIVarsSingleForce::SetUngroupedPriorType(priortype ptype, long int pindex)
{
    if (GetUseDefaultPrior(pindex))
    {
        m_useDefaultPrior[pindex] = false;
        m_priors[pindex] = m_defaultPrior;
    }
    m_priors[pindex].SetPriorType(ptype);
}

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
void UIVarsSingleForce::SetRelativeSampling(long rate, long int pindex)
{
    checkIndexValue(pindex);
    long int gindex = ParamInGroup(pindex);
    if (gindex == FLAGLONG)
    {
        SetUngroupedRelativeSampling(rate, pindex);
    }
    else
        for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
        {
            SetUngroupedRelativeSampling(rate, m_groups[gindex].second[gpindex]);
        }
}

void UIVarsSingleForce::SetUngroupedRelativeSampling(long rate, long int pindex)
{
    checkIndexValue(pindex);
    if (GetUseDefaultPrior(pindex))
    {
        m_useDefaultPrior[pindex] = false;
        m_priors[pindex] = m_defaultPrior;
    }
    m_priors[pindex].SetRelativeSampling(rate);
}

#endif

void UIVarsSingleForce::SetLowerBound(double bound, long int pindex)
{
    checkIndexValue(pindex);
    long int gindex = ParamInGroup(pindex);
    if (gindex == FLAGLONG)
    {
        SetUngroupedLowerBound(bound, pindex);
    }
    else
        for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
        {
            SetUngroupedLowerBound(bound, m_groups[gindex].second[gpindex]);
        }
}

void UIVarsSingleForce::SetUngroupedLowerBound(double bound, long int pindex)
{
    checkIndexValue(pindex);
    if (GetUseDefaultPrior(pindex))
    {
        m_useDefaultPrior[pindex] = false;
        m_priors[pindex] = m_defaultPrior;
    }
    m_priors[pindex].SetLowerBound(bound);
}

void UIVarsSingleForce::SetUpperBound(double bound, long int pindex)
{
    checkIndexValue(pindex);
    long int gindex = ParamInGroup(pindex);
    if (gindex == FLAGLONG)
    {
        SetUngroupedUpperBound(bound, pindex);
    }
    else for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
         {
             SetUngroupedUpperBound(bound, m_groups[gindex].second[gpindex]);
         }
}

void UIVarsSingleForce::SetUngroupedUpperBound(double bound, long int pindex)
{
    checkIndexValue(pindex);
    if (GetUseDefaultPrior(pindex))
    {
        m_useDefaultPrior[pindex] = false;
        m_priors[pindex] = m_defaultPrior;
    }
    m_priors[pindex].SetUpperBound(bound);
}

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

const UIVarsPrior& UIVarsSingleForce::GetDefaultPrior() const
{
    return m_defaultPrior;
};

bool UIVarsSingleForce::GetUseDefaultPrior(long int pindex) const
{
    checkIndexValue(pindex);
    return m_useDefaultPrior[pindex];
};

priortype UIVarsSingleForce::GetDefaultPriorType() const
{
    return m_defaultPrior.GetPriorType();
};

double UIVarsSingleForce::GetDefaultLowerBound() const
{
    return m_defaultPrior.GetLowerBound();
};

double UIVarsSingleForce::GetDefaultUpperBound() const
{
    return m_defaultPrior.GetUpperBound();
};

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
long UIVarsSingleForce::GetDefaultRelativeSampling() const
{
    return m_defaultPrior.GetRelativeSampling();
};
#endif

void UIVarsSingleForce::SetUseDefaultPrior(bool val, long int pindex)
{
    checkIndexValue(pindex);
    m_useDefaultPrior[pindex] = val;
};

void UIVarsSingleForce::SetUseAllDefaultPriors()
{
    for (unsigned long int pindex = 0; pindex < m_useDefaultPrior.size(); pindex++)
    {
        m_useDefaultPrior[pindex] = true;
    }
};

void UIVarsSingleForce::SetDefaultPriorType(priortype ptype)
{
    m_defaultPrior.SetPriorType(ptype);
}

void UIVarsSingleForce::SetDefaultLowerBound(double bound)
{
    m_defaultPrior.SetLowerBound(bound);
}

void UIVarsSingleForce::SetDefaultUpperBound(double bound)
{
    m_defaultPrior.SetUpperBound(bound);
}

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
void UIVarsSingleForce::SetDefaultRelativeSampling(long rate)
{
    m_defaultPrior.SetRelativeSampling(rate);
}
#endif

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

void
UIVarsSingleForce::SetMaxEvents(long int events)
{
    if (m_ftype != force_GROW && m_ftype != force_COAL && events <= 0)
    {
        throw data_error("The maximum number of events for force " + ToString(m_ftype) + " must be greater than zero.");
    }
    m_maxEvents = events;
}

void
UIVarsSingleForce::SetOnOff(bool onOffVal)
{
    assert( m_canSetOnOff || ( onOffVal == m_onOff ) );
    m_onOff = onOffVal;
}

bool
UIVarsSingleForce::GetParamValid(long int id) const
{
    return (GetOnOff() && (id >= 0) && (id < m_numParameters));
}

bool
UIVarsSingleForce::GetParamUnique(long int id) const
{
    long int gindex = ParamInGroup(id);
    if (gindex == FLAGLONG) return true;
    if (m_groups[gindex].second[0] == id) return true;
    return false;
}

bool UIVarsSingleForce::SomeVariableParams() const
{
    for (long int pnum = 0; pnum < m_numParameters; pnum++)
    {
        if (GetParamstatus(pnum).Varies()) return true;
    }
    return false;  // nothing ever varied
}

void UIVarsSingleForce::AddGroup(ParamStatus mystatus, LongVec1d indices)
{
    assert(indices.size() > 1);

    //Make a new ParamGroup and put it into m_groups.  Don't just put 'indices'
    // into it because we need to check each entry in it.
    LongVec1d blanklist;
    ParamGroup newgroup = make_pair(mystatus,blanklist);
    m_groups.push_back(newgroup);
    long int gindex = m_groups.size()-1; //<-OK because we just added one.

    //AddGroup is at the interface between the user and the program, whether
    // through the menu or through the XML. As such, the indexes used start
    // at one instead of zero--this is where that is changed for the program.
    for (unsigned long int vecindex = 0; vecindex < indices.size(); vecindex++)
    {
        AddParamToGroup(indices[vecindex] - 1, gindex);
    }
}

void UIVarsSingleForce::AddParamToGroup(long int pindex, long int gindex)
{
    checkGIndexValue(gindex);
    checkIndexValue(pindex);
    if (ParamInGroup(pindex) != FLAGLONG)
    {
        string e = "Invalid parameter indexes in group constraints for force "
            + ToString(m_ftype)
            + ":  The same parameter index is included in more than one group.";
        throw data_error(e);
    }
    if (m_groups[gindex].second.size() > 0)
    {
        //copy the prior to the new parameter.
        long int onegroupparam = m_groups[gindex].second[0];
        m_useDefaultPrior[pindex] = m_useDefaultPrior[onegroupparam];
        m_priors[pindex] = m_priors[onegroupparam];
    }
    m_groups[gindex].second.push_back(pindex);
    std::sort(m_groups[gindex].second.begin(), m_groups[gindex].second.end());
}

void UIVarsSingleForce::AddParamToNewGroup(long int pindex)
{
    checkIndexValue(pindex);
    if (ParamInGroup(pindex) != FLAGLONG)
    {
        throw data_error("Cannot add a parameter to a new group because it is already in a group.");
    }
    LongVec1d pindices;
    pindices.push_back(pindex);
    ParamGroup newgroup = make_pair(defaults::groupPstat, pindices);
    m_groups.push_back(newgroup);
}

void UIVarsSingleForce::RemoveParamIfInAGroup(long int pindex)
{
    long int gnum = ParamInGroup(pindex);
    if (gnum != FLAGLONG)
    {
        for (vector<long int>::iterator gpindex = m_groups[gnum].second.begin();
             gpindex != m_groups[gnum].second.end(); )
        {
            if ((*gpindex) == pindex)
            {
                ParamStatus old_pstat = GetGroupParamstatus(gnum);
                gpindex = m_groups[gnum].second.erase(gpindex);
                if (m_groups[gnum].second.size() == 0)
                {
                    //Delete the group--it's empty.
                    m_groups.erase(m_groups.begin() + gnum);
                }
                if (old_pstat.Grouped())
                {
                    SetParamstatus(ParamStatus(pstat_unconstrained), pindex);
                }
                else
                {
                    SetParamstatus(old_pstat, pindex);
                }
                return;
            }
            else
            {
                gpindex++;
            }
        }
    }
}

void UIVarsSingleForce::SetGroupParamstatus(ParamStatus pstat, long int gindex)
{
    // MDEBUG not correct for multiplicative!
    if (pstat.Status() == pstat_unconstrained)
    {
        pstat = ParamStatus(pstat_identical);
        GetConstUIVars().GetUI()->AddWarning("Warning:  Groups may not be set 'unconstrained'; setting to 'identical' instead.");
    }
    checkGIndexValue(gindex);
    m_groups[gindex].first = pstat;
}

ParamStatus UIVarsSingleForce::GetGroupParamstatus(long int gindex) const
{
    checkGIndexValue(gindex);
    return m_groups[gindex].first;
}

LongVec1d UIVarsSingleForce::GetGroupParamList(long int gindex) const
{
    checkGIndexValue(gindex);
    return m_groups[gindex].second;
}

long int UIVarsSingleForce::ParamInGroup(long int pindex) const
{
    checkIndexValue(pindex);
    for (long int gindex = 0; gindex < static_cast<long int>(m_groups.size()); gindex++)
    {
        for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
        {
            if (pindex == m_groups[gindex].second[gpindex])
            {
                return static_cast<long int>(gindex);
            }
        }
    }
    return FLAGLONG;
}

void UIVarsSingleForce::SetDoProfilesForGroup(bool doIt, long int gindex)
{
    //Since this is a private function that runs from SetDoProfile(), all the
    // error checking has already been done, and we need only set the
    // appropriate values.
    checkGIndexValue(gindex);
    for (unsigned long int gpindex = 0; gpindex < m_groups[gindex].second.size(); gpindex++)
    {
        m_doProfile[m_groups[gindex].second[gpindex]] = doIt;
    }
}

void UIVarsSingleForce::checkGIndexValue(long int gindex) const
{
    if (gindex==FLAGLONG) return;
    assert (gindex >= 0);
    if (gindex<0)
        throw data_error("Error:  group index value less than zero.");
    assert(gindex < static_cast<long int>(m_groups.size()));
    if (gindex >= static_cast<long int>(m_groups.size()))
        throw data_error("Error:  group index value greater than the number of groups.");
}

bool UIVarsSingleForce::AreGroupsValid() const
{
    LongVec1d allpindices;
    for (unsigned long int gindex = 0; gindex < m_groups.size(); gindex++)
    {
        if (m_groups[gindex].second.size() < 2) return false;
        allpindices.insert(allpindices.end(), m_groups[gindex].second.begin(), m_groups[gindex].second.end());
    }
    std::sort(allpindices.begin(), allpindices.end());
    for (unsigned long int gpindex=1; gpindex<allpindices.size(); gpindex++)
    {
        if (allpindices[gpindex] == allpindices[gpindex-1])
            return false;
    }
    return true;
}

void UIVarsSingleForce::FixGroups()
{
    //First, we remove any parameter that existed in more than one group. This
    // shouldn't ever be the case due to checks in AddParamToGroup, but hey.
    // First we have to find them again.
    LongVec1d allpindices;
    for (unsigned long int gindex = 0; gindex < m_groups.size(); gindex++)
    {
        allpindices.insert(allpindices.end(), m_groups[gindex].second.begin(), m_groups[gindex].second.end());
    }
    std::sort(allpindices.begin(), allpindices.end());
    for (unsigned long int gpindex=1; gpindex<allpindices.size(); gpindex++)
    {
        if (allpindices[gpindex] == allpindices[gpindex-1])
        {
            RemoveParamIfInAGroup(allpindices[gpindex]);
            //This removes the param from all groups.  It's not the ideal solution,
            // and the real solution is to never get in this situation in the first
            // place.  Which should be covered in AddParamToGroup.  So we're set.
        }
    }

    //Now, any parameter status flags that are 'joint' or 'identical' but do
    // not correspond with any group are set to 'unconstrained'.
    for (long int pindex = 0; pindex < m_numParameters; pindex++)
    {
        ParamStatus mystatus = GetParamstatus(pindex);
        if (mystatus.Grouped())
        {
            if (ParamInGroup(pindex) == FLAGLONG)
            {
                SetParamstatus(ParamStatus(pstat_unconstrained), pindex);
            }
        }
        //Also, any truly invalid parameters (i.e. the migration diagonals) may
        // not be in a group, nor may they be set to some status other than
        // 'invalid'.
        if (!GetParamValid(pindex))
        {
            RemoveParamIfInAGroup(pindex);
            SetParamstatus(ParamStatus(pstat_invalid), pindex);
        }
    }

    //Now, any group with only one parameter in it is deleted, setting the
    // corresponding ParamStatus appropriately.
    for (vector<ParamGroup>::iterator giter = m_groups.begin();
         giter != m_groups.end();)      // <- iterator incremented below.
    {
        if ((*giter).second.size() == 0)
        {
            giter = m_groups.erase(giter);
        }
        else if ((*giter).second.size() == 1)
        {
            long int pindex = (*giter).second[0];
            ParamStatus mystatus((*giter).first);
            if (mystatus.Grouped())
            {
                SetParamstatus(ParamStatus(pstat_unconstrained), pindex);
            }
            else
            {
                SetParamstatus(ParamStatus((*giter).first), pindex);
            }
            giter = m_groups.erase(giter);
        }
        else
        {
            ++giter;
        }
    }
}

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

UIVars2DForce::UIVars2DForce
(   UIVars * myUIVars,
    force_type ftype,
    long int numPartitions,
    double defaultVal,
    method_type methodType,
    bool canBeSetUnset,
    bool isOn,
    long int eventMax,
    UIVarsPrior defaultPrior
    )
    :
    UIVarsSingleForce(
        myUIVars,
        ftype,
        numPartitions*numPartitions,
        defaultVal,
        methodType,
        canBeSetUnset,         // can't change original on/off setting
        isOn,
        eventMax,
        defaultPrior),
    m_npartitions(numPartitions)
{
    // EWFIX.P5 DIMENSIONS -- we need the 1D <=> 2D translation objects
    //Make the diagonals of the parameter statuses (stati?) invalid.
    for (long int diagonal = 1; diagonal <= numPartitions; diagonal++)
    {
        SetParamstatus(ParamStatus(pstat_invalid), ((diagonal-1)*numPartitions + diagonal)-1);
    }
}

UIVars2DForce::UIVars2DForce
(   UIVars * myUIVars, const UIVars2DForce & twoD )
    :
    UIVarsSingleForce(myUIVars,twoD),
    m_npartitions(twoD.m_npartitions)
{
}

void UIVars2DForce::AssertOnIllegalStartMethod(method_type method)
{
    assert( (method == method_FST) || (method == method_PROGRAMDEFAULT) );
}

bool
UIVars2DForce::GetParamValid(long int id) const
{
    bool possible = UIVarsSingleForce::GetParamValid(id);
    if(possible)
    {
        // NB I believe this means "is it off-diagonal"?
        // MDEBUG This may need to be overridden in DivMig, but I don't
        // know how as it requires the Epoch structure to determine which
        // off-diagonal entries are invalid
        const UIVars& vars = GetConstUIVars();
        long int numPops = vars.datapackplus.GetNPartitionsByForceType(m_ftype);
        long int intoIndex = id / numPops;
        long int fromIndex = id % numPops;
        return (intoIndex != fromIndex);
    }
    return false;
}

//AreZeroesValid checks to make sure that we haven't set so many migration
// events to zero (and constant) that we can no longer coalesce.  Note that
// setting a parameter to 'pstat_invalid' means that we are setting it to zero
// and constant, too.
//
// This function could make suggestions about what to fix, but it suggestions
// are unlikely to be biologically relevant for the populations involved,
// so it seems better to just force the user to make their own changes, or just
// remove all 0/constant constraints.
//
// If it is to make suggestions, it would need to be non-const.

bool
UIVars2DForce::AreZeroesValid() const
{
    //The goal of this algorithm is to find 'Rome', a mythical land from which
    // all of your current populations descended.  In a migration system with
    // no constraints, all partitions could be Rome, but we only need to find
    // one to satisfy our test.  If we find at least one Rome, we return true,
    // and if not, we return false.
    //
    // First, we make two lists that will hold partition indexes in them,
    // 'reached' and 'unreached'.  At the beginning, everything is in
    // 'unreached'.

    list<long int> reached, unreached;
    for (long int partind = 0; partind < m_npartitions; partind++)
        unreached.push_back(partind);

    // Now we seize a partition at random (0) to see if it's Rome.
    long int testRome = 0;

    // Now we loop until we return either true or false.
    while (true)
    {
        //First, we don't have to worry about getting to testRome.  We also
        // don't need to add testRome to 'reached', since that list is only
        // used to back away from testRome if need be.
        unreached.remove(testRome);

        //Check to see who we can get to in 'unreached' from Rome.  We only
        // test routes that go through the unreached list, since everything in
        // the 'reached' list has already been exhaustively searched.
        UpdateReachList(testRome, unreached, reached);

        //Check to see if we're done.
        if (unreached.size() == 0) return true;

        //Otherwise, we need a new Rome candidate.  If our system is valid, Rome
        // *must* be in the 'unreached' group.  (If it wasn't, our old testRome
        // would have to have also been a valid Rome:  If everyone can come from
        // Rome by definition, and people in Rome can come from testRome, testRome
        // is also a valid 'Rome'.)
        testRome = PickANewRome(testRome, unreached, reached);

        //If PickANewRome failed, our system is invalid--return false;
        if (testRome == FLAGLONG) return false;

        //Otherwise, we can go ahead and loop with our new Rome.  Before we do,
        // though, we can clear out our 'reached' list.  We can do this because
        // not only do we know that migrants from our new testRome can get to
        // everywhere in 'reached', we *also* know that nothing in 'reached' can
        // get to 'unreached'--testRome and the other 'unreached' entries are the
        // only possibilities left.
        //
        // So, we clear 'reached'.  This has the effect of making our next loop
        // still more efficient.
        reached.clear();
    }
}

void UIVars2DForce::UpdateReachList(long int testRome, list<long int>& unreached,
                                    list<long int>& reached) const
{
    list<long int> oneStepAway;
    oneStepAway.push_back(testRome);

    while (oneStepAway.size() > 0)
    {
        list<long int> twoStepsAway;
        for (list<long int>::iterator partfrom = oneStepAway.begin();
             partfrom != oneStepAway.end(); partfrom++)
        {
            for (list<long int>::iterator partto = unreached.begin();
                 partto != unreached.end(); partto++)
            {
                if (CanReach((*partfrom), (*partto)))
                {
                    twoStepsAway.push_back((*partto));
                }
            }
        }
        for (list<long int>::iterator stepit = oneStepAway.begin();
             stepit != oneStepAway.end(); stepit++ )
        {
            reached.push_back((*stepit));
            unreached.remove((*stepit));
        }
        oneStepAway = twoStepsAway;
    }
}

long int UIVars2DForce::PickANewRome(long int oldRome, list<long int> unreached,
                                     list<long int> reached) const
{
    list<long int> intermediates;
    intermediates.push_back(oldRome);
    reached.remove(oldRome);
    while (intermediates.size() > 0)
    {
        list<long int> newintermediates;
        for (list<long int>::iterator partto = intermediates.begin();
             partto != intermediates.end(); partto++)
        {
            //Check to see if we can get to anything in the 'unreached' list.
            // If so, that's our new Rome candidate, and we're done.
            for (list<long int>::iterator partfrom = unreached.begin();
                 partfrom != unreached.end(); partfrom++)
            {
                if (CanReach((*partfrom), (*partto)))
                {
                    return (*partfrom);
                }
            }
            //We didn't get to 'unreached' directly; try going through the 'reached'
            // list.  If we find anything, we can remove it from 'reached'.
            for (list<long int>::iterator partfrom = reached.begin();
                 partfrom != reached.end();)
            {
                if (CanReach((*partfrom), (*partto)))
                {
                    newintermediates.push_back(*partfrom);
                    reached.erase(partfrom++);
                }
                else
                {
                    partfrom++;
                }
            }
        }
        intermediates = newintermediates;
    }

    //We never found any way to get to oldRome from anything in unreached.
    // Hence, our network is invalid.
    return FLAGLONG;
}

bool UIVars2DForce::CanReach(long int partfrom, long int partto) const
{
    long int pindex = partto*m_npartitions + partfrom;
    if (GetParamstatus(pindex).Valid())
    {
        double pval = GetStartValue(pindex);
        // If the start value is zero, the populations are unconnected at least
        // for the first tree, and that is a problem even in a Bayesian run.
        if (pval == 0.0) return false;
        else return true;
    }
    else
        return false;
}

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

UIVarsCoalForce::UIVarsCoalForce(UIVars * myUIVars, long int numCrossPartitions)
    :   UIVarsSingleForce(
        myUIVars,
        force_COAL,
        numCrossPartitions,
        defaults::theta,        //
        defaults::thetaMethod,
        false,                  // Coal force can't be turned off ever
        true,                   // Coal force always ON
        defaults::coalEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_COAL,
                    defaults::priortypeTheta,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundTheta,
                    defaults::upperboundTheta)
        )
{
}

UIVarsCoalForce::UIVarsCoalForce(UIVars * myUIVars, const UIVarsCoalForce& coalForce)
    :
    UIVarsSingleForce(myUIVars,coalForce)
{
}

void UIVarsCoalForce::AssertOnIllegalStartMethod(method_type method)
{
    assert( (method == method_FST)          ||
            (method == method_WATTERSON)    ||
            (method == method_PROGRAMDEFAULT));
    if ((method==method_FST) && !(GetConstUIVars().forces.GetForceLegal(force_MIG)))
    {
        throw data_error("The FST method is only available for data sets with more than one population.");
    }
}

void
UIVarsCoalForce::FillCalculatedStartValues()
{
    // do FST values
    if(GetConstUIVars().datapackplus.GetNCrossPartitions() > 1)
    {
        DoubleVec1d estimates;
        std::deque<bool> isCalculated;
        const UIVars& vars = GetConstUIVars();
        DoubleVec1d regionalthetascalars(vars.datapackplus.GetEffectivePopSizes());
        DoubleVec2d regionalmurates(vars.datamodel.GetRelativeMuRates());
        ThetaFSTMultiregion(registry.GetDataPack(), regionalthetascalars, regionalmurates, estimates, isCalculated);
        m_calculatedStartValues.erase(method_FST);
        m_calculatedStartValues.insert(make_pair(method_FST,estimates));

        std::set<string> defaultedEstimates;
        for (long int pindex = 0; pindex < static_cast<long int>(isCalculated.size()); pindex++)
        {
            if (!isCalculated[pindex])
            {
                string name = GetParamName(pindex);
                defaultedEstimates.insert(name);
            }
        }
        if (defaultedEstimates.size())
        {
            std::set<string>::iterator onename;
            string warning = "Warning:  calculating FST estimates for ";
            for (onename = defaultedEstimates.begin(); onename != defaultedEstimates.end(); onename++)
            {
                warning += "\"" + *onename + "\", ";
            }
            warning += "is impossible due to the data for the populations involved.  If the FST method is invoked to obtain starting values for those parameters, defaults will be used instead.";
            GetConstUIVars().GetUI()->AddWarning(warning);
        }
    }
    // do Watterson values
    {
        DoubleVec1d estimates;
        std::deque<bool> isCalculated;
        const UIVars& vars = GetConstUIVars();
        DoubleVec1d regionalthetascalars(vars.datapackplus.GetEffectivePopSizes());
        DoubleVec2d regionalmurates(vars.datamodel.GetRelativeMuRates());
        ThetaWattersonMultiregion(registry.GetDataPack(), regionalthetascalars,
                                  regionalmurates, estimates, isCalculated);
        m_calculatedStartValues.erase(method_WATTERSON);
        m_calculatedStartValues.insert(make_pair(method_WATTERSON,estimates));

        std::set<string> defaultedEstimates;
        for (long int pindex = 0; pindex < static_cast<long int>(isCalculated.size()); pindex++)
        {
            if (!isCalculated[pindex])
            {
                string name = GetParamName(pindex);
                defaultedEstimates.insert(name);
            }
        }
        if (defaultedEstimates.size())
        {
            std::set<string>::iterator onename;
            string warning = "Warning:  calculating Watterson estimates for ";
            for (onename = defaultedEstimates.begin();
                 onename != defaultedEstimates.end(); onename++)
            {
                warning += "\"" + *onename + "\", ";
            }
            warning += "is impossible because one or more genetic regions have no data or only one data point in the corresponding population(s).  If the Watterson method is invoked to obtain starting values for theta, defaults will be used in that calculation.";
            GetConstUIVars().GetUI()->AddWarning(warning);
        }
    }
}

bool
UIVarsCoalForce::AreZeroesValid() const
{
    //No force can be invalid or zero for coalescence.  It can't even start at
    // zero, regardless of whether it's held there.
    for (long int pindex = 0; pindex<GetNumParameters(); pindex++)
    {
        if (!GetParamstatus(pindex).Valid())
        {
            return false;
        }
        if (GetStartValue(pindex) == 0.0)
        {
            return false;
        }
    }
    return true;
}

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

UIVarsMigForce::UIVarsMigForce (UIVars * myUIVars, long int numPopulations, bool onOrOff)
    :
    UIVars2DForce(
        myUIVars,
        force_MIG,
        numPopulations,
        defaults::migration,
        defaults::migMethod,
        false,                  // can't change original on/off setting
        onOrOff,
        defaults::migEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_MIG,
                    defaults::priortypeMig,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundMig,
                    defaults::upperboundMig)
        )
{
}

UIVarsMigForce::UIVarsMigForce
(   UIVars * myUIVars, const UIVarsMigForce & migForce)
    :
    UIVars2DForce(myUIVars,migForce)
{
}

void
UIVarsMigForce::FillCalculatedStartValues()
{
    DoubleVec1d estimates;
    std::deque<bool> isCalculated;
    const UIVars& vars = GetConstUIVars();
    DoubleVec2d regionalmurates(vars.datamodel.GetRelativeMuRates());
    MigFSTMultiregion(registry.GetDataPack(), regionalmurates,
                      estimates, isCalculated);
    m_calculatedStartValues.erase(method_FST);
    m_calculatedStartValues.insert(make_pair(method_FST,estimates));

    std::set<string> defaultedEstimates;
    for (long int pindex = 0; pindex < static_cast<long int>(isCalculated.size()); pindex++)
    {
        if (!isCalculated[pindex])
        {
            string name = GetParamName(pindex);
            defaultedEstimates.insert(name);
        }
    }
    if (defaultedEstimates.size())
    {
        std::set<string>::iterator onename;
        string warning = uiwarn::calcFST_0;
        for (onename = defaultedEstimates.begin(); onename != defaultedEstimates.end(); onename++)
        {
            warning += "\"" + *onename + "\", ";
        }
        warning += uiwarn::calcFST_1;
        GetConstUIVars().GetUI()->AddWarning(warning);
        //LS DEBUG:  This warning is always printed out when the user
        // first fires up LAMARC, and never afterwards, due to the fact
        // that FST is the default starting method for Migration.  If
        // this changes, the text of this method should probably change,
        // too.
    }
}

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

UIVarsDivMigForce::UIVarsDivMigForce (UIVars * myUIVars, long int numPopulations, bool onOrOff)
    :
    UIVars2DForce(
        myUIVars,
        force_DIVMIG,
        numPopulations,
        defaults::divMigration,
        defaults::divMigMethod,
        false,                  // can't change original on/off setting
        onOrOff,
        defaults::migEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_DIVMIG,
                    defaults::priortypeDivMig,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundDivMig,
                    defaults::upperboundDivMig)
        )
{
}

UIVarsDivMigForce::UIVarsDivMigForce
(   UIVars * myUIVars, const UIVarsDivMigForce & divmigForce)
    :
    UIVars2DForce(myUIVars,divmigForce)
{
}

void UIVarsDivMigForce::AssertOnIllegalStartMethod(method_type method)
{
    assert(method == method_PROGRAMDEFAULT); // can't use FST here
}

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

UIVarsDivergenceForce::UIVarsDivergenceForce(UIVars * myUIVars, long int nDivPopulations)
    :   UIVarsSingleForce(
        myUIVars,
        force_DIVERGENCE,
        (nDivPopulations-1)/2, // truncation deliberate
        defaults::epochtime,
        defaults::divMethod,
        // MDEBUG not sure about the values of these two bools
        true,                   // Divergence force can't be turned off ever
        false,                   // Divergence force always ON
        defaults::epochEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_DIVERGENCE,
                    defaults::priortypeEpoch,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundEpoch,
                    defaults::upperboundEpoch)
        )
{
    // the Divergence parameters, which are epoch times, must be specially marked
    long i;
    for (i = 0; i < (nDivPopulations-1)/2; ++i)
    {
        // JDEBUG, need to do something here? MakeParamAnEpochTime(i);
    }
}

UIVarsDivergenceForce::UIVarsDivergenceForce(UIVars * myUIVars, const UIVarsDivergenceForce& divForce)
    :
    UIVarsSingleForce(myUIVars,divForce)
{
    // copy epoch info
    newpops = divForce.newpops;
    ancestors = divForce.ancestors;
}

void UIVarsDivergenceForce::AssertOnIllegalStartMethod(method_type method)
{
    assert (method == method_PROGRAMDEFAULT);
}

bool
UIVarsDivergenceForce::AreZeroesValid() const
{
    //No parameter can be invalid or zero for divergence.  It can't even start at
    // zero, regardless of whether it's held there.
    for (long int pindex = 0; pindex<GetNumParameters(); pindex++)
    {
        if (!GetParamstatus(pindex).Valid())
        {
            return false;
        }
        if (GetStartValue(pindex) == 0.0)
        {
            return false;
        }
    }
    return true;
}

void
UIVarsDivergenceForce::AddNewPops(const vector<string>& newp)
{
    newpops.push_back(newp);
}

void
UIVarsDivergenceForce::AddAncestor(const string& anc)
{
    ancestors.push_back(anc);
}

vector<vector<string> >
UIVarsDivergenceForce::GetNewPops() const
{
    return newpops;
}

vector<string>
UIVarsDivergenceForce::GetAncestors () const
{
    return ancestors;
}

string
UIVarsDivergenceForce::GetAncestor (long index) const
{
    return ancestors[index];
}

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

UIVarsDiseaseForce::UIVarsDiseaseForce(
    UIVars * myUIVars,
    long int numDiseaseStates,
    bool isLegal)
    :
    UIVars2DForce(
        myUIVars,
        force_DISEASE,
        numDiseaseStates,
        defaults::disease,
        defaults::diseaseMethod,
        false,              // EWFIX.P5 DISEASE change false => isLegal later
        // when we figure out how to accomodate
        // change in number of theta and growth
        // params when disease forces are turned
        // on and off in the menu
        isLegal,
        defaults::diseaseEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_DISEASE,
                    defaults::priortypeDisease,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundDisease,
                    defaults::upperboundDisease)
        ),
    location(defaults::diseaseLocation)
{
    //LS NOTE:  When we get multiple diseases states in, this is the place to
    // set up the invalid and identical parameters (re: the tic-tac-toe setup
    // we talked about on Joe's blackboard)
    //
    /* For example:  say we have two traits, Aa, and XYZ.  You can set things up:

       AX  AY  AZ    aX  aY  aZ
       AX --  ++  ++    ++  --  --
       AY ++  --  ++    --  ++  --
       AZ ++  ++  --    --  --  ++

       aX ++  --  --    --  ++  ++
       aY --  ++  --    ++  --  ++
       aZ --  --  ++    ++  ++  --

       So all the ++'s are possible, and all the --'s are set 'invalid', like the
       diagonal.  Additionally, if A->a is constant, and not influenced by
       the current XYZ, you can set AX->aX, AY->aY, and AZ->aZ to be equal to
       each other, and likewise for the other sets of rates.
    */
}

UIVarsDiseaseForce::UIVarsDiseaseForce
(   UIVars * myUIVars, const UIVarsDiseaseForce & diseaseForce)
    :
    UIVars2DForce(myUIVars,diseaseForce),
    location(diseaseForce.location)
{
}

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

UIVarsRecForce::UIVarsRecForce(UIVars * myUIVars, bool canTurnOn)
    :
    UIVarsSingleForce(
        myUIVars,
        force_REC,
        1L,                     // only one param value for recombine
        defaults::recombinationRate,
        defaults::recMethod,
        canTurnOn,              // recombine illegal for unlinked data
        false,                  // recombine force not required when possible
        defaults::recEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_REC,
                    defaults::priortypeRec,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundRec,
                    defaults::upperboundRec)
        )
{
}

UIVarsRecForce::UIVarsRecForce(UIVars * myUIVars, const UIVarsRecForce& recForce)
    :
    UIVarsSingleForce(myUIVars,recForce)
{
}

bool
UIVarsRecForce::AreZeroesValid() const
{
    //Recombination can be zero.  You might as well not have it on if so,
    // but who are we to quibble?  Also, if we allow different populations to
    // have different rec values at some point, any of them could be zero.
    return true;
}

bool UIVarsRecForce::GetOnOff() const
{
    if (GetConstUIVars().traitmodels.AnyMappingAnalyses()) return true;
    return UIVarsSingleForce::GetOnOff();
}

void UIVarsRecForce::SetOnOff(bool onOffVal)
{
    if (GetConstUIVars().traitmodels.AnyMappingAnalyses() && !onOffVal)
    {
        throw data_error("Cannot turn off recombination when there are traits to be mapped.");
    }
    if (GetConstUIVars().datapackplus.AnySNPDataWithDefaultLocations() && onOffVal)
    {
        throw data_error("We cannot accurately estimate recombination as there is SNP data whose locations are unknown.  "
                         "If you want to estimate recombination, go back to the converter and use a map file to tell it where the sequenced SNPs are.  "
                         "See the documentation: the \"Length and Spacing Information with Segment Coordinates\" section in "
                         "genetic_map.html for how to do this within the GUI converter, or the \"Segments\" section "
                         "of converter_cmd.html for how to do this using the converter command file.");
    }
    UIVarsSingleForce::SetOnOff(onOffVal);
}

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

UIVarsRegionGammaForce::UIVarsRegionGammaForce(UIVars * myUIVars)
    :
    UIVarsSingleForce(
        myUIVars,
        force_REGION_GAMMA,
        1L,                     // only one param value for region gamma
        defaults::gammaOverRegions,
        defaults::recMethod,
        true,              // legal to turn on--ONLY for multi-region data
        false,             // this "force" is not required
        0,                 // no such thing as region gamma "events"
        UIVarsPrior(myUIVars->GetUI(),
                    force_REC,               //BUGBUG!!! no prior allowed!
                    defaults::priortypeRec,  //what to do???
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundRec,
                    defaults::upperboundRec)
        )
{
    if (GetConstUIVars().datapackplus.GetNumRegions() == 1)
    {
        m_canSetOnOff = false;
    }
}

UIVarsRegionGammaForce::UIVarsRegionGammaForce(UIVars * myUIVars, const UIVarsRegionGammaForce& regionGammaForce)
    :
    UIVarsSingleForce(myUIVars,regionGammaForce)
{
}

bool UIVarsRegionGammaForce::GetOnOff() const
{
    if (GetConstUIVars().datapackplus.GetNumRegions() == 1) return false;
    if (GetConstUIVars().chains.GetDoBayesianAnalysis()) return false;
    if (GetConstUIVars().forces.GetForceOnOff(force_GROW)) return false;
    return UIVarsSingleForce::GetOnOff();
}

void
UIVarsRegionGammaForce::SetOnOff(bool onOffVal)
{
    if (GetConstUIVars().datapackplus.GetNumRegions() == 1 && onOffVal)
    {
        throw data_error("Cannot analyze the region gamma when there is only one region to analyze.");
    }
    if (GetConstUIVars().chains.GetDoBayesianAnalysis() && onOffVal)
    {
        throw data_error("Cannot analyze the region gamma during a Bayesian analysis.");
    }
    if (GetConstUIVars().forces.GetForceOnOff(force_GROW) && onOffVal)
    {
        throw data_error("Cannot analyze the region gamma if you are simultaneously estimating the growth rate.");
        return;
    }
    UIVarsSingleForce::SetOnOff(onOffVal);
}

bool
UIVarsRegionGammaForce::AreZeroesValid() const
{
    if (!GetParamstatus(0).Valid())
        return false;
    if (GetStartValue(0) == 0.0)
        return false;

    return true;
}

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

UIVarsGrowForce::UIVarsGrowForce(UIVars * myUIVars, long int numCrossPartitions)
    :
    UIVarsSingleForce(
        myUIVars,
        force_GROW,
        numCrossPartitions,
        defaults::growth,
        defaults::growMethod,
        true,                   // growth force always legal to turn on
        false,                  // growth force not required when possible
        defaults::growEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_GROW,
                    defaults::priortypeGrowth,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundGrowth,
                    defaults::upperboundGrowth)
        ),
    growthType(defaults::growType)
{
}

UIVarsGrowForce::UIVarsGrowForce(UIVars * myUIVars, const UIVarsGrowForce& growForce)
    :
    UIVarsSingleForce(myUIVars,growForce),
    growthType(growForce.growthType)
{
}

bool
UIVarsGrowForce::AreZeroesValid() const
{
    //Any growth can be zero; even up to all of them.
    return true;
}

// We don't currently have any STAIRSTEP models, so the case statement
// is probably wrong for when we have them... Jon 2006/09/15
growth_scheme
UIVarsGrowForce::GetGrowthScheme() const
{
    switch(growthType)
    {
        case growth_CURVE:
            return growth_EXP;
            break;
        case growth_STICK:  // not sure what this case is, assuming
            // that we can use the same as STICKEXP for
            // now
        case growth_STICKEXP:
            return growth_EXP;
            break;
    }
    return growth_EXP;
}

void
UIVarsGrowForce::SetGrowthScheme(growth_scheme g)
{
    // this is currently a no-op, as soon as we get STAIRSTEP
    // we'll need something else, probably ickly complicated here
    // in order to account for interlocking nature with growth_type
    // management...refactor growth_type?  Jon 2006/09/15
}

force_type
UIVarsGrowForce::GetPhase2Type(force_type f) const
{
    if (growthType == growth_STICKEXP) return force_EXPGROWSTICK;

    return f;
}

//LS NOTE:  Cannot have a parallel 'GetOnOff' for growth here that checks
// whether gamma is on or off, because gamma checks if *growth* is on or off.
// If you try it, you'll get an infinite loop.  Yes, this has been... verified.

void UIVarsGrowForce::SetOnOff(bool onOffVal)
{
    if (GetConstUIVars().forces.GetForceOnOff(force_REGION_GAMMA))
    {
        throw data_error("Cannot estimate growth if you are simulataneously estimating the region gamma.");
        return;
    }
    UIVarsSingleForce::SetOnOff(onOffVal);
}

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

UIVarsLogisticSelectionForce::UIVarsLogisticSelectionForce
(   UIVars * myUIVars, long numCrossPartitions)
    :
    UIVarsSingleForce(
        myUIVars,
        force_LOGISTICSELECTION,
        1, // there's only one selection coefficient
        defaults::logisticSelection,
        defaults::lselectMethod,
        true,                   // always legal to turn on
        false,                  // but it's not required
        defaults::lselectEvents,
        UIVarsPrior(myUIVars->GetUI(),
                    force_LOGISTICSELECTION,
                    defaults::priortypeLSelect,
#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
                    defaults::samplingRate,
#endif
                    defaults::lowerboundLSelect,
                    defaults::upperboundLSelect)
        )
    ,selectionType(defaults::selectionType) // MCHECK major allele freq here?
{
    if (2 != numCrossPartitions)
        m_canSetOnOff = false;
}

UIVarsLogisticSelectionForce::UIVarsLogisticSelectionForce(UIVars * myUIVars,
                                                           const UIVarsLogisticSelectionForce& logisticSelectionForce) :
    UIVarsSingleForce(myUIVars, logisticSelectionForce)
{
}

bool
UIVarsLogisticSelectionForce::AreZeroesValid() const
{
    // Sure; zero means selection isn't acting.
    return true;
}

force_type
UIVarsLogisticSelectionForce::GetPhase2Type(force_type f) const
{
    if (selectionType == selection_STOCHASTIC) return force_LOGSELECTSTICK;

    return f;
}

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

bool
UIVarsForces::IsInactive::operator()(force_type ft) const
{
    const UIVarsSingleForce & theForce
        = m_vars_forces.getForceRegardlessOfLegality(ft);
    bool isInactive = !(theForce.GetOnOff());
    return isInactive;
}

UIVarsForces::IsInactive::IsInactive(const UIVarsForces& vars_forces)
    : m_vars_forces(vars_forces)
{
}

UIVarsForces::IsInactive::~IsInactive()
{
}

bool
UIVarsForces::IsIllegal::operator()(force_type ft) const
{
    const UIVarsSingleForce & theForce
        = m_vars_forces.getForceRegardlessOfLegality(ft);
    bool isIllegal = (!theForce.GetLegal());
    return isIllegal;
}

UIVarsForces::IsIllegal::IsIllegal(const UIVarsForces& vars_forces)
    : m_vars_forces(vars_forces)
{
}

UIVarsForces::IsIllegal::~IsIllegal()
{
}

UIVarsForces::UIVarsForces(UIVars * myUIVars,long int nCrossPartitions, long int nMigPopulations,
                           long int nDivPopulations, long int nDiseaseStates,bool canMeasureRecombination)
    :
    UIVarsComponent(myUIVars),
    coalForce(myUIVars,nCrossPartitions),
    diseaseForce(myUIVars,nDiseaseStates,(nDiseaseStates >= 2)),
    growForce(myUIVars,nCrossPartitions),
    migForce(myUIVars,nMigPopulations,(nMigPopulations >= 2)),
    divMigForce(myUIVars,nDivPopulations,(nDivPopulations > 2)),
    divForce(myUIVars,nDivPopulations),   // yes, integer division with truncation is desired
    recForce(myUIVars,canMeasureRecombination),
    regionGammaForce(myUIVars),
    logisticSelectionForce(myUIVars,nCrossPartitions)
{
    //requires that both datapackplus and datamodel be set up in uivars.
    FillCalculatedStartValues();
}

UIVarsForces::UIVarsForces(UIVars * myUIVars, const UIVarsForces& forcesRef)
    :
    UIVarsComponent(myUIVars),
    coalForce(myUIVars,forcesRef.coalForce),
    diseaseForce(myUIVars,forcesRef.diseaseForce),
    growForce(myUIVars,forcesRef.growForce),
    migForce(myUIVars,forcesRef.migForce),
    divMigForce(myUIVars,forcesRef.divMigForce),
    divForce(myUIVars,forcesRef.divForce),
    recForce(myUIVars,forcesRef.recForce),
    regionGammaForce(myUIVars,forcesRef.regionGammaForce),
    logisticSelectionForce(myUIVars,forcesRef.logisticSelectionForce)
{
}

UIVarsForces::~UIVarsForces()
{
}

void
UIVarsForces::FillCalculatedStartValues()
{
    coalForce.FillCalculatedStartValues();
    migForce.FillCalculatedStartValues();
}

LongVec1d UIVarsForces::GetForceSizes() const
{
    LongVec1d retVec;
    vector<force_type> forces = GetActiveForces();
    for (unsigned long int fnum = 0; fnum < forces.size(); fnum++)
    {
        retVec.push_back(getLegalForce(forces[fnum]).GetNumParameters());
    }
    return retVec;
}

ForceTypeVec1d
UIVarsForces::GetActiveForces() const
{
    ForceTypeVec1d returnVec = GetLegalForces();
    returnVec.erase(
        std::remove_if(returnVec.begin(),returnVec.end(),IsInactive(*this)),
        returnVec.end());
    return returnVec;
}

ForceTypeVec1d
UIVarsForces::GetPhase2ActiveForces() const
{
    ForceTypeVec1d tempVec = GetActiveForces();
    ForceTypeVec1d returnVec;
    for (unsigned long fnum=0; fnum<tempVec.size(); fnum++)
    {
        returnVec.push_back(getLegalForce(tempVec[fnum]).GetPhase2Type(tempVec[fnum]));
    }
    return returnVec;
}

ForceTypeVec1d
UIVarsForces::GetLegalForces() const
{
    ForceTypeVec1d returnVec = GetPossibleForces();
    returnVec.erase(
        std::remove_if(returnVec.begin(),returnVec.end(),IsIllegal(*this)),
        returnVec.end());
    return returnVec;
}

ForceTypeVec1d
UIVarsForces::GetPossibleForces() const
{
    //This is where the so-called 'canonical order' of forces in LAMARC
    // gets set up.  If we add a new force, add it here appropriately.
    ForceTypeVec1d returnVec;
    returnVec.push_back(force_COAL);
    returnVec.push_back(force_MIG);
    returnVec.push_back(force_DIVMIG);
    returnVec.push_back(force_DISEASE);
    returnVec.push_back(force_REC);
    returnVec.push_back(force_GROW);
    returnVec.push_back(force_DIVERGENCE);
    returnVec.push_back(force_LOGISTICSELECTION);
    returnVec.push_back(force_REGION_GAMMA);
    return returnVec;
}

long int UIVarsForces::GetNumGroups(force_type ftype) const
{
    return getLegalForce(ftype).GetNumGroups();
}

long int UIVarsForces::ParamInGroup(force_type ftype, long int pindex) const
{
    return getLegalForce(ftype).ParamInGroup(pindex);
}

bool
UIVarsForces::GetForceCanTurnOnOff(force_type ftype) const
{
    const UIVarsSingleForce & thisForce = getForceRegardlessOfLegality(ftype);
    return thisForce.GetCanSetOnOff();
}

bool
UIVarsForces::GetForceLegal(force_type ftype) const
{
    const UIVarsSingleForce & thisForce = getForceRegardlessOfLegality(ftype);
    return thisForce.GetLegal();
}

bool
UIVarsForces::GetForceZeroesValidity(force_type ftype) const
{
    const UIVarsSingleForce & thisForce = getForceRegardlessOfLegality(ftype);
    return thisForce.AreZeroesValid();
}

void
UIVarsForces::FixGroups(force_type ftype)
{
    UIVarsSingleForce & thisForce = getLegalForce(ftype);
    thisForce.FixGroups();
}

const UIVarsSingleForce &
UIVarsForces::getLegalForce(force_type ftype) const
{
    // if you change this method, make sure you change the
    // non-const version and getForceRegardlessOfLegality too!!
    switch(ftype)
    {
        case force_COAL:
            assert(coalForce.GetLegal());
            return coalForce;
            break;
        case force_MIG:
            assert(migForce.GetLegal());
            return migForce;
            break;
        case force_DIVMIG:
            assert(divMigForce.GetLegal());
            return divMigForce;
            break;
        case force_EXPGROWSTICK:
        case force_GROW:
            assert(growForce.GetLegal());
            return growForce;
            break;
        case force_REC:
            assert(recForce.GetLegal());
            return recForce;
            break;
        case force_DISEASE:
            assert(diseaseForce.GetLegal());
            return diseaseForce;
            break;
        case force_LOGSELECTSTICK:
        case force_LOGISTICSELECTION:
            assert(logisticSelectionForce.GetLegal());
            return logisticSelectionForce;
            break;
        case force_REGION_GAMMA:
            assert(regionGammaForce.GetLegal());
            return regionGammaForce;
            break;
        case force_DIVERGENCE:
            assert(divForce.GetLegal());
            return divForce;
            break;
        default:
            assert(false);              //Uncaught force type.
    }
    throw implementation_error("UIVarsForces::getLegalForce(force_type) given an unknown force type.");
}

UIVarsSingleForce &
UIVarsForces::getLegalForce(force_type ftype)
{
    // if you change this method, make sure you change the
    // const version and getForceRegardlessOfLegality too!!
    switch(ftype)
    {
        case force_COAL:
            assert(coalForce.GetLegal());
            return coalForce;
            break;
        case force_MIG:
            assert(migForce.GetLegal());
            return migForce;
            break;
        case force_DIVMIG:
            assert(divMigForce.GetLegal());
            return divMigForce;
            break;
        case force_EXPGROWSTICK:
        case force_GROW:
            assert(growForce.GetLegal());
            return growForce;
            break;
        case force_REC:
            assert(recForce.GetLegal());
            return recForce;
            break;
        case force_DISEASE:
            assert(diseaseForce.GetLegal());
            return diseaseForce;
            break;
        case force_LOGSELECTSTICK:
        case force_LOGISTICSELECTION:
            assert(logisticSelectionForce.GetLegal());
            return logisticSelectionForce;
            break;
        case force_REGION_GAMMA:
            assert(regionGammaForce.GetLegal());
            return regionGammaForce;
            break;
        case force_DIVERGENCE:
            assert(divForce.GetLegal());
            return divForce;
            break;
        default:
            assert(false);              //Uncaught force type.
    }
    throw implementation_error("UIVarsForces::getLegalForce(force_type) given an unknown force type.");
}

const UIVarsSingleForce &
UIVarsForces::getForceRegardlessOfLegality(force_type ftype) const
{
    // if you change this method, make sure you change the
    // const and non-const versions of getLegalForce
    switch(ftype)
    {
        case force_COAL:
            return coalForce;
            break;
        case force_MIG:
            return migForce;
            break;
        case force_DIVMIG:
            return divMigForce;
            break;
        case force_EXPGROWSTICK:
        case force_GROW:
            return growForce;
            break;
        case force_REC:
            return recForce;
            break;
        case force_DISEASE:
            return diseaseForce;
            break;
        case force_LOGSELECTSTICK:
        case force_LOGISTICSELECTION:
            return logisticSelectionForce;
            break;
        case force_REGION_GAMMA:
            return regionGammaForce;
            break;
        case force_DIVERGENCE:
            return divForce;
            break;
        default:
            assert(false);              //Uncaught force type.
    }
    throw implementation_error("UIVarsForces::getForceRegardlessOfLegality(force_type) given an unknown force type.");
}

long int
UIVarsForces::GetNumParameters(force_type force) const
{
    return getLegalForce(force).GetNumParameters();
}

bool
UIVarsForces::GetForceOnOff   (force_type force) const
{
    const UIVarsSingleForce & thisForce = getForceRegardlessOfLegality(force);
    return thisForce.GetOnOff();
}

long int
UIVarsForces::GetMaxEvents    (force_type force) const
{
    return getLegalForce(force).GetMaxEvents();
}

proftype
UIVarsForces::GetProfileType  (force_type force) const
{
    return getLegalForce(force).GetProfileType();
}

proftype
UIVarsForces::GetProfileType    (force_type force, long int id) const
{
    return getLegalForce(force).GetProfileType(id);
}

string
UIVarsForces::GetProfileTypeSummaryDescription(force_type force)  const
{
    return getLegalForce(force).GetProfileTypeSummaryDescription();
}

ParamStatus
UIVarsForces::GetParamstatus  (force_type force, long int id) const
{
    if (id == uiconst::GLOBAL_ID)
    {
        //LS NOTE:  This is called this way from uiParameter::Min and ::Max in force_interface.cpp
        return ParamStatus(pstat_unconstrained);
    }
    return getLegalForce(force).GetParamstatus(id);
}

ParamStatus
UIVarsForces::GetGroupParamstatus (force_type force, long int id) const
{
    return getLegalForce(force).GetGroupParamstatus(id);
}

LongVec1d
UIVarsForces::GetGroupParamList (force_type force, long int id) const
{
    return getLegalForce(force).GetGroupParamList(id);
}

// MFIX This looked bad in merge; functions GetIdentGroups and GetMultGroups may
// not be functional or correct.
// JDEBUG -- this is a short term kludge to get something/anything in,
// it is wrong WRONG wrong!
vector<ParamGroup>
UIVarsForces::GetIdentGroups (force_type force) const
{
    return getLegalForce(force).GetGroups();
}

// JDEBUG -- this is a short term kludge to get something/anything in,
// it is wrong WRONG wrong!
vector<ParamGroup>
UIVarsForces::GetMultGroups (force_type force) const
{
    return getLegalForce(force).GetGroups();
}

double
UIVarsForces::GetStartValue   (force_type force, long int id) const
{
    return getLegalForce(force).GetStartValue(id);
}

double
UIVarsForces::GetTrueValue   (force_type force, long int id) const
{
    return getLegalForce(force).GetTrueValue(id);
}

force_type
UIVarsForces::GetPhase2Type   (force_type force) const
{
    return getLegalForce(force).GetPhase2Type(force);
}

method_type
UIVarsForces::GetStartMethod  (force_type force, long int id) const
{
    return getLegalForce(force).GetStartMethod(id);
}

DoubleVec1d
UIVarsForces::GetStartValues (force_type force) const
{
    return getLegalForce(force).GetStartValues();
}

long int
UIVarsForces::GetDiseaseLocation() const
{
    assert(diseaseForce.GetLegal());
    return diseaseForce.GetLocation();
}

growth_type
UIVarsForces::GetGrowthType() const
{
    return growForce.GetGrowthType();
}

growth_scheme
UIVarsForces::GetGrowthScheme() const
{
    return growForce.GetGrowthScheme();
}

bool
UIVarsForces::GetDoProfile(force_type force, long int id) const
{
    return getLegalForce(force).GetDoProfile(id);
}

bool
UIVarsForces::GetParamValid(force_type force, long int id) const
{
    return getLegalForce(force).GetParamValid(id);
}

bool
UIVarsForces::GetParamUnique(force_type force, long int id) const
{
    return getLegalForce(force).GetParamUnique(id);
}

//Getters/setters for Bayesian information
const UIVarsPrior& UIVarsForces::GetPrior(force_type force, long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return getLegalForce(force).GetDefaultPrior();
    }
    return getLegalForce(force).GetPrior(pindex);
}

const UIVarsPrior& UIVarsForces::GetDefaultPrior(force_type force) const
{
    return getLegalForce(force).GetDefaultPrior();
}

bool UIVarsForces::GetUseDefaultPrior(force_type force, long int pindex) const
{
    return getLegalForce(force).GetUseDefaultPrior(pindex);
}

priortype UIVarsForces::GetPriorType(force_type force, long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return getLegalForce(force).GetDefaultPriorType();
    }
    return getLegalForce(force).GetPriorType(pindex);
}

double    UIVarsForces::GetLowerBound(force_type force, long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return getLegalForce(force).GetDefaultLowerBound();
    }
    return getLegalForce(force).GetLowerBound(pindex);
}

selection_type
UIVarsForces::GetSelectionType() const
{
    return logisticSelectionForce.GetSelectionType();
}

double    UIVarsForces::GetUpperBound(force_type force, long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return getLegalForce(force).GetDefaultUpperBound();
    }
    return getLegalForce(force).GetUpperBound(pindex);
}

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
long    UIVarsForces::GetRelativeSampling(force_type force, long int pindex) const
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        return getLegalForce(force).GetDefaultRelativeSampling();
    }
    return getLegalForce(force).GetRelativeSampling(pindex);
}
#endif

void UIVarsForces::SetUseDefaultPrior(bool use, force_type force, long int pindex)
{
    getLegalForce(force).SetUseDefaultPrior(use, pindex);
}

void UIVarsForces::SetUseDefaultPriorsForForce(force_type force)
{
    getLegalForce(force).SetUseAllDefaultPriors();
}

void UIVarsForces::SetPriorType(priortype ptype, force_type force, long int pindex)
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        getLegalForce(force).SetDefaultPriorType(ptype);
    }
    else getLegalForce(force).SetPriorType(ptype, pindex);
}

#ifdef LAMARC_NEW_FEATURE_RELATIVE_SAMPLING
void UIVarsForces::SetRelativeSampling(long rate, force_type force, long int pindex)
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        getLegalForce(force).SetDefaultRelativeSampling(rate);
    }
    else getLegalForce(force).SetRelativeSampling(rate, pindex);
}
#endif


void UIVarsForces::SetLowerBound(double bound, force_type force, long int pindex)
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        getLegalForce(force).SetDefaultLowerBound(bound);
    }
    else getLegalForce(force).SetLowerBound(bound, pindex);
}

void UIVarsForces::SetUpperBound(double bound, force_type force, long int pindex)
{
    if (pindex == uiconst::GLOBAL_ID)
    {
        getLegalForce(force).SetDefaultUpperBound(bound);
    }
    else getLegalForce(force).SetUpperBound(bound, pindex);
}

string UIVarsForces::GetPriorTypeSummaryDescription(force_type force) const
{
    return getLegalForce(force).GetPriorTypeSummaryDescription();
}

string UIVarsForces::GetPriorTypeSummaryDescription(force_type force, long int pindex, bool sayDefault) const
{
    if (pindex != uiconst::GLOBAL_ID)
    {
        ParamStatus mystatus = getLegalForce(force).GetParamstatus(pindex);
        if (!mystatus.Inferred()) return "<" + ToString(mystatus.Status()) + ">";
        // otherwise do the stuff below
    }
    string desc = getLegalForce(force).GetPrior(pindex).GetSummaryDescription();
    if (pindex != uiconst::GLOBAL_ID)
    {
        if (sayDefault)
        {
            if (getLegalForce(force).GetUseDefaultPrior(pindex))
            {
                desc = "<default>";
            }
        }
    }
    return desc;
}

void
UIVarsForces::SetForceOnOff  (bool onOffVal, force_type force)
{
    if(GetForceCanTurnOnOff(force))
    {
        return getLegalForce(force).SetOnOff(onOffVal);
    }
    else
    {
        if (onOffVal != GetForceOnOff(force))
        {
            assert(false);
            string msg;
            msg = "Unable to turn ";
            if (onOffVal)
            {
                msg += "on";
            }
            else
            {
                msg += "off";
            }
            msg += " the force " + ToString(force)
                + " for this dataset.  Try removing the force entirely from the "
                "XML input file.";
            throw data_error(msg);
        }
    }
}

void
UIVarsForces::SetMaxEvents   (long int maxEvents,    force_type force)
{
    return getLegalForce(force).SetMaxEvents(maxEvents);
}

void
UIVarsForces::SetProfileType (proftype ptype,    force_type force)
{
    return getLegalForce(force).SetProfileType(ptype);
}

void
UIVarsForces::SetProfileType (proftype ptype)
{
    // OK, so it's pretty weird to do this for all of the
    // forces, since some of them will not be legal, but
    // hey, if we're never accessing them, we can certainly
    // write to them, huh?
    coalForce.SetProfileType(ptype);
    diseaseForce.SetProfileType(ptype);
    growForce.SetProfileType(ptype);
    migForce.SetProfileType(ptype);
    recForce.SetProfileType(ptype);
    logisticSelectionForce.SetProfileType(ptype);
    regionGammaForce.SetProfileType(ptype);
}

void
UIVarsForces::SetDoProfile   (bool doProfile,    force_type force, long int id)
{
    getLegalForce(force).SetDoProfile(doProfile,id);
}

void
UIVarsForces::SetDoProfile   (bool doProfile,    force_type force)
{
    getLegalForce(force).SetDoProfile(doProfile);
}

void
UIVarsForces::SetDoProfile   (bool doProfile)
{
    // OK, so it's pretty wierd to do this for all of the
    // forces, since some of them will not be legal, but
    // hey, if we're never accessing them, we can certainly
    // write to them, huh?
    coalForce.SetDoProfile(doProfile);
    diseaseForce.SetDoProfile(doProfile);
    growForce.SetDoProfile(doProfile);
    migForce.SetDoProfile(doProfile);
    recForce.SetDoProfile(doProfile);
    logisticSelectionForce.SetDoProfile(doProfile);
    regionGammaForce.SetDoProfile(doProfile);
}

void
UIVarsForces::SetStartMethod(method_type meth, force_type ftype, long int id)
{
    getLegalForce(ftype).SetStartMethod(meth,id);
}

void
UIVarsForces::SetParamstatus (const ParamStatus& mystatus, force_type ftype, long int id)
{
    getLegalForce(ftype).SetParamstatus(mystatus,id);
}

void
UIVarsForces::SetGroupParamstatus (ParamStatus pstat, force_type ftype, long int id)
{
    getLegalForce(ftype).SetGroupParamstatus(pstat,id);
}

void
UIVarsForces::AddGroup (LongVec1d params, force_type ftype, long int id)
{
    getLegalForce(ftype).AddGroup(defaults::groupPstat, params);
}

void
UIVarsForces::RemoveParamFromGroup(force_type ftype, long int id)
{
    getLegalForce(ftype).RemoveParamIfInAGroup(id);
}

void
UIVarsForces::AddParamToGroup(force_type ftype, long int pindex, long int gindex)
{
    getLegalForce(ftype).AddParamToGroup(pindex, gindex);
}

void
UIVarsForces::AddParamToNewGroup(force_type ftype, long int pindex)
{
    getLegalForce(ftype).AddParamToNewGroup(pindex);
}

void
UIVarsForces::SetUserStartValue(double startValue, force_type ftype, long int id)
{
    getLegalForce(ftype).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetTrueValue(double startValue, force_type ftype, long int id)
{
    getLegalForce(ftype).SetTrueValue(startValue,id);
}

void
UIVarsForces::SetAllThetaStartValues(double startValue)
{
    getLegalForce(force_COAL).SetUserStartValues(startValue);
}

void
UIVarsForces::SetAllThetaStartValuesFST()
{
    getLegalForce(force_COAL).SetStartMethods(method_FST);
}

void
UIVarsForces::SetAllThetaStartValuesWatterson()
{
    getLegalForce(force_COAL).SetStartMethods(method_WATTERSON);
}

void
UIVarsForces::SetThetaStartValue  (double startValue, long int id)
{
    getLegalForce(force_COAL).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetAllMigrationStartValues(double startValue)
{
    getLegalForce(force_MIG).SetUserStartValues(startValue);
}

void
UIVarsForces::SetAllMigrationStartValuesFST()
{
    getLegalForce(force_MIG).SetStartMethods(method_FST);
}

void
UIVarsForces::SetMigrationStartValue  (double startValue, long int id)
{
    getLegalForce(force_MIG).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetDivergenceEpochStartTime  (double startValue, long int id)
{
    getLegalForce(force_DIVERGENCE).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetAllDivMigrationStartValues(double startValue)
{
    getLegalForce(force_DIVMIG).SetUserStartValues(startValue);
}

#if 0 // JREMOVE, JDEBUG, passes testing then remove this function
void
UIVarsForces::SetAllDivMigrationStartValuesFST()
{
    assert(false);
    getLegalForce(force_DIVMIG).SetStartMethods(method_FST);
}
#endif

void
UIVarsForces::SetDivMigrationStartValue  (double startValue, long int id)
{
    getLegalForce(force_DIVMIG).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetAllDiseaseStartValues(double startValue)
{
    getLegalForce(force_DISEASE).SetUserStartValues(startValue);
}

void
UIVarsForces::SetDiseaseStartValue(double startValue, long int id)
{
    getLegalForce(force_DISEASE).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetAllGrowthStartValues(double startValue)
{
    getLegalForce(force_GROW).SetUserStartValues(startValue);
}

void
UIVarsForces::SetGrowthStartValue(double startValue, long int id)
{
    getLegalForce(force_GROW).SetUserStartValue(startValue,id);
}

void
UIVarsForces::SetGrowthType(growth_type gType)
{
    growForce.SetGrowthType(gType);
}

void
UIVarsForces::SetGrowthScheme(growth_scheme gScheme)
{
    growForce.SetGrowthScheme(gScheme);
}

void
UIVarsForces::SetRecombinationStartValue(double startValue)
{
    getLegalForce(force_REC).SetUserStartValues(startValue);
}

void
UIVarsForces::SetLogisticSelectionCoefficientStartValue(double startValue)
{
    getLegalForce(force_LOGISTICSELECTION).SetUserStartValues(startValue);
}

void
UIVarsForces::SetRegionGammaStartValue(double startValue)
{
    getLegalForce(force_REGION_GAMMA).SetUserStartValues(startValue);
}

void
UIVarsForces::SetSelectionType(selection_type sType)
{
    logisticSelectionForce.SetSelectionType(sType);
}

void
UIVarsForces::SetDiseaseLocation(long int location)
{
    assert(diseaseForce.GetLegal());
    diseaseForce.SetLocation(location);
}

bool UIVarsForces::AreZeroesValid(force_type forceId)
{
    return getLegalForce(forceId).AreZeroesValid();
}

bool UIVarsForces::SomeVariableParams() const
{
    vector<force_type> forces = GetActiveForces();
    for (unsigned long int fnum = 0; fnum < forces.size(); fnum++)
    {
        if (getLegalForce(forces[fnum]).SomeVariableParams())
        {
            return true;
        }
    }
    return false;
}

std::string UIVarsForces::GetEpochAncestorName(long id) const
{
    assert(divForce.GetLegal());
    std::vector<string>  names = divForce.GetAncestors();
    //std::vector<string> mine = names[id];
    std::string retval = names[id];
    return retval;
}

std::string UIVarsForces::GetEpochDescendentNames(long id) const
{
    assert(divForce.GetLegal());
    std::vector< std::vector <string> > names = divForce.GetNewPops();
    std::vector<string> mine = names[id];
    std::string retval = mine[0] + " " + mine[1];
    return retval;
}

//____________________________________________________________________________________