xref: /dports/biology/lamarc/lamarc-2.1.8/src/postlike/plforces.cpp

// $Id: plforces.cpp,v 1.74 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.
*/

// Posterior Likelihood forces classes
// - CoalescePL
//   - coalesce with no-growth [using compressed summaries]
// - CoalesceGrowPL
//   - coalesce with exponential growth [using non-compressed summaries]
// - GrowPL
//   - coalesce with exponential growth [using non-compressed summaries]
// - MigratePL
// - RecombinePL
// - SelectPL [stubs]

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

#include <cassert>

#include "local_build.h"

#include "mathx.h"
#include "plforces.h"
#include "runreport.h"
#include "summary.h"
#include "vectorx.h"
#include "force.h"
#include "timemanager.h"                // for CoalesceGrowPL::lnPoint()

#ifdef DMALLOC_FUNC_CHECK
#include <dmalloc.h>
#endif

// DEBUG CODE HELPERS
#include <iostream>
#include <fstream>
#include <sstream>
std::ofstream numfile;

using namespace std;

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

const double LOG_ONEPLUSEPSILON = log(1.0 + numeric_limits<double>::epsilon());

//------------------------------------------------------------------------------------
// DiseasePL: disease force specific waiting time and
//            point probabilites plus derivatives
//------------------------------------------------------------------------------------

double DiseasePL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    vector < double >::iterator i;
    vector < double >::const_iterator pstart;
    vector < double >::const_iterator pend;

    const vector<double>& mWait = treedata->GetDiseaseSummary()->GetShortWait();
    // precomputing the transition rates to a status Sum[d_ji] into
    // membervariable m_msum
    for (i = m_msum.begin (),
             pstart = param.begin () + m_start,
             pend = param.begin () + m_start + m_nPop;
         i != m_msum.end (); ++i, pstart += m_nPop, pend += m_nPop)
    {
        (*i) = accumulate (pstart, pend, 0.0);
    }

    // sum_pop(SM_status *kt)
    return  -1.0 * inner_product (mWait.begin (),
                                  mWait.end (), m_msum.begin (), 0.0);
}

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

double DiseasePL::lnPoint (const vector < double >&param, const vector < double >&lparam,
                           const TreeSummary * treedata)
{
    //  Sum_j(Sum_i(disevent[j,i] * log(M[j,i]))
    const vector<double>& nevents = treedata->GetDiseaseSummary()->GetShortPoint();
    return  inner_product (nevents.begin (), nevents.end (),
                           lparam.begin () + m_start, 0.0);
}

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

double DiseasePL::DlnWait (const vector < double >&param, const TreeSummary * treedata, const long int &whichparam)
{
    long int which = (whichparam - m_start) / m_nPop;
    const vector<double>& mWait = treedata->GetDiseaseSummary()->GetShortWait();
    return -mWait[which];
}

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

double DiseasePL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long int &whichparam)
{
    long int which = whichparam - m_start;
    const vector<double>& nmig = treedata->GetDiseaseSummary()->GetShortPoint();
    return SafeDivide (nmig[which], param[whichparam]);
}

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

double DiseaseLogisticSelectionPL::lnWait(const vector<double>& param, const TreeSummary * treedata)
{
    return 0.0; // contained in CoalesceLogisticSelectionPL::lnWait(), for speed
}

//------------------------------------------------------------------------------------
// Note:  DlnPoint() is selection-independent, because the point terms
// separate into a sum of logarithms.  So, the DiseasePL method is used for that.

double DiseaseLogisticSelectionPL::lnPoint(const vector<double>& param, const vector<double>& lparam,
                                           const TreeSummary * treedata)
{
    const Interval *treesum = treedata->GetDiseaseSummary()->GetLongPoint();
    const Interval *pTreesum;

    double logTheta_A0(lparam[0]), logTheta_a0(lparam[1]), s(param[m_s_is_here]);
    double log_mu_into_A_from_a(lparam[3]), log_mu_into_a_from_A(lparam[4]);

    if (0 != param[2] || 0 != param[5] || 0 == param[3] || 0 == param[4])
        throw implementation_error("Error parsing disease rates in DiseaseLogisticSelectionPL::lnPoint()");

    double result = 0.0;

    for(pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
    {
        if (0L == pTreesum->oldstatus)
            result += log_mu_into_A_from_a + logTheta_a0 - logTheta_A0 + s*pTreesum->endtime;
        else if (1L == pTreesum->oldstatus)
            result += log_mu_into_a_from_A + logTheta_A0 - logTheta_a0 - s*pTreesum->endtime;
        else
        {
            string msg = "DiseaseLogisticSelectionPL::lnPoint(), received a TreeSummary ";
            msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
            msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescence in either ";
            msg += "the subpopulation with allele A or the subpopulation with allele a.";
            throw implementation_error(msg);
        }
    }

    return result;
}

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

double DiseaseLogisticSelectionPL::DlnWait(const vector<double>&param, const TreeSummary * treedata,
                                           const long int &whichparam)
{
    if (2 != m_nPop)
    {
        string msg = "DiseaseLogisticSelectionPL::DlnWait() called with m_nPop = ";
        msg += ToString(m_nPop);
        msg += ".  m_nPop must equal 2, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }

    if (3 != whichparam && 4 != whichparam)
        throw implementation_error("DiseaseLogisticSelectionPL::DlnWait(), bad \"whichparam\"");

    const list<Interval>& treesum = treedata->GetDiseaseSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;

    double result(0.0), s(param[m_s_is_here]);
    double theta_A0(param[m_start]), theta_a0(param[m_start+1]);

    if (theta_a0 <= 0.0 || theta_A0 <= 0.0)
    {
        string msg = "DiseaseLogisticSelectionPL::DlnWait(), received an invalid Theta value ";
        msg += "(theta_A0 = " + ToString(theta_A0) + ", theta_a0 = " + ToString(theta_a0) + ").";
        throw impossible_error(msg);
    }

    if (fabs(s) < LOGISTIC_SELECTION_COEFFICIENT_EPSILON)
        return DlnWaitForTinyLogisticSelectionCoefficient(param, treedata, whichparam);

    if (fabs(s) >= DBL_BIG)
        return -DBL_BIG; // unavoidable overflow

    double t_e, t_s(0.0), dt, e_toThe_sts(1.0), e_toThe_ste, term(0.0);
    double term_A(0.0), term_a(0.0);
    double factor_A(theta_a0/(s*theta_A0)), factor_a(theta_A0/(-s*theta_a0));
    bool derivativeWithRespectTo_mu_into_A_from_a = (3 == whichparam ? true : false);

    for (treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        dt = t_e - t_s;
        term_A = term_a = 0.0;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        if (s > 0.0)
        {
            // Use overflow protection for term_A; if term_a underflows, that's okay.
            if (s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor_A,s*dt);
                if (term >= DBL_BIG && derivativeWithRespectTo_mu_into_A_from_a)
                    return -DBL_BIG; // unavoidable overflow
                e_toThe_ste = (term/factor_A)*e_toThe_sts;
                if (k_A > 0.0 && derivativeWithRespectTo_mu_into_A_from_a)
                {
                    term_A = k_A*(term - factor_A)*e_toThe_sts;
                    if (term_A >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_A > 0.0 && derivativeWithRespectTo_mu_into_A_from_a)
                {
                    term_A = k_A*factor_A*s*dt*e_toThe_sts;
                    if (term_A >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            if (k_a > 0.0 && !derivativeWithRespectTo_mu_into_A_from_a)
                term_a = k_a*factor_a*(1.0/e_toThe_ste - 1.0/e_toThe_sts); // note: factor_a includes minus sign
        }
        else // s < 0
        {
            // Use overflow protection for term_a; if term_A underflows, that's okay.

            if (-s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor_a,-s*dt);
                if (term >= DBL_BIG && !derivativeWithRespectTo_mu_into_A_from_a)
                    return -DBL_BIG; // unavoidable overflow
                e_toThe_ste = e_toThe_sts/(term/factor_a);
                if (k_a > 0.0 && !derivativeWithRespectTo_mu_into_A_from_a)
                {
                    term_a = k_a*(term - factor_a)/e_toThe_sts;
                    if (term_a >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (-s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_a > 0.0 && !derivativeWithRespectTo_mu_into_A_from_a)
                {
                    term_a = k_a*factor_a*(-s)*dt/e_toThe_sts;
                    if (term_a >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            if (k_A > 0.0 && derivativeWithRespectTo_mu_into_A_from_a)
                term_A = k_A*factor_A*(e_toThe_ste - e_toThe_sts);
        }

        if (derivativeWithRespectTo_mu_into_A_from_a)
            result -= term_A;
        else
            result -= term_a;
        if (result <= -DBL_BIG)
            return -DBL_BIG;

        t_s = t_e;
        e_toThe_sts = e_toThe_ste;
    }

    return result;
}

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

double DiseaseLogisticSelectionPL::DlnWaitForTinyLogisticSelectionCoefficient(const vector<double>& param,
                                                                              const TreeSummary *treedata,
                                                                              const long int &whichparam)
{
    const list<Interval>& treesum = treedata->GetDiseaseSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;
    double result(0.0), s(param[m_s_is_here]);
    double theta_A0(param[m_start]), theta_a0(param[m_start+1]);
    double t_e, t_s(0.0);
    bool derivativeWithRespectTo_mu_into_A_from_a = (4 == whichparam ? true : false);

    for (treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        if (derivativeWithRespectTo_mu_into_A_from_a)
            result -= (k_A*theta_a0/theta_A0)*((t_e - t_s) + 0.5*s*(t_e*t_e - t_s*t_s));
        else
            result -= (k_a*theta_A0/theta_a0)*((t_e - t_s) - 0.5*s*(t_e*t_e - t_s*t_s));

        t_s = t_e;
    }

    return result;
}

//------------------------------------------------------------------------------------
// MigratePL: migration forces specific waiting time and
//            point probabilites and its derivatives

double MigratePL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    vector < double >::iterator i;
    vector < double >::const_iterator pstart;
    vector < double >::const_iterator pend;

    const vector<double>& mWait = treedata->GetMigSummary()->GetShortWait();

    // precomputing the immigration rates Sum[m_ji] into membervariable m_msum
    for (i = m_msum.begin (),
             pstart = param.begin () + m_start,
             pend = param.begin () + m_start + m_nPop;
         i != m_msum.end (); ++i, pstart += m_nPop, pend += m_nPop)
    {
        (*i) = accumulate (pstart, pend, 0.0);
    }

    // sum_pop(SM_pop *kt)
    return  -1.0 * inner_product (mWait.begin (),
                                  mWait.end (), m_msum.begin (), 0.0);
}

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

double MigratePL::lnPoint (const vector < double >&param, const vector < double >&lparam,
                           const TreeSummary * treedata)
{
    //  Sum_j(Sum_i(migevent[j,i] * log(M[j,i]))
    const vector<double>& nmig = treedata->GetMigSummary()->GetShortPoint();

    return  inner_product (nmig.begin(), nmig.end(), lparam.begin() + m_start,
                           0.0);
}

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

double MigratePL::DlnWait (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    long int which = (whichparam - m_start) / m_nPop;
    const vector<double>& mWait = treedata->GetMigSummary()->GetShortWait();
    return -mWait[which];
}

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

double MigratePL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    long int which = whichparam - m_start;
    const vector<double>& nmig = treedata->GetMigSummary()->GetShortPoint();

    return SafeDivide (nmig[which], param[whichparam]);
}

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

// CoalescePL: coalescence forces specific waiting time and
//             point probabilites and its derivatives

// CalculateScaledLPCounts returns the log-likelihood of all
// of the choices of recombination-branch partition assignment
// made in a genealogy.  It does this by summing the log of (the
// number of times each partition was chosen times the relative
// frequency of that partition).
//
// JDEBUG--:  This code cannot handle migration.  To do so, it
// will be necessary to (a) store the partition in which each
// recombination happened, and (b) make the relative frequency
// relative to the partition, not the whole.  For example, the
// term for choice of a disease state in Boston should be
// Theta(disease & Boston) over Theta(Boston).  Currently it
// is Theta(disease) over Theta(total).

double CoalescePL::CalculateScaledLPCounts(const DoubleVec1d& params, const LongVec2d& picks) const
{
    double total_theta(std::accumulate(params.begin()+m_start,
                                       params.begin()+m_end,0.0));
    DoubleVec1d myths(params.begin()+m_start,params.begin()+m_end);
    double answ(0.0);
    LongVec2d::size_type lpforce;
    for(lpforce = 0; lpforce < picks.size(); ++lpforce)
    {
        LocalPartitionForce* pforce(dynamic_cast<LocalPartitionForce*>
                                    (m_localpartforces[lpforce]));
        assert(pforce);  // it wasn't a local partition force?
        DoubleVec1d part_thetas(pforce->SumXPartsToParts(myths));
        LongVec1d::size_type partition;
        for(partition = 0; partition < picks[lpforce].size(); ++partition)
        {
            answ += log(picks[lpforce][partition] * part_thetas[partition] / total_theta);
        }
    }
    return answ;
}

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

void CoalescePL::SetLocalPartForces(const ForceSummary& fs)
{
    m_localpartforces = fs.GetLocalPartitionForces();
    if (m_localpartforces.empty()) return;

    // set up vector which answers the following question:
    // for each local partition force
    //   for any xpartition theta (identified by index),
    //   what partition (of the given force) is it?

    // This code is copied approximately from PartitionForce::
    // SumXPartsToParts().
    ForceVec::const_iterator pforce;
    const ForceVec& partforces = fs.GetPartitionForces();
    LongVec1d indicator(partforces.size(), 0L);
    LongVec1d nparts(registry.GetDataPack().GetAllNPartitions());
    DoubleVec1d::size_type xpart;
    long partindex;
    for (pforce = partforces.begin(), partindex = 0;
         pforce != partforces.end(); ++pforce, ++partindex)
    {
        vector<DoubleVec1d::size_type> indices(m_nPop, 0);
        for(xpart = 0; xpart < static_cast<DoubleVec1d::size_type>(m_nPop); ++xpart)
        {
            indices[xpart] = indicator[partindex];
            long int part;
            for (part = nparts.size() - 1; part >= 0; --part)
            {
                ++indicator[part];
                if (indicator[part] < nparts[part]) break;
                indicator[part] = 0;
            }
        }
        // initialize xparts vectors
        vector<DoubleVec1d::size_type> emptyvec;
        vector<vector<DoubleVec1d::size_type> > partvec(nparts[partindex], emptyvec);

        if ((*pforce)->IsLocalPartitionForce())
        {
            m_whichlocalpart.push_back(indices);
            m_whichlocalxparts.push_back(partvec);
        }
        m_whichpart.push_back(indices);
        m_whichxparts.push_back(partvec);
    }

    // now construct the vector mapping partition to a set of xpartitions.
    //
    // we will do this by going through the vector mapping xpartition to
    // partition we just constructed, letting its contents tell us which
    // partition to add the parameter (xpartition) to and keeping a
    // counter to let us know which xpartition to add.
    long lpindex = 0;
    for (pforce = partforces.begin(), partindex = 0;
         pforce != partforces.end(); ++pforce, ++partindex)
    {
        for(xpart = 0; xpart < m_whichpart[partindex].size(); ++xpart)
            m_whichxparts[partindex][m_whichpart[partindex][xpart]].push_back(xpart);

        if ((*pforce)->IsLocalPartitionForce())
        {
            for(xpart = 0; xpart < m_whichpart[partindex].size(); ++xpart)
            {
                m_whichlocalxparts[lpindex][m_whichpart[partindex][xpart]].
                    push_back(xpart);
            }
            ++lpindex;
        }
    }
} // SetLocalPartForces

//------------------------------------------------------------------------------------
// COMPRESSED DATA (SHORT) SUMMARIES

double CoalescePL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    // result = sum_pop(k(k-1)t/theta_pop)
    const vector<double>& cWait = treedata->GetCoalSummary()->GetShortWait();

    return -1.0 * inner_product (cWait.begin (), cWait.end (),
                                 param.begin () + m_start,
                                 0.0, plus < double >(),
                                 divides < double >());
}

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

double CoalescePL::lnPoint (const vector < double >&param, const vector < double >&lparam,
                            const TreeSummary * treedata)
{
    //  in general it is Sum_j(coalesceevent[j] * (log(2) - log(theta_j)))
    //  sumcpoint = Sum_j(coalesceevent[j] * log(2))
    //   - Sum_j(coalesceevent[j] * log(theta_j))
    const vector<double>& ncoal = treedata->GetCoalSummary()->GetShortPoint();

    double point = (LOG2 * accumulate (ncoal.begin (), ncoal.end (), 0.0) -
                    inner_product (ncoal.begin (), ncoal.end (),
                                   lparam.begin () + m_start, 0.0));
    const LongVec2d& picks = treedata->GetCoalSummary()->GetShortPicks();
    if (picks.empty()) return point;
    else return point + CalculateScaledLPCounts(param, picks);
}

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

double CoalescePL::DlnWait (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    long int which = whichparam - m_start;
    const vector<double>& cWait = treedata->GetCoalSummary()->GetShortWait();

    return cWait[which] / (param[whichparam] * param[whichparam]);
}

//------------------------------------------------------------------------------------
// erynes 2004/03/11 -- This method presently returns the following:
//
//          -ncoal/theta,
//
// where theta is the current estimate of parameter theta for the
// given population ("whichparam"), and ncoal is the total number
// of coalescent events for this population (i.e., one less than the
// total number of tips that the population has in this genealogy tree,
// unless recombinations are present).
//
// Note:  CoalescePL::DlnPoint calculates the partial derivative of log(Point)
// with respect to theta for the given population ("whichparam"),
// in both the absence and presence of growth.
// When growth is present, CoalesceGrowPL::DlnWait does this for log(Wait).
// Partial derivatives with respect to growth are found in GrowPL.
//
// jay 2007/03/22 -- added support for disease w/ recombination
//
// JDEBUG--:  only in absence of migration

double CoalescePL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    long int which = whichparam - m_start;
    const vector<double>& ncoal = treedata->GetCoalSummary()->GetShortPoint();

    double answ(-ncoal[which] / param[whichparam]);

    const LongVec2d& picks = treedata->GetCoalSummary()->GetShortPicks();
    if (picks.empty()) return answ;
    else  {
        // count recs in crosspartition corresponding to whichparam
        const RecSummary* rsum = dynamic_cast<const RecSummary*>(treedata->GetSummary(force_REC));
        long recs_in_xpart = rsum->GetNRecsByXPart()[whichparam];
        // divide by whichparam
        double xpart_theta = param[whichparam];
        answ += recs_in_xpart / xpart_theta;
        // from that, subtract number of recs (should be "in that partition")
        long nrecs = rsum->GetNRecs();
        // divided by sum of thetas (should be "for that partition")
        double total_theta(std::accumulate(param.begin()+m_start,
                                           param.begin()+m_end,0.0));
        answ -= nrecs/total_theta;
    }

#if 0 // this was what CVS merged into here.....
    if (!m_localpartforces.empty())
    {
        double total_theta(std::accumulate(param.begin()+m_start,
                                           param.begin()+m_end,0.0));
        DoubleVec1d myths(param.begin()+m_start,param.begin()+m_end);
        long int total_npartitions(0);
        LongVec2d::size_type lpforce;
        for(lpforce = 0; lpforce < m_localpartforces.size(); ++lpforce)
        {
            LocalPartitionForce* pforce(dynamic_cast<LocalPartitionForce*>
                                        (m_localpartforces[lpforce]));
            assert(pforce);  // it wasn't a local partition force?
            total_npartitions += pforce->GetNPartitions();

            const vector<DoubleVec1d::size_type>& indices(
                m_whichxparts[lpforce][m_whichpart[lpforce][which]]);
            vector<DoubleVec1d::size_type>::const_iterator index;
            double part_theta(0.0);
            for(index = indices.begin(); index != indices.end(); ++index)
            {
                part_theta += myths[*index];
            }
            answ += 1.0 / part_theta;

        }
        answ -= total_npartitions / total_theta;
    }
#endif

    return answ;
}

//------------------------------------------------------------------------------------
// CoalesceGrowPL: coalescence forces specific waiting time and
//             point probabilites and its derivatives
// LONG DATA SUMMARIES
// [only in effect when GROWTH force is turned on]

//------------------------------------------------------------------------------------
// erynes 2004/03/11 -- This method presently returns the following:
//
// sum_over_populations(sum_over_each_population's_coalescent_time_intervals(
//  (k_i*(k_i - 1)/(theta_p * g_p)) * (exp(g_p * t_start)- exp(g_p * t_end))
//                                                                          )
//                     )
//
// where
//        k_i is the number of lineages of population p in time interval i
//        theta_p is the present estimate of parameter theta for population p
//        g_p is the present estimate of the growth parameter for population p
//        t_start is the timepoint at the start of time interval i
//        t_end is the timepoint at the end of time interval i
//
// This quantity is equal to log(WaitProb(G|P)), where WaitProb(G|P) is the
// "waiting" probability of the input genealogy G (parameter "treedata")
// given the parameters P (parameter "param").
// Prob(G|P) = PointProb(G|P) * WaitProb(G|P).
// See CoalesceGrowPL::Point() for the point probability term.

double CoalesceGrowPL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    const list<Interval>& treesum=treedata->GetCoalSummary()->GetLongWait();
    list <Interval> :: const_iterator tit;
    vector < double > :: const_iterator pit = param.begin() + m_start;
    vector < double > :: const_iterator pend = param.begin() + m_start + m_nPop;
    vector < double > :: const_iterator git = param.begin() + m_growthstart;

    unsigned long int xpartition = 0L;
    double growth, theta, gts, gte;
    double coeff_s, coeff_e, arg1_s, arg1_e;
    double result = 0.0;
    double starttime;

    for( ; pit != pend; ++pit, ++git, ++xpartition)
    {
        growth = *git;
        theta = *pit;

        if (fabs(growth) < GROWTH_EPSILON)
        {
            result += lnWaitForTinyGrowth(theta, growth, xpartition, treesum);
            continue;
        }

        starttime = 0.0;
        for(tit = treesum.begin(); tit != treesum.end(); ++tit)
        {
            gts = growth * starttime;
            gte = growth * (*tit).endtime;

            coeff_s = coeff_e = arg1_s = arg1_e = 1.0;
            const double& k = (*tit).xpartlines[xpartition];

            // IMPORTANT NOTE:  The following over/underflow checking
            // assumes that each starttime and endtime is always >= 0.
            // Hence, e.g., gts < 0 means growth < 0.

            if (gts > 1.0)
            {
                coeff_s = k * (k - 1.0);
                if (growth > 1.0)
                    arg1_s /= growth;
                else
                    coeff_s /= growth;
                if (theta > 1.0)
                    arg1_s /= theta;
                else
                    coeff_s /= theta;
            }
            else if (gts < -1.0)
            {
                arg1_s *= k * (k - 1.0);
                if (growth < -1.0)
                    coeff_s /= growth;
                else
                    arg1_s /= growth;
                if (theta > 1.0)
                    coeff_s /= theta;
                else
                    arg1_s /= theta;
            }
            else // no over/underflow in exp()
                coeff_s = k * (k - 1.0)/(theta * growth); // arg1_s already set

            // BUGBUG erynes -- what follows is easy to understand,
            // but we should revisit this and reuse some of what we
            // learned about growth and theta above, while considering gts.

            if (gte > 1.0)
            {
                coeff_e = k * (k - 1.0);
                if (growth > 1.0)
                    arg1_e /= growth;
                else
                    coeff_e /= growth;
                if (theta > 1.0)
                    arg1_e /= theta;
                else
                    coeff_e /= theta;
            }
            else if (gte < -1.0)
            {
                arg1_e *= k * (k - 1.0);
                if (growth < -1.0)
                    coeff_e /= growth;
                else
                    arg1_e /= growth;
                if (theta > 1.0)
                    coeff_e /= theta;
                else
                    arg1_e /= theta;
            }
            else // no over/underflow in exp()
                coeff_e = k * (k - 1.0)/(theta * growth); // arg1_e already set

            double incrementalResult =
                coeff_s * SafeProductWithExp(arg1_s, gts) -
                coeff_e * SafeProductWithExp(arg1_e, gte);

            if (incrementalResult > 0.0)
            {
                // We have precision problems,
                // because we expect te > ts always,
                // and hence gte > gts always.
                // During debugging, we have encountered precision problems
                // in which gdb claims gts == gte and arg1_s == arg1_e
                // but incrementalResult > 0, typically reflecting
                // a difference near the 15th digit.
                // (Note that te == ts corresponds to two events
                // occurring simultaneously--e.g., k = 1 for a really long
                // time in each of two populations, then one finally
                // migrates to the other, then they coalesce instantaneously.)
                // It should be safe to set incrementalResult to zero in cases of
                // roundoff error.  If incrementalResult is positive,
                // but _not_ tiny compared with either of the terms
                // whose difference yielded incrementalResult,
                // then something is _very_ wrong--probably te < ts,
                // which should never happen.  We assert in this case;
                // this will enable us to see the problem, but the end user
                // will not see the assertion, unless s/he's running in debug mode.
                if (incrementalResult/fabs(coeff_s * SafeProductWithExp(arg1_s, gts)) > 1e-10)
                {
                    string msg = "incrementalResult in plforces.cpp is greater "
                        "than zero.  coeff_s = "
                        + ToString(coeff_s) + ", coeff_e = " + ToString(coeff_e)
                        + ", theta = " + ToString(theta) + ", growth = "
                        + ToString(growth) + ", and incrementalResult is "
                        + ToString(incrementalResult);
                    registry.GetRunReport().ReportDebug(msg);
                    //assert(0);
                }
                incrementalResult = 0.0;
            }

            starttime = (*tit).endtime;
            result += incrementalResult;

        }

    }

    return result;
}

//------------------------------------------------------------------------------------
// erynes 2004/05/20 -- This method presently returns the following:
//
// sum_over_a_population's_coalescent_time_intervals(
//  (k_i*(k_i - 1)/theta_p) * ( (tis - tie) + (1/2)*g_p*(tis^2 - tie^2) )
//                                                  )
//
// where
//        k_i is the number of lineages of population p in time interval i
//        theta_p is the present estimate of parameter theta for population p
//        g_p is the present estimate of the growth parameter for population p
//        tis is the timepoint at the start of time interval i
//        tie is the timepoint at the end of time interval i
//
// The population "p" is determined in Wait() by iterating over all
// populations.  That method calls this method with the appropriate
// parameters for the population in question.
//
// This non-public method is used to handle the case in which fabs(g)
// is tiny--that is, less than GROWTH_EPSILON.  It is derived from the
// formula used in Wait by expanding exp(g*t) in a Taylor series
// about g == 0, multiplying this by the term k(k-1)/(g*theta),
// and dropping all terms of order g^2 and above.
// Brief tests with Mathematica suggest that the original formula
// converges very nicely to this linear formula for normal values
// of the timepoints, the difference being less than one part in one
// billion, dropping down to zero difference for g == 0 (using
// L'Hopital's rule).
//
// Please also see the comments for DWait.

double CoalesceGrowPL::lnWaitForTinyGrowth(const double theta, const double growth,
                                           const unsigned long int xpartition,
                                           const list<Interval>& treesummary)
{
    assert(fabs(growth) < GROWTH_EPSILON); // otherwise use Wait

    list<Interval>::const_iterator tit;
    double k, te;

    double result = 0.0;
    // ts starts at 0 and subsequently is the previous cycle's te
    double ts = 0.0;

    for(tit = treesummary.begin(); tit != treesummary.end(); ++tit)
    {
        te = (*tit).endtime;
        k  = (*tit).xpartlines[xpartition];

        result += (k * (k - 1.0) / theta) *
            ((ts - te) + 0.5 * growth * (ts * ts - te * te));
        ts = te;
    }

    return result;
}

//------------------------------------------------------------------------------------
// erynes 2004/03/11 -- This method presently returns the following:
//
// sum_over_populations(sum_over_each_population's_coalescent_time_intervals(
//                              log(2) + g_p*t_i_end - log(theta_p)
//                                                                          )
//                     )
//
// where
//        theta_p is the present estimate of parameter theta for population p
//        g_p is the present estimate of the growth parameter for population p
//        t_i_end is the timepoint at the end of time interval i
//
// This quantity is equal to log(PointProb(G|P)), where PointProb(G|P) is the
// "point" probability of the input genealogy G (parameter "treedata")
// given the parameters P (parameter "param").
// Prob(G|P) = PointProb(G|P) * WaitProb(G|P).
// See CoalesceGrowPL::lnWait() for the wait probability term.

double CoalesceGrowPL::lnPoint (const vector < double >&param, const vector < double >&lparam,
                                const TreeSummary * treedata)

{
    const Interval *  treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval * tit;
    vector < double > :: const_iterator pit = lparam.begin() + m_start;
    vector < double > :: const_iterator pend = lparam.begin () + m_start + m_nPop;
    vector < double > :: const_iterator git = param.begin() + m_growthstart;

    long int xpartition = 0L;
    double result = 0.0;

    for( ; pit != pend; ++pit, ++git, ++xpartition)
    {
        for(tit= treesum; tit != NULL; tit = tit->next)
        {
            if(xpartition==(*tit).oldstatus)
                result +=  LOG2 + (*tit).endtime * (*git) - (*pit);
        }
    }

    // we pull this array only to check for emptyness, it is not used!
    // JDEBUG--:  not correct in presence of migration
    const LongVec2d& picks = treedata->GetCoalSummary()->GetShortPicks();
    if (!picks.empty())                 // partner-picks must be accomodated
    {
        DoubleVec1d timesizes(m_nPop, 0.0);
        ForceParameters fp(unknown_region);
        DoubleVec1d myths(param.begin()+m_start,param.begin()+m_end);
        fp.SetRegionalThetas(myths);
        DoubleVec1d mygs(param.begin()+m_growthstart,
                         param.begin()+m_growthstart+m_nPop);
        fp.SetGrowthRates(mygs);
        // now for every recombination event....
        treesum = treedata->GetRecSummary()->GetLongPoint();
        for (tit = treesum; tit != NULL; tit = tit->next)
        {
            // compute partner-picks denominator
            timesizes = m_timesize->XpartThetasAtT(tit->endtime,fp);
            assert(timesizes.size() == static_cast<DoubleVec1d::size_type>(m_nPop));
            double denominator(std::accumulate(timesizes.begin(), timesizes.end(),0.0));

            // compute partner-picks numerator
            ForceVec::size_type lpforce = 0;
            // assuming no migration and only one local partition force
            long whichtheta =
                m_whichlocalxparts[lpforce][tit->partnerpicks[lpforce]][0];
            result += log(timesizes[whichtheta] / denominator);
        }
    }

#if 0 // JDEBUG--JWARNING this is the logic that was in Jim's original commit
    // compute partner-picks numerator
    for(lpforce = 0; lpforce < m_localpartforces.size(); ++lpforce)
    {
        LocalPartitionForce* pforce(dynamic_cast<LocalPartitionForce*>
                                    (m_localpartforces[lpforce]));
        assert(pforce);  // it wasn't a local partition force?
        DoubleVec1d part_thetas(pforce->SumXPartsToParts(myths, mygs, tit->endtime));
        DoubleVec1d::size_type partition;
        for(partition = 0; partition < part_thetas.size(); ++partition)
        {
            result += log(part_thetas[partition] / denominator);
        }
    }
#endif

    return result;
}

//------------------------------------------------------------------------------------
// erynes 2004/03/11 -- This method currently returns the following:
//
// sum_over_coalescent_time_intervals_of_population_whichparam(
//   -(k_i(k_i - 1)/(g*theta*theta))*(exp(g*t_i_start) - exp(g*t_i_end))
//                                                            )
// where
//   k_i is the number of lineages of this population in time interval i
//   g   is the current estimate of the growth parameter for this pop.
//   theta is the current estimate of the theta parameter for this pop.
//   t_i_start is the timepoint at the beginning of time interval i
//   t_i_end is the timepoint at the end of time interval i
//
// IMPORTANT NOTE:
// Wait() and Point() return log(WaitProb(G|P)) and log(PointProb(G|P)),
// and PostLike::DCalculate() currently assumes that DWait() and DPoint()
// return the derivatives of log(WaitProb(G|P)) and log(PointProb(G|P)),
// instead of returning the derivatives of WaitProb/PointProb directly.
// As long as everything is kept consistent, this is all very good
// for reducing complexity and improving speed.
//
// Note:  CoalesceGrowPL::DWait calculates the partial derivative of Wait
// with respect to theta for the given population ("whichparam").
// CoalescePL::DPoint does the same for Point.
// Partial derivatives with respect to growth are found in GrowPL.

double CoalesceGrowPL::DlnWait (const vector < double >&param, const TreeSummary * treedata,
                                const long int & whichparam)
{
    const double growth = param[whichparam + m_growthstart];
    if(fabs(growth) < GROWTH_EPSILON)
        return DlnWaitForTinyGrowth(param, treedata, whichparam);

    double result = 0.0;
    long int which = whichparam - m_start;
    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list <Interval> :: const_iterator tit=treesum.begin();
    const double theta = param[whichparam];
    double coeff_s, coeff_e, arg1_s, arg1_e;
    double starttime = 0.0;

    for( ; tit != treesum.end(); ++tit)
    {
        double gts = growth * starttime;
        double gte = growth * (*tit).endtime;

        const long int & k = (*tit).xpartlines[which];

        coeff_s = coeff_e = arg1_s = arg1_e = 1.0;

        // IMPORTANT NOTE:  The following over/underflow checking
        // assumes that each starttime and endtime is always >= 0.
        // Hence, e.g., gts < 0 means growth < 0.

        if (gts > 1.0)
        {
            coeff_s = -k * (k - 1.0);
            if (growth > 1.0)
                arg1_s /= growth;
            else
                coeff_s /= growth;
            if (theta > 1.0)
                arg1_s /= theta*theta;
            else
                coeff_s /= theta*theta;
        }
        else if (gts < -1.0)
        {
            arg1_s *= -k * (k - 1.0);
            if (growth < -1.0)
                coeff_s /= growth;
            else
                arg1_s /= growth;
            if (theta > 1.0)
                coeff_s /= theta*theta;
            else
                arg1_s /= theta*theta;
        }
        else // no over/underflow in exp()
            coeff_s = -k * (k - 1.0)/(theta * theta * growth); // arg1_s already set

        // BUGBUG erynes -- what follows is easy to understand,
        // but we should revisit this and reuse some of what we
        // learned about growth and theta above, while considering gts.

        if (gte > 1.0)
        {
            coeff_e = -k * (k - 1.0);
            if (growth > 1.0)
                arg1_e /= growth;
            else
                coeff_e /= growth;
            if (theta > 1.0)
                arg1_e /= theta*theta;
            else
                coeff_e /= theta*theta;
        }
        else if (gte < -1.0)
        {
            arg1_e *= -k * (k - 1.0);
            if (growth < -1.0)
                coeff_e /= growth;
            else
                arg1_e /= growth;
            if (theta > 1.0)
                coeff_e /= theta*theta;
            else
                arg1_e /= theta*theta;
        }
        else // no over/underflow in exp()
            coeff_e = -k * (k - 1.0)/(theta * theta * growth); // arg1_e already set

        double incrementalResult =
            coeff_s * SafeProductWithExp(arg1_s, gts) -
            coeff_e * SafeProductWithExp(arg1_e, gte);

        if (incrementalResult < 0.0)
        {
            // We have precision problems,
            // because we expect te > ts always,
            // and hence gte > gts always.
            // During debugging, we have encountered precision problems
            // in which gdb claims gts == gte and arg1_s == arg1_e
            // but incrementalResult < 0, typically reflecting
            // a difference near the 15th digit.
            // (Note that te == ts corresponds to two events
            // occurring simultaneously--e.g., k = 1 for a really long
            // time in each of two populations, then one finally
            // migrates to the other, then they coalesce instantaneously.)
            // It should be safe to set incrementalResult to zero in cases of
            // roundoff error.  If incrementalResult is negative,
            // but _not_ tiny compared with either of the terms
            // whose difference yielded incrementalResult,
            // then something is _very_ wrong--probably te < ts,
            // which should never happen.  We assert in this case;
            // this will enable us to see the problem, but the end user
            // will not see the assertion, unless s/he's running in debug mode.
            if (incrementalResult/fabs(coeff_s * SafeProductWithExp(arg1_s, gts))
                < -1e-10)
                assert(0);
            incrementalResult = 0.0;
        }

        starttime = (*tit).endtime;
        result += incrementalResult;

    }

    return result;
}

//------------------------------------------------------------------------------------
// erynes 2004/05/20 -- This method currently returns the following:
//
// sum_over_coalescent_time_intervals_of_population_whichparam(
//   -(k_i(k_i - 1)/(theta*theta))*( (t_i_start - t_i_end) +
//                                   (1/2)*g*(t_i_start^2 - t_i_end^2) )
//                                                            )
// where
//   k_i is the number of lineages of this population in time interval i
//   g   is the current estimate of the growth parameter for this pop.
//   theta is the current estimate of the theta parameter for this pop.
//   t_i_start is the timepoint at the beginning of time interval i
//   t_i_end is the timepoint at the end of time interval i
//
// This non-public method is used to handle the case in which fabs(g)
// is tiny--that is, less than GROWTH_EPSILON.  It is derived from the
// formula used in DWait by expanding exp(g*t) in a Taylor series
// about g == 0, multiplying this by the term -k(k-1)/(g*theta*theta),
// and dropping all terms of order g^2 and above.
// Brief tests with Mathematica suggest that the original formula
// converges very nicely to this linear formula for normal values
// of the timepoints, the difference being less than one part in one
// billion, dropping down to zero difference for g == 0 (using
// L'Hopital's rule).
//
// Please also see the comments for DWait.

double CoalesceGrowPL::DlnWaitForTinyGrowth(const vector < double >&param, const TreeSummary * treedata,
                                            const long int & whichparam)
{
    const double growth = param[whichparam + m_growthstart];
    assert(fabs(growth) < GROWTH_EPSILON); // otherwise use DWait()
    const double theta = param[whichparam];
    double k, te;
    long int which = whichparam - m_start;
    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator tit = treesum.begin();

    double result = 0.0;
    double ts = 0.0;

    for( ; tit != treesum.end(); ++tit)
    {
        k  = (*tit).xpartlines[which];
        te = (*tit).endtime;

        result += -(k * (k - 1.0) / (theta * theta)) *
            ((ts - te) + 0.5 * growth * (ts * ts - te * te));
        ts = te;
    }

    return result;
}

//------------------------------------------------------------------------------------
// Mary 2007/08/13 -- added support for disease w/ recombination
//
// JDEBUG--:  only in absence of migration

double CoalesceGrowPL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long &whichparam)
{
    long which = whichparam - m_start;
    const vector<double>& ncoal = treedata->GetCoalSummary()->GetShortPoint();

    double answ(-ncoal[which] / param[whichparam]);

    const LongVec2d& picks = treedata->GetCoalSummary()->GetShortPicks();
    if (picks.empty()) return answ;
    else  {
        DoubleVec1d timesizes(m_nPop, 0.0);
        ForceParameters fp(unknown_region);
        DoubleVec1d myths(param.begin()+m_start,param.begin()+m_end);
        fp.SetRegionalThetas(myths);
        DoubleVec1d mygs(param.begin()+m_growthstart,
                         param.begin()+m_growthstart+m_nPop);
        fp.SetGrowthRates(mygs);
        // now for every recombination event....
        const Interval* reclist  = treedata->GetRecSummary()->GetLongPoint();
        const Interval* tit;
        for (tit = reclist; tit != NULL; tit = tit->next)
        {
            // compute partner-picks denominator
            timesizes = m_timesize->XpartThetasAtT(tit->endtime,fp);
            assert(timesizes.size() == static_cast<DoubleVec1d::size_type>(m_nPop));
            double denominator(std::accumulate(timesizes.begin(), timesizes.end(),0.0));

            // compute partner-picks numerator
            // ALL assuming no migration and only one local partition force
            // therefore xpart == part
            answ -= 1.0/denominator;
            if (tit->partnerpicks[0] == which)
            {
                answ += 1.0/timesizes[which];
            }
        }
    }

    return answ;
} // CoalesceGrowPL::DlnPoint

//------------------------------------------------------------------------------------
// GrowPL: exponential growth forces specific waiting time (=0) and
//              point probabilites and its derivatives
// LONG DATA SUMMARIES

double GrowPL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    // Growth is a modifier of theta, hence its contribution
    // is included in CoalesceGrowPL::Wait.
    return 0.0;
}

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

double GrowPL::lnPoint (const vector < double >&param, const vector < double >&lparam, const TreeSummary * treedata)
{
    // Growth is a modifier of theta, hence its contribution
    // is included in CoalesceGrowPL::Point.
    return 0.0;
}

//------------------------------------------------------------------------------------
// erynes 2004/04/09 -- This method currently returns the following:
//
// sum_over_coalescent_time_intervals_of_population_whichparam(
//   -(k_i(k_i - 1)/(g*g*theta))*(
//             (1 - g*t_i_start)*exp(g*t_i_start) -
//             (1 - g*t_i_end)  *exp(g*t_i_end)
//                               )
//                                                            )
// where
//   k_i is the number of lineages of this population in time interval i
//   g   is the current estimate of the growth parameter for this pop.
//   theta is the current estimate of the theta parameter for this pop.
//   t_i_start is the timepoint at the beginning of time interval i
//   t_i_end is the timepoint at the end of time interval i

double GrowPL::DlnWait (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    const double growth = param[whichparam];
    if (fabs(growth) < GROWTH_EPSILON)
        return DlnWaitForTinyGrowth(param, treedata, whichparam);

    const long int whicht = m_thetastart + whichparam - m_start;
    const double theta = param[whicht];
    assert(theta); // we divide by it below
    double gte, gts, arg1_s, arg1_e, coeff_s, coeff_e;
    long int k;

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list <Interval> :: const_iterator tit = treesum.begin();

    double result = 0.0;
    double starttime = 0.0;

    for( ; tit != treesum.end(); ++tit)
    {
        gts = growth * starttime;
        gte = growth * (*tit).endtime;

        k = (*tit).xpartlines[whicht];

        arg1_s = 1.0 - gts;
        arg1_e = 1.0 - gte;
        coeff_s = 1.0;
        coeff_e = 1.0;

        // IMPORTANT NOTE:  The following over/underflow checking
        // assumes that each starttime and endtime is always >= 0.
        // Hence, e.g., gts < 0 means growth < 0.

        if (gts > 1.0)
        {
            coeff_s = -k * (k - 1.0);
            if (growth > 1.0)
                arg1_s /= growth*growth;
            else
                coeff_s /= growth*growth;
            if (theta > 1.0)
                arg1_s /= theta;
            else
                coeff_s /= theta;
        }
        else if (gts < -1.0)
        {
            arg1_s *= -k*(k-1.0);
            if (growth < -1.0)
                coeff_s /= growth*growth;
            else
                arg1_s /= growth*growth;
            if (theta > 1.0)
                coeff_s /= theta;
            else
                arg1_s /= theta;
        }
        else // no over/underflow in exp()
            coeff_s = -k*(k-1.0)/(theta*growth*growth); // arg1_s already set

        // BUGBUG erynes -- this is easy to follow,
        // but we should revisit it and possibly reuse
        // the growth and theta checks that we made for the gts case.

        if (gte > 1.0)
        {
            coeff_e = -k*(k-1.0);
            if (growth > 1.0)
                arg1_e /= growth*growth;
            else
                coeff_e /= growth*growth;
            if (theta > 1.0)
                arg1_e /= theta;
            else
                coeff_e /= theta;
        }
        else if (gte < -1.0)
        {
            arg1_e *= -k*(k-1.0);
            if (growth < -1.0)
                coeff_e /= growth*growth;
            else
                arg1_e /= growth*growth;
            if (theta > 1.0)
                coeff_e /= theta;
            else
                arg1_e /= theta;
        }
        else // no over/underflow in exp()
            coeff_e = -k*(k-1.0)/(theta*growth*growth); // arg1_e already set

        starttime = (*tit).endtime;
        result +=  coeff_s * SafeProductWithExp(arg1_s, gts)
            - coeff_e * SafeProductWithExp(arg1_e, gte);
    }

    return result;
}

//------------------------------------------------------------------------------------
// erynes 2004/05/20 -- This method currently returns the following:
//
// sum_over_coalescent_time_intervals_of_population_whichparam(
//    (k_i(k_i - 1)/(2*theta))*(
//                  tis^2 * (1 + g*tis)  -  tie^2 * (1 + g*tie)
//                             )
//                                                            )
// where
//   k_i is the number of lineages of this population in time interval i
//   g   is the current estimate of the growth parameter for this pop.
//   theta is the current estimate of the theta parameter for this pop.
//   tis is the timepoint at the beginning of time interval i
//   tie is the timepoint at the end of time interval i
//
// This non-public method is used to handle the case in which fabs(g)
// is tiny--that is, less than GROWTH_EPSILON.  It is derived from the
// formula used in DWait by expanding exp(g*t) in a Taylor series
// about g == 0, multiplying through in the original formula,
// and dropping all terms of order g^2 and above.
// Brief tests with Mathematica suggest that the original formula
// converges very nicely to this linear formula for normal values
// of the timepoints, the difference being less than one part in one
// billion, dropping down to zero difference for g == 0 (using
// L'Hopital's rule).
//
// Please also see the comments for DWait.

double GrowPL::DlnWaitForTinyGrowth(const vector < double >&param, const TreeSummary * treedata,
                                    const long int & whichparam)
{
    const double growth = param[whichparam];
    assert(fabs(growth) < GROWTH_EPSILON); // otherwise use DWait

    const long int whicht = m_thetastart + whichparam - m_start;
    const double theta = param[whicht];
    assert(theta); // we divide by it below
    double k, te;
    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list <Interval> :: const_iterator tit = treesum.begin();

    double result = 0.0;
    double ts = 0.0;

    for( ; tit != treesum.end(); ++tit)
    {
        te = (*tit).endtime;
        k = (*tit).xpartlines[whicht];

        result += (k * (k - 1.0) / (2.0 * theta)) *
            (ts * ts * (1.0 + growth * ts) - te * te * (1.0 + growth * te));
        ts = te;
    }

    return result;
}

//------------------------------------------------------------------------------------
// This method presently returns:
//
// sum_over_time_intervals_of_this_population's_coalescences(endtime)
//
// This is the derivative of log(Point) with respect to growth.

double GrowPL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    const Interval * treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval * tit;
    double result = 0.0 ;
    // indicator to Theta, we need this to get the correct
    // time intervals in treesum
    long int whicht = m_thetastart + whichparam - m_start;

    for(tit = treesum; tit != NULL; tit = tit->next)
    {
        if(whicht == (*tit).oldstatus)
            result += (*tit).endtime;
    }

    return result;
}

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

double DivMigPL::lnWait(const vector<double>& param, const TreeSummary* treedata)
{
    const vector<double>& mWait = treedata->GetDivMigSummary()->GetShortWait();

    // sum_pop(SM_pop *kt)
    return  -1.0 * inner_product (mWait.begin (),
                                  mWait.end (), param.begin() + m_start, 0.0);

}

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

double DivMigPL::lnPoint(const vector<double>& param, const vector<double>& lparam, const TreeSummary* treedata)
{
    //  Sum_j(Sum_i(migevent[j,i] * log(M[j,i]))
    const vector<double>& nmig = treedata->GetDivMigSummary()->GetShortPoint();

    return  inner_product (nmig.begin(), nmig.end(), lparam.begin() + m_start,
                           0.0);
}

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

// RecombinePL: recombination forces specific waiting time and
//              point probabilites and its derivatives

double RecombinePL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    const vector<double>& rWait = treedata->GetRecSummary()->GetShortWait();
    return  -rWait[0] * *(param.begin () + m_start);
}

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

double RecombinePL::lnPoint (const vector < double >&param, const vector < double >&lparam,
                             const TreeSummary * treedata)
{
    //  Sum_j(recevents * log(r))
    const vector<double>& nrec = treedata->GetRecSummary()->GetShortPoint();
    return nrec[0] * *(lparam.begin () + m_start);
}

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

double RecombinePL::DlnWait (const vector < double >&param, const TreeSummary * treedata,
                             const long int & whichparam)
{
    const vector<double>& rWait = treedata->GetRecSummary()->GetShortWait();
    return -rWait[0];
}

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

double RecombinePL::DlnPoint (const vector < double >&param, const TreeSummary * treedata,
                              const long int & whichparam)
{
    const vector<double>& nrecs = treedata->GetRecSummary()->GetShortPoint();
    if(nrecs[0] == 0)
    {
        return 0.0;
    }
    else
    {
        return SafeDivide(nrecs[0], *(param.begin() + m_start));
    }
}

//------------------------------------------------------------------------------------
// SelectPL: selection forces specific waiting time and
//              point probabilites and its derivatives

double SelectPL::lnWait (const vector < double >&param, const TreeSummary * treedata)
{
    return 0.0;
}

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

double SelectPL::lnPoint (const vector < double >&param, const vector < double >&lparam, const TreeSummary * treedata)
{
    return 0.0;
}

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

double SelectPL::DlnWait (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    return 0.0;
}

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

double SelectPL::DlnPoint (const vector < double >&param, const TreeSummary * treedata, const long int & whichparam)
{
    return 0.0;
}

//------------------------------------------------------------------------------------
// This function computes the "waiting time" contribution to the log-likelihood.
// It computes and returns
//
// sum_over_all_timeIntervals_i( -kA_i*(kA_i - 1)*dtauA_i/thetaA0 -ka_i*(ka_i - 1)*dtaua_i/thetaa0 ),
//
// where "A" is the favored allele, "a" is the disfavored allele, "k" is the number of lineages
// in population "A" or "a," and dtau is the length of the time interval in the "Kingman reference frame"
// (i.e., time running as usual and pop. size remaining constant).  The transformations from
// "magic time" dtau to real time "dt," where dt = te - ts = endtime - starttime, are:
//
// dtau_A = (theta_A0*dt + (theta_a0/s)*(exp(s*te) - exp(s*ts)))/(theta_A0 + theta_a0),
// dtau_a = (theta_a0*dt - (theta_A0/s)*(exp(-s*te) - exp(-s*ts)))/(theta_A0 + theta_a0).

double CoalesceLogisticSelectionPL::lnWait(const vector<double>& param, const TreeSummary *treedata)
{
    if (2 != m_nPop)
    {
        string msg = "CoalesceLogisticSelectionPL::lnWait() called with m_nPop = ";
        msg += ToString(m_nPop);
        msg += ".  m_nPop must equal 2, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }
    double mu_into_A_from_a(param[3]), mu_into_a_from_A(param[4]);

    if (0 != param[2] || 0 != param[5] || 0 == param[3] || 0 == param[4])
        throw implementation_error("Error parsing disease rates in CoalesceLogisticSelectionPL::lnWait()");

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;

    double result(0.0), s(param[m_s_is_here]);
    double theta_A0(param[m_start]), theta_a0(param[m_start+1]);
    double totalTheta_0(theta_A0 + theta_a0);

    if (theta_a0 <= 0.0 || theta_A0 <= 0.0)
    {
        string msg = "CoalesceLogisticSelectionPL::lnWait(), received an invalid Theta value ";
        msg += "(theta_A0 = " + ToString(theta_A0) + ", theta_a0 = " + ToString(theta_a0) + ").";
        throw impossible_error(msg);
    }

    if (fabs(s) < LOGISTIC_SELECTION_COEFFICIENT_EPSILON)
        return lnWaitForTinyLogisticSelectionCoefficient(param, treedata);

    if (fabs(s) >= DBL_BIG)
        return -DBL_BIG; // unavoidable overflow

    double factor_A(theta_a0/(s*theta_A0)), factor_a(theta_A0/(-s*theta_a0));
    double t_e, t_s(0.0), dt, e_toThe_sts(1.0), e_toThe_ste, term;
    double term_A(0.0), term_a(0.0);

    for (treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        dt = t_e - t_s;
        term_A = term_a = 0.0;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        if (s > 0.0)
        {
            // Use overflow protection for term_A; if term_a underflows, that's okay.

            if (s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor_A, s*dt); // coefficient guaranteed finite and positive
                if (term >= DBL_BIG)
                    return -DBL_BIG; // unavoidable overflow
                e_toThe_ste = (term/factor_A)*e_toThe_sts;
                if (k_A > 0.0)
                {
                    term_A = k_A*(term - factor_A)*e_toThe_sts*((k_A - 1.0)/totalTheta_0 + mu_into_A_from_a)
                        + k_A*(k_A - 1.0)*dt/totalTheta_0;
                    // Note:  A factor of theta_A0 cancels in the "dt" term.
                    if (term_A >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_A > 0.0)
                {
                    term_A = k_A*factor_A*s*dt*e_toThe_sts*((k_A - 1.0)/totalTheta_0 + mu_into_A_from_a)
                        + k_A*(k_A - 1.0)*dt/totalTheta_0;
                    // Note:  A factor of theta_A0 cancels in the "dt" term.
                    if (term_A >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            if (k_a > 0.0)
                term_a = k_a * factor_a
                    * (1.0/e_toThe_ste - 1.0/e_toThe_sts)
                    * ((k_a - 1.0)/totalTheta_0 + mu_into_a_from_A)
                    + k_a*(k_a - 1.0)*dt/totalTheta_0;
            // Note:  A factor of theta_a0 cancels in the "dt" term.
        }
        else // s < 0
        {
            // Use overflow protection for term_a; if term_A underflows, that's okay.

            if (-s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor_a,-s*dt); // coefficient guaranteed finite and positive
                if (term >= DBL_BIG)
                    return -DBL_BIG; // unavoidable overflow
                e_toThe_ste = e_toThe_sts/(term/factor_a);
                if (k_a > 0.0)
                {
                    term_a = k_a*((term - factor_a)/e_toThe_sts)*((k_a - 1.0)/totalTheta_0 + mu_into_a_from_A)
                        + k_a*(k_a - 1.0)*dt/totalTheta_0;
                    // Note:  A factor of theta_a0 cancels in the "dt" term.
                    if (term_a >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (-s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_a > 0.0)
                {
                    term_a = k_a*(factor_a*(-s)*dt/e_toThe_sts)*((k_a - 1.0)/totalTheta_0 + mu_into_a_from_A)
                        + k_a*(k_a - 1.0)*dt/totalTheta_0;
                    // Note:  A factor of theta_a0 cancels in the "dt" term.
                    if (term_a >= DBL_BIG)
                        return -DBL_BIG; // unavoidable overflow
                }
            }
            if (k_A > 0.0)
                term_A = k_A*factor_A*(e_toThe_ste - e_toThe_sts)*((k_A - 1.0)/totalTheta_0 + mu_into_A_from_a)
                    + k_A*(k_A - 1.0)*dt/totalTheta_0;
            // Note:  A factor of theta_a0 cancels in the "dt" term.
        }

        result -= term_A + term_a;
        if (result <= -DBL_BIG)
            return -DBL_BIG;

        t_s = t_e;
        e_toThe_sts = e_toThe_ste;
    }

    return result;
}

//------------------------------------------------------------------------------------
// See comment for lnWait(), above.
// This function expands lnWait() in a Taylor series about s = 0,
// so that we can smoothly handle the factor of 1/s in lnWait().
// This version retains only the first two terms in
// (exp(s*te) - exp(s*ts))/s = (te-ts) + (s/2!)(te^2 - ts^2) + (s^2/3!)(te^3 - ts^3) + ....

double CoalesceLogisticSelectionPL::lnWaitForTinyLogisticSelectionCoefficient(const vector<double>& param,
                                                                              const TreeSummary* treedata)
{
    double mu_into_A_from_a(param[3]), mu_into_a_from_A(param[4]);

    if (0 != param[2] || 0 != param[5] || 0 == param[3] || 0 == param[4])
        throw implementation_error("Error parsing disease rates in "
                                   "CoalesceLogisticSelectionPL::lnWaitForTinyLogisticSelectionCoefficient()");

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;

    double result(0.0), t_e, t_s(0.0), s(param[m_s_is_here]), dt, term;
    double theta_A0(param[m_start]), theta_a0(param[m_start+1]), totalTheta_0(theta_A0+theta_a0);

    for (treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a
        dt = t_e - t_s;

        term = 0.5*s*(t_e*t_e - t_s*t_s);

        result += -(k_A*dt/theta_A0)*((k_A - 1.0) + theta_a0*mu_into_A_from_a)
            - (k_A*theta_a0*term/theta_A0)*((k_A - 1.0)/totalTheta_0 + mu_into_A_from_a);

        result += -(k_a*dt/theta_a0)*((k_a - 1.0) + theta_A0*mu_into_a_from_A)
            + (k_a*theta_A0*term/theta_a0)*((k_a - 1.0)/totalTheta_0 + mu_into_a_from_A);

        t_s = t_e;
    }

    return result;
}

//------------------------------------------------------------------------------------
// This function computes the "point" contribution to the log-likelihood.
// It computes and returns
//
// sum_over_all_coalescent_times_in_A( ln(2/theta_A(te_A)) )
//  + sum_over_all_coalescent_times_in_a( ln(2/theta_a(te_a)) ),
//
// where te_A is the time (te = "endtime") at the bottom of a time interval
// which ends with a coalescence in the population with allele "A," and
// similarly for the pop. with allele "a."
// theta_A(t) + theta_a(t) = theta_A(0) + theta_a(0) for all times t;
// that is, the overall population size is assumed constant.  We have
//
// theta_A(t) = ((theta_A0+theta_a0)/(theta_A0+theta_a0*exp(s*t)))*theta_A0,
// theta_a(t) = ((theta_A0+theta_a0)/(theta_A0+theta_a0*exp(s*t)))*theta_a0*exp(s*t).

double CoalesceLogisticSelectionPL::lnPoint(const vector<double>& param, const vector<double>& lparam,
                                            const TreeSummary* treedata)
{
    if (2 != m_nPop)
    {
        string msg = "CoalesceLogisticSelectionPL::lnPoint() called with m_nPop = ";
        msg += ToString(m_nPop);
        msg += ".  m_nPop must equal 2, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }

    const Interval *treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval *pTreesum;
    double theta_A0(param[0]), theta_a0(param[1]), logTheta_A0(lparam[0]),
        logTheta_a0(lparam[1]), logThetaTotal_0, s(param[m_s_is_here]);

    if (theta_A0 > 0.0 && theta_a0 > 0.0)
        logThetaTotal_0 = log(theta_A0 + theta_a0);
    else
    {
        // This really, really should never ever happen.
        string msg = "CoalesceLogisticSelectionPL::lnPoint(), received an invalid value ";
        msg += "(theta_A0 = ";
        msg += ToString(theta_A0);
        msg += " and theta_a0 = " + ToString(theta_a0) + "); s = ";
        msg += ToString(s) + ".";
        throw impossible_error(msg);
    }

    double result = 0.0;

    for(pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
    {
        double thetaDenominator = theta_A0 +
            SafeProductWithExp(theta_a0, s * pTreesum->endtime); // sum is guaranteed to be > 0

        if (thetaDenominator < DBL_BIG)
            result += LOG2 - logThetaTotal_0 + log(thetaDenominator);
        else
            result += LOG2 - logThetaTotal_0 + logTheta_a0 + s*pTreesum->endtime;
        if (0L == pTreesum->oldstatus)
            result -= logTheta_A0;
        else if (1L == pTreesum->oldstatus)
            result -= logTheta_a0 + (pTreesum->endtime * s);
        else
        {
            string msg = "CoalesceLogisticSelectionPL::lnPoint(), received a TreeSummary ";
            msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
            msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescene in either ";
            msg += "the population with allele A or the population with allele a.";
            throw implementation_error(msg);
        }
    }

    return result;
}

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

double CoalesceLogisticSelectionPL::DlnWait(const vector<double>& param, const TreeSummary * treedata,
                                            const long int & whichparam)
{
    if (2 != m_nPop)
    {
        string msg = "CoalesceLogisticSelectionPL::DlnWait() called with m_nPop = ";
        msg += ToString(m_nPop);
        msg += ".  m_nPop must equal 2, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }

    double s(param[m_s_is_here]);

    if (fabs(s) < LOGISTIC_SELECTION_COEFFICIENT_EPSILON)
        return DlnWaitForTinyLogisticSelectionCoefficient(param, treedata, whichparam);

    bool derivativeWithRespectTo_A;

    if (0 == whichparam)
        derivativeWithRespectTo_A = true;
    else if (1 == whichparam)
        derivativeWithRespectTo_A = false;
    else
    {
        string msg = "CoalesceLogisticSelectionPL::DlnWait() called with whichparam = ";
        msg += ToString(whichparam);
        msg += ".  \"whichparam\" must equal 0 or 1, reflecting one population with the favored ";
        msg += "allele and one population with the disfavored allele.";
        throw implementation_error(msg);
    }

    if (s >= DBL_BIG)
        return derivativeWithRespectTo_A ? DBL_BIG : -DBL_BIG;
    if (s <= -DBL_BIG)
        return derivativeWithRespectTo_A ? -DBL_BIG : DBL_BIG;

    if (0 != param[2] || 0 != param[5] || 0 == param[0] || 0 == param[1] ||
        0 == param[3] || 0 == param[4])
        throw implementation_error("Bad param received by CoalesceLogisticSelectionPL::DlnWait().");

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;
    double theta_A0(param[0]), theta_a0(param[1]), totalTheta_0(theta_A0+theta_a0);
    double term_A(0.0), term_a(0.0), term, result(0.0);
    double t_s(0.0),t_e,dt,e_toThe_sts(1.0),e_toThe_ste;
    double factor = 1.0/(totalTheta_0*totalTheta_0);
    double thetaFactor_A = theta_a0/(theta_A0*theta_A0);
    double thetaFactor_a = theta_A0/(theta_a0*theta_a0);
    double mu_into_A_from_a(param[3]), mu_into_a_from_A(param[4]);

    for (treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        dt = t_e - t_s;
        term_A = term_a = 0.0;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        if (s > 0.0)
        {
            // Use overflow protection for term_A; if term_a underflows, that's okay.

            if (s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor/s,s*dt); // quotient guaranteed finite and positive
                if (term >= DBL_BIG)
                    return derivativeWithRespectTo_A ? DBL_BIG : -DBL_BIG; // unavoidable overflow
                e_toThe_ste = (term/(factor/s))*e_toThe_sts;
                if (k_A > 0)
                {
                    if (derivativeWithRespectTo_A)
                        term_A = (term - factor/s)*thetaFactor_A*e_toThe_sts;
                    else
                        term_A = (term - factor/s)*e_toThe_sts;
                    if (term_A >= DBL_BIG)
                        return derivativeWithRespectTo_A ? DBL_BIG : -DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_A > 0)
                {
                    if (derivativeWithRespectTo_A)
                        term_A = factor*dt*thetaFactor_A*e_toThe_sts;
                    else
                        term_A = factor*dt*e_toThe_sts;
                    if (term_A >= DBL_BIG)
                        return derivativeWithRespectTo_A ? DBL_BIG : -DBL_BIG; // unavoidable overflow
                }
            }
            if (k_a > 0)
            {
                if (derivativeWithRespectTo_A)
                    term_a = (factor/(-s))*(1.0/e_toThe_ste - 1.0/e_toThe_sts);
                else
                    term_a = (factor/(-s))*thetaFactor_a*(1.0/e_toThe_ste - 1.0/e_toThe_sts);
            }
        }

        else // s < 0
        {
            // Use overflow protection for dtau_a; if dtau_A underflows, that's okay.

            if (-s*dt > LOG_ONEPLUSEPSILON)
            {
                term = SafeProductWithExp(factor/(-s),-s*dt);// quotient guaranteed finite and positive
                if (term >= DBL_BIG)
                    return derivativeWithRespectTo_A ? -DBL_BIG : DBL_BIG; // unavoidable overflow
                e_toThe_ste = e_toThe_sts/(term/(factor/(-s)));
                if (k_a > 0)
                {
                    if (derivativeWithRespectTo_A)
                        term_a = (term - factor/(-s))/e_toThe_sts;
                    else
                        term_a = (term - factor/(-s))*thetaFactor_a/e_toThe_sts;
                    if (term_a >= DBL_BIG)
                        return derivativeWithRespectTo_A ? -DBL_BIG : DBL_BIG; // unavoidable overflow
                }
            }
            else
            {
                if (-s*dt >= numeric_limits<double>::epsilon())
                    e_toThe_ste = (1.0 + s*dt)*e_toThe_sts;
                else
                    e_toThe_ste = SafeProductWithExp(1.0,s*t_e);
                if (k_a > 0)
                {
                    if (derivativeWithRespectTo_A)
                        term_a = factor*dt/e_toThe_sts;
                    else
                        term_a = factor*dt*thetaFactor_a/e_toThe_sts;
                    if (term_a >= DBL_BIG)
                        return derivativeWithRespectTo_A ? -DBL_BIG : DBL_BIG; // unavoidable overflow
                }
            }
            if (k_A > 0)
            {
                if (derivativeWithRespectTo_A)
                    term_A = (factor/s)*thetaFactor_A*(e_toThe_ste - e_toThe_sts);
                else
                    term_A = (factor/s)*(e_toThe_ste - e_toThe_sts);
            }
        }

        if (derivativeWithRespectTo_A)
            result += k_A*term_A*((k_A-1.0)*(theta_A0+totalTheta_0) + mu_into_A_from_a/factor)
                - k_a*term_a*((k_a-1.0) + mu_into_a_from_A/(factor*theta_a0))
                + dt*factor*(k_A*(k_A-1.0) + k_a*(k_a-1.0));
        else
            result += -k_A*term_A*((k_A-1.0) + mu_into_A_from_a/(factor*theta_A0))
                + k_a*term_a*((k_a-1.0)*(theta_a0+totalTheta_0) + mu_into_a_from_A/factor)
                + dt*factor*(k_A*(k_A-1.0) + k_a*(k_a-1.0));

        if (result >= DBL_BIG)
            return DBL_BIG;
        if (result <= -DBL_BIG)
            return -DBL_BIG;

        t_s = t_e;
        e_toThe_sts = e_toThe_ste;
    }

    return result;
}

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

double CoalesceLogisticSelectionPL::DlnWaitForTinyLogisticSelectionCoefficient(const vector<double>& param,
                                                                               const TreeSummary *treedata,
                                                                               const long int & whichparam)
{
    bool derivativeWithRespectTo_A;

    if (0 == whichparam)
        derivativeWithRespectTo_A = true;
    else if (1 == whichparam)
        derivativeWithRespectTo_A = false;
    else
    {
        string msg = "CoalesceLogisticSelectionPL::DlnWaitForTinyLogisticSelectionCoefficient() ";
        msg += "called with whichparam = ";
        msg += ToString(whichparam);
        msg += ".  \"whichparam\" must equal 0 or 1, reflecting one subpopulation with the favored ";
        msg += "allele and one subpopulation with the disfavored allele.";
        throw implementation_error(msg);
    }

    if (0 != param[2] || 0 != param[5] || 0 == param[0] || 0 == param[1] ||
        0 == param[3] || 0 == param[4])
        throw implementation_error
            ("Bad param received by CoalesceLogisticSelectionPL::DlnWaitForTinyLogisticSelectionCoefficient().");

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;
    double theta_A0(param[0]), theta_a0(param[1]), totalTheta_0(theta_A0+theta_a0);
    double t_e, t_s(0.0), resultTimesDenominator(0.0), s(param[m_s_is_here]);
    double mu_into_A_from_a(param[3]), mu_into_a_from_A(param[4]);
    double denominator(totalTheta_0*totalTheta_0), term_A, term_a;
    double factor_A(theta_a0/(theta_A0*theta_A0)), factor_a(theta_A0/(theta_a0*theta_a0));

    for(treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        double dt = t_e - t_s;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        if (derivativeWithRespectTo_A)
        {
            term_A = k_A*factor_A*(dt + 0.5*s*(t_e*t_e - t_s*t_s))*((k_A-1.0)*(theta_A0+totalTheta_0)
                                                                    + mu_into_A_from_a*denominator);
            term_a = -k_a*(dt - 0.5*s*(t_e*t_e - t_s*t_s))*((k_a-1.0) + mu_into_a_from_A*denominator);
        }
        else
        {
            term_A = -k_A*(dt + 0.5*s*(t_e*t_e - t_s*t_s))*((k_A-1.0) + mu_into_A_from_a*denominator);
            term_a = k_a*factor_a*(dt = 0.5*s*(t_e*t_e - t_s*t_s))*((k_a-1.0)*(theta_a0+totalTheta_0)
                                                                    + mu_into_a_from_A*denominator);
        }

        resultTimesDenominator += term_A + term_a + dt*(k_A*(k_A-1.0) + k_a*(k_a-1.0));

        t_s = t_e;
    }

    return resultTimesDenominator/denominator;
}

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

double CoalesceLogisticSelectionPL::DlnPoint(const vector<double>& param, const TreeSummary *treedata,
                                             const long int & whichparam)
{
    if (2 != m_nPop)
    {
        string msg = "CoalesceLogisticSelectionPL::DlnPoint() called with m_nPop = ";
        msg += ToString(m_nPop);
        msg += ".  m_nPop must equal 2, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }

    bool derivativeWithRespectTo_A;

    if (0 == whichparam)
        derivativeWithRespectTo_A = true;
    else if (1 == whichparam)
        derivativeWithRespectTo_A = false;
    else
    {
        string msg = "CoalesceLogisticSelectionPL::DlnPoint() called with whichparam = ";
        msg += ToString(whichparam);
        msg += ".  \"whichparam\" must equal 0 or 1, reflecting one population with the major ";
        msg += "allele and one population with the minor allele.";
        throw implementation_error(msg);
    }

    const Interval *treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval *pTreesum;
    double theta_A0(param[0]), theta_a0(param[1]), s(param[m_s_is_here]),
        totalTheta_0(theta_A0+theta_a0), theta_A0squared(theta_A0*theta_A0),
        theta_a0squared(theta_a0*theta_a0), result(0.0), x;

    if (theta_a0 <= 0.0 || theta_A0 <= 0.0)
    {
        string msg = "CoalesceLogisticSelectionPL::DlnPoint(), received an invalid Theta value ";
        msg += "(theta_A0 = " + ToString(theta_A0) + ", theta_a0 = " + ToString(theta_a0) + ").";
        throw impossible_error(msg);
    }

    // First, add up the contributions from the coalescent events.

    if (s >= 0.0)
    {
        for (pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
        {
            double t_e = pTreesum->endtime;
            if (derivativeWithRespectTo_A)
            {
                if (0L == pTreesum->oldstatus) // coal. in allele A
                {
                    x = SafeProductWithExp(theta_a0*(theta_A0+totalTheta_0),s*t_e);
                    if (x < DBL_BIG)
                        result -= (theta_A0squared + x) /
                            ((theta_A0*totalTheta_0)*(theta_A0 + x/(theta_A0+totalTheta_0)));
                    else
                        result -= (theta_A0 + totalTheta_0)/(theta_A0*totalTheta_0);
                }
                else if (1L == pTreesum->oldstatus) // coal. in allele a
                {
                    x = SafeProductWithExp(theta_a0,s*t_e);
                    if (x < DBL_BIG)
                        result -= (x - theta_a0)/(totalTheta_0*(theta_A0 + x));
                    else
                        result -= 1.0/totalTheta_0;
                }
                else
                {
                    string msg = "CoalesceLogisticSelectionPL::DlnPoint(), received a TreeSummary ";
                    msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
                    msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescene in either ";
                    msg += "the population with allele A or the population with allele a.";
                    throw implementation_error(msg);
                }
            }
            else // derivative with repect to theta_a0
            {
                if (0L == pTreesum->oldstatus) // coal. in allele A
                {
                    x = SafeProductWithExp(theta_A0,s*t_e);
                    if (x < DBL_BIG)
                        result += (x - theta_A0) /
                            (totalTheta_0*(theta_A0 + x*theta_a0/theta_A0));
                    else
                        result += theta_A0/(theta_a0*totalTheta_0);
                }
                else if (1L == pTreesum->oldstatus) // coal. in allele a
                {
                    x = SafeProductWithExp(theta_a0squared,s*t_e);
                    if (x < DBL_BIG)
                        result -= (theta_A0*(theta_a0+totalTheta_0) + x) /
                            (totalTheta_0*(theta_A0*theta_a0 + x));
                    else
                        result -= 1.0/totalTheta_0;
                }
                else
                {
                    string msg = "CoalesceLogisticSelectionPL::DlnPoint(), received a TreeSummary ";
                    msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
                    msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescene in either ";
                    msg += "the population with allele A or the population with allele a.";
                    throw implementation_error(msg);
                }
            }
        }
    }
    else // s < 0
    {
        for (pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
        {
            double t_e = pTreesum->endtime;
            if (derivativeWithRespectTo_A)
            {
                if (0L == pTreesum->oldstatus) // coal. in allele A
                {
                    x = SafeProductWithExp(theta_A0squared,-s*t_e);
                    if (x < DBL_BIG)
                        result -= (theta_a0*(theta_A0+totalTheta_0) + x) /
                            (totalTheta_0*(theta_A0*theta_a0 + x));
                    else
                        result -= 1.0/totalTheta_0;
                }
                else if (1L == pTreesum->oldstatus) // coal. in allele a
                {
                    x = SafeProductWithExp(theta_a0,-s*t_e);
                    if (x < DBL_BIG)
                        result += (x - theta_a0) / (totalTheta_0*(theta_a0 + x*theta_A0/theta_a0));
                    else
                        result += theta_a0/(theta_A0*totalTheta_0);
                }
                else
                {
                    string msg = "CoalesceLogisticSelectionPL::DlnPoint(), received a TreeSummary ";
                    msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
                    msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescene in either ";
                    msg += "the population with allele A or the population with allele a.";
                    throw implementation_error(msg);
                }
            }
            else // derivative with respect to theta_a0
            {
                if (0L == pTreesum->oldstatus) // coal. in allele A
                {
                    x = SafeProductWithExp(theta_A0,-s*t_e);
                    if (x < DBL_BIG)
                        result -= (x - theta_A0) / (totalTheta_0*(theta_a0 + x));
                    else
                        result -= 1.0/totalTheta_0;
                }
                else if (1L == pTreesum->oldstatus) // coal. in allele a
                {
                    x = SafeProductWithExp(theta_A0*(theta_a0+totalTheta_0),-s*t_e);
                    if (x < DBL_BIG)
                        result -= (x + theta_a0squared) /
                            (theta_a0*totalTheta_0*(theta_a0 + x/(theta_a0+totalTheta_0)));
                    else
                        result -= (theta_a0+totalTheta_0)/(theta_a0*totalTheta_0);
                }
                else
                {
                    string msg = "CoalesceLogisticSelectionPL::DlnPoint(), received a TreeSummary ";
                    msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
                    msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescene in either ";
                    msg += "the population with allele A or the population with allele a.";
                    throw implementation_error(msg);
                }
            }
        }
    }

    // Now, add the contributions from the mutation events a --> A and A --> a.
    treesum = treedata->GetDiseaseSummary()->GetLongPoint();

    double one__over__theta_A0(1.0/theta_A0), one__over__theta_a0(1.0/theta_a0);
    for(pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
    {
        if (0L == pTreesum->oldstatus)
        {
            if (derivativeWithRespectTo_A)
                result -= one__over__theta_A0;
            else
                result += one__over__theta_a0;
        }
        else if (1L == pTreesum->oldstatus)
        {
            if (derivativeWithRespectTo_A)
                result += one__over__theta_A0;
            else
                result -= one__over__theta_a0;
        }
        else
        {
            string msg = "DiseaseLogisticSelectionPL::DlnPoint(), received a TreeSummary ";
            msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
            msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a mutation to either ";
            msg += "allele A or allele a.";
            throw implementation_error(msg);
        }
    }

    return result;

}

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

double LogisticSelectionPL::lnWait(const vector<double>& param, const TreeSummary *treedata)
{
    return 0.0;
}

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

double LogisticSelectionPL::lnPoint(const vector<double>& param, const vector<double>& lparam,
                                    const TreeSummary *treedata)
{
    return 0.0;
}

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

double LogisticSelectionPL::DlnWait(const vector<double>& param, const TreeSummary *treedata,
                                    const long int & whichparam)
{
    if (whichparam != m_s_is_here)
    {
        string msg = "LogisticSelectionPL::DlnWait() called with whichparam = ";
        msg += ToString(whichparam);
        msg += ".  \"whichparam\" must equal " + ToString(m_s_is_here);
        msg += ", which corresponds to the selection coefficient \"s.\"";
        throw implementation_error(msg);
    }

    double s(param[m_s_is_here]);

    if (fabs(s) < LOGISTIC_SELECTION_COEFFICIENT_EPSILON)
        return DlnWaitForTinyLogisticSelectionCoefficient(param, treedata, whichparam);

    double theta_A0(param[0]), theta_a0(param[1]), mu_into_A_from_a(param[3]),
        mu_into_a_from_A(param[4]);

    if (0.0 == theta_A0 || 0.0 == theta_a0)
    {
        string msg = "LogisticSelectionPL::DlnWait(), theta_A0 = ";
        msg += ToString(theta_A0) + ", theta_a0 = " + ToString(theta_a0);
        msg += "; attempted to divide by zero.";
        throw impossible_error(msg);
    }

    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;
    double t_s(0.0), t_e, term_s, term_e, denominator(s*s*(theta_A0+theta_a0)),
        dtau_A__ds(0.0), dtau_a__ds(0.0), result(0.0);
    double mu_term_A(mu_into_A_from_a*(theta_A0+theta_a0)), mu_term_a(mu_into_a_from_A*(theta_A0+theta_a0));
    term_s = s > 0.0 ? -theta_a0/denominator : theta_A0/denominator;

    for(treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        t_e = treesum_it->endtime;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a
        dtau_A__ds = dtau_a__ds = 0;

        if (s > 0.0)
        {
            term_e = SafeProductWithExp(theta_a0*(s*t_e - 1.0)/denominator,s*t_e);
            if (term_e >= DBL_BIG)
                return -DBL_BIG;
            if (k_A > 0)
                dtau_A__ds = term_e - term_s;
            if (k_a > 0)
                dtau_a__ds = ((theta_A0*theta_a0)/(denominator*denominator)) *
                    ((s*s*t_e*t_e - 1.0)/term_e - (s*s*t_s*t_s - 1.0)/term_s);
        }
        else
        {
            term_e = SafeProductWithExp(theta_A0*(s*t_e + 1.0)/denominator,-s*t_e);
            if (term_e <= -DBL_BIG)
                return DBL_BIG;
            if (k_a > 0)
                dtau_a__ds = term_e - term_s;
            if (k_A > 0)
                dtau_A__ds = ((theta_A0*theta_a0)/(denominator*denominator)) *
                    ((s*s*t_e*t_e - 1.0)/term_e - (s*s*t_s*t_s - 1.0)/term_s);
        }

        result -= (k_A*(k_A - 1.0 + mu_term_A)/theta_A0)*dtau_A__ds
            + (k_a*(k_a - 1.0 + mu_term_a)/theta_a0)*dtau_a__ds;
        if (result <= -DBL_BIG)
            return -DBL_BIG;
        if (result >= DBL_BIG)
            return DBL_BIG;

        t_s = t_e;
        term_s = term_e;
    }

    return result;

}

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

double LogisticSelectionPL::DlnWaitForTinyLogisticSelectionCoefficient(const vector<double> & param,
                                                                       const TreeSummary * treedata,
                                                                       const long int & whichparam)
{
    const list<Interval>& treesum = treedata->GetCoalSummary()->GetLongWait();
    list<Interval>::const_iterator treesum_it;
    double theta_A0(param[0]), theta_a0(param[1]), theta_A0_over_theta_a0;
    double ts(0.0), te, result(0.0), s(param[m_s_is_here]);
    double ts_squared(0.0), te_squared;

    if (0.0 == theta_A0 || 0.0 == theta_a0 || 0.0 == theta_A0 + theta_a0)
    {
        string msg = "LogisticSelectionPL::DlnWaitForTinySelectionCoefficient(), theta_A0 = ";
        msg += ToString(theta_A0) + ", theta_a0 = " + ToString(theta_a0);
        msg += "; attempted to divide by zero.";
        throw impossible_error(msg);
    }

    theta_A0_over_theta_a0 = theta_A0/theta_a0;
    double mu_into_A_from_a(param[3]), mu_into_a_from_A(param[4]),
        mu_term_A(mu_into_A_from_a*(theta_A0+theta_a0)), mu_term_a(mu_into_a_from_A*(theta_A0+theta_a0));

    for(treesum_it = treesum.begin(); treesum_it != treesum.end(); treesum_it++)
    {
        te = treesum_it->endtime;
        te_squared = te*te;
        const double& k_A = treesum_it->xpartlines[0]; // num. lineages, allele A
        const double& k_a = treesum_it->xpartlines[1]; // num. lineages, allele a

        result -= 0.5*(te_squared - ts_squared)*(k_A*(k_A-1.0+mu_term_A)/theta_A0_over_theta_a0 -
                                                 k_a*(k_a-1.0+mu_term_a)*theta_A0_over_theta_a0);
        result -= (s/3.0)*(te_squared*te - ts_squared*ts)*(k_A*(k_A-1.0+mu_term_A)/theta_A0_over_theta_a0 +
                                                           k_a*(k_a-1.0+mu_term_a)*theta_A0_over_theta_a0);
        ts= te;
        ts_squared = te_squared;
    }

    return result/(theta_A0 + theta_a0);
}

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

double LogisticSelectionPL::DlnPoint(const vector<double> & param, const TreeSummary * treedata,
                                     const long int & whichparam)
{
    if (whichparam != m_s_is_here)
    {
        string msg = "LogisticSelectionPL::DlnPoint() called with whichparam = ";
        msg += ToString(whichparam);
        msg += ".  \"whichparam\" must equal " + ToString(m_s_is_here);
        msg += ", which corresponds to the selection coefficient \"s.\"";
        throw implementation_error(msg);
    }

    const Interval *treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval *pTreesum;
    double theta_A0(param[0]), theta_a0(param[1]), s(param[m_s_is_here]);
    double result(0.0), term_a;

    // First, add up the contributions from coalescent events.

    for (pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
    {
        const double& te = pTreesum->endtime;
        term_a = SafeProductWithExp(theta_a0, s*te);

        if (0L == pTreesum->oldstatus) // coal. in allele A
        {
            if (term_a < DBL_BIG)
                result += te*term_a/(theta_A0 + term_a); // rhs = 0 for s << 0
            else
                result += te;
        }

        else if (1L == pTreesum->oldstatus) // coal. in allele a
            result -= te/(term_a/theta_A0 + 1.0);

        else
        {
            string msg = "LogisticSelectionPL::DlnPoint(), received a TreeSummary ";
            msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
            msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a coalescence in either ";
            msg += "the subpopulation with allele A or the subpopulation with allele a.";
            throw implementation_error(msg);
        }
    }

    // Now, add the contributions from the mutation events a --> A and A --> a.
    treesum = treedata->GetDiseaseSummary()->GetLongPoint();

    for(pTreesum = treesum; pTreesum != NULL; pTreesum = pTreesum->next)
    {
        const double& te = pTreesum->endtime;

        if (0L == pTreesum->oldstatus)
            result += te;
        else if (1L == pTreesum->oldstatus)
            result -= te;
        else
        {
            string msg = "LogisticSelectionPL::DlnPoint(), received a TreeSummary ";
            msg += "containing an event with oldstatus = " + ToString(pTreesum->oldstatus);
            msg += ".  \"oldstatus\" must be either 0 or 1, reflecting a mutation to either ";
            msg += "allele A or allele a.";
            throw implementation_error(msg);
        }
    }

    return result;
}

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

StickSelectPL::StickSelectPL(const ForceSummary& fs)
    : PLForces(0),  // MDEBUG arbitrary value, is that okay?
      m_thetastart(0),
      m_thetaend(0),
      m_toBigA_here(FLAGLONG),
      m_toSmallA_here(FLAGLONG),
      m_s_here(0),
      m_r_here(FLAGLONG),
      m_disstart(FLAGLONG)
{
    m_minuslnsqrt2pi = -log(sqrt(2*PI));
    m_nxparts = fs.GetNParameters(force_COAL);
    // JDEBUG--when switch to allowing migration remove the assert
    assert(fs.GetNonLocalPartitionIndexes().size() < 2);
    m_partindex = fs.GetNonLocalPartitionIndexes();
    // JDEBUG--this doesn't work due to invalid parameters
    // m_ndis = fs.GetNParameters(force_DISEASE);
    m_ndis = 2L;
    m_disindex = fs.GetPartIndex(force_DISEASE);
    assert(m_ndis == 2);
    assert(m_nxparts == 2);
}

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

StickSelectPL::~StickSelectPL()
{
    // intentionally blank
}

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

void StickSelectPL::SetToBigAIndex(long index)
{
    m_toBigA_here = index;
    if (m_disstart == FLAGLONG)
        // if toBig or toSmall isn't set, m_disstart will end
        // up as FLAGLONG again, as that's smaller than all legal values
        m_disstart = std::min(m_toBigA_here,m_toSmallA_here);
}

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

void StickSelectPL::SetToSmallAIndex(long index)
{
    m_toSmallA_here = index;
    if (m_disstart == FLAGLONG)
        // if toBig or toSmall isn't set, m_disstart will end
        // up as FLAGLONG again, as that's smaller than all legal values
        m_disstart = std::min(m_toBigA_here,m_toSmallA_here);
}

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

void StickSelectPL::SetLocalPartForces(const ForceSummary& fs)
{
    // JDEBUG--:  This code is cut and paste from CoalescePL
    // and smells pretty bad as a result.

    m_localpartforces = fs.GetLocalPartitionForces();
    assert(!m_localpartforces.empty());

    // set up vector which answers the following question:
    // for each local partition force
    //   for any xpartition theta (identified by index),
    //   what partition (of the given force) is it?

    // This code is copied approximately from PartitionForce::
    // SumXPartsToParts().

    ForceVec::const_iterator pforce;
    const ForceVec& partforces = fs.GetPartitionForces();
    LongVec1d indicator(partforces.size(), 0L);
    LongVec1d nparts(registry.GetDataPack().GetAllNPartitions());
    DoubleVec1d::size_type xpart;
    long partindex;
    for (pforce = partforces.begin(), partindex = 0;
         pforce != partforces.end(); ++pforce, ++partindex)
    {
        vector<DoubleVec1d::size_type> indices(m_nxparts, 0);
        for(xpart = 0; xpart < static_cast<DoubleVec1d::size_type>(m_nxparts); ++xpart)
        {
            indices[xpart] = indicator[partindex];
            long int part;
            for (part = nparts.size() - 1; part >= 0; --part)
            {
                ++indicator[part];
                if (indicator[part] < nparts[part]) break;
                indicator[part] = 0;
            }
        }
        // initialize xparts vectors
        vector<DoubleVec1d::size_type> emptyvec;
        vector<vector<DoubleVec1d::size_type> > partvec(nparts[partindex], emptyvec);

        if ((*pforce)->IsLocalPartitionForce())
        {
            m_whichlocalpart.push_back(indices);
            m_whichlocalxparts.push_back(partvec);
        }
        m_whichpart.push_back(indices);
        m_whichxparts.push_back(partvec);
    }

    // now construct the vector mapping partition to a set of xpartitions.
    //
    // we will do this by going through the vector mapping xpartition to
    // partition we just constructed, letting its contents tell us which
    // partition to add the parameter (xpartition) to and keeping a
    // counter to let us know which xpartition to add.
    long lpindex = 0;
    for (pforce = partforces.begin(), partindex = 0;
         pforce != partforces.end(); ++pforce, ++partindex)
    {
        for(xpart = 0; xpart < m_whichpart[partindex].size(); ++xpart)
            m_whichxparts[partindex][m_whichpart[partindex][xpart]].
                push_back(xpart);
        if ((*pforce)->IsLocalPartitionForce())
        {
            for(xpart = 0; xpart < m_whichpart[partindex].size(); ++xpart)
            {
                m_whichlocalxparts[lpindex][m_whichpart[partindex][xpart]].push_back(xpart);
            }
            ++lpindex;
        }
    }
} // SetLocalPartForces

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

double StickSelectPL::StickMean(double freqA, double s, double toA, double toa) const
{
    return s*freqA*(1-freqA) + toA*(1-freqA) - toa*freqA;
}

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

double StickSelectPL::StickVar(double freqA, double tipfreqA, double thetaA) const
{
    return 2.0*freqA*(1-freqA)*tipfreqA / thetaA;
}

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

double StickSelectPL::CommonFunctional(double freqBigA, double prevfreqBigA, double s,
                                       double toBigA, double toSmallA, double length) const
{
    // this implements the un-squared numerator of the posterior on stick selection:
    return freqBigA - prevfreqBigA - length *
        StickMean(prevfreqBigA,s,toBigA,toSmallA);
}

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

double StickSelectPL::lnPTreeStick(const DoubleVec1d& param, const DoubleVec1d& lparam, const TreeSummary* treedata)
{
    double answ(lnWait(param,treedata));

    answ += lnPointTreeStick(param,lparam,treedata);

    return answ;
}

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

double StickSelectPL::lnWait(const DoubleVec1d& param, const TreeSummary* treedata)
{
    assert(m_disstart != FLAGLONG);
    double answ(0.0);

    assert(m_ndis == 2);
    assert(m_ndis == m_nxparts);

    // compute waiting times for theta and disease
    // all LongWait() are the same....
    const list<Interval>& coalsum(treedata->GetCoalSummary()->GetLongWait());
    list<Interval>::const_iterator interval(coalsum.begin());
    const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
    const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
    double tipthetaA(param[m_thetastart]);
    double jointend(stairlengths[0]),intervalend(interval->endtime);
    double jointstart(0.0),intervalstart(0.0);
    long joint(0);

    while (interval != coalsum.end())
    {
        // how long is the bit of tree where nothing changes?
        double lengthininterval = std::min(jointend, intervalend) -
            std::max(jointstart, intervalstart);

        // adjust for coalescent contribution
        long xpart;
        for(xpart = 0; xpart < m_nxparts; ++xpart)
        {
            answ -= lengthininterval * stairfreqs[0][0] *
                (interval->xpartlines[xpart])*(interval->xpartlines[xpart]-1) /
                // because there is no migration we can use xpart, otherwise
                // we'll need to add a loop in local partitions
                (stairfreqs[joint][xpart] * tipthetaA);
        }

        // adjust for trait mutation contribution
        long part;
        for(part = 0; part < m_ndis; ++part)
        {
            answ -= lengthininterval * interval->partlines[m_disindex][part] *
                param[(part ? m_disstart+1 : m_disstart)] *
                // both macros rely on only two disease states present
                // equation relies on nxparts == nparts!
                stairfreqs[joint][(part ? 0 : 1)]/stairfreqs[joint][part];
        }

        if (jointend < intervalend)     // joint ends first
        {
            ++joint;
            jointstart = jointend;
            jointend += stairlengths[joint];
        }
        else
            if (intervalend < jointend) // else interval ends first
            {
                ++interval;
                intervalstart = intervalend;
                intervalend = interval->endtime;
            }
            else                          // they both ended at the exact same time!
            {
                ++joint;
                jointstart = jointend;
                jointend += stairlengths[joint];
                ++interval;
                intervalstart = intervalend;
                intervalend = interval->endtime;
            }
    }

    // compute rec-rate contrib
    if (!m_localpartforces.empty())
    {
        const vector<double>& rWait = treedata->GetRecSummary()->GetShortWait();
        answ -= rWait[0] * param[m_r_here];
    }

    // there aren't separate "wait" and "point" terms wrt. sticks,
    // therefore all the actual calculations are handled in "point".

    return answ;
}

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

double StickSelectPL::lnPointTreeStick(const DoubleVec1d& param, const DoubleVec1d& lparam,
                                       const TreeSummary* treedata)
// there aren't separate "wait" and "point" terms with stick selection,
// all the actual calculations are handled in "point".
//
// The actual equation for the log tree posterior likelihood:
// Sum_Over_All_Stick_Intervals[
//    -log(sqrt(2*pi))-log(sqrt(Var(curr_freq)*int_length))-Exponent]
//
// Exponent =
//    (curr_freq-(prev_freq*Mean(prev_freq)*int_length))**2
//    ------------------------------------------------------------
//             2*Var(prev_freq)*int_length
//
// Mean(freq) = sel_coeff*freq*(1-freq)+mu_to*(1-freq)-mu_from*freq
// Var(freq) = freq*(1-freq) / theta
{
    assert(m_disstart != FLAGLONG);
    double answ(0.0);

    // stuff useful for all stick terms
    const Interval * treesum = treedata->GetCoalSummary()->GetLongPoint();
    const Interval * tit(treesum);
    const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
    const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
    const DoubleVec2d& stairlnfreqs(treedata->GetStickSummary().lnfreqs);
    double tipprobA(stairfreqs[0][0]);
    double lnthA(lparam[m_thetastart]), lnprA(log(tipprobA));
    double jointend(stairlengths[0]);
    long joint(0);

    // compute theta contrib
    do {
        while (jointend < tit->endtime)
        {
            ++joint;
            jointend += stairlengths[joint];
        }

        answ +=  LOG2 + lnprA - lnthA - stairlnfreqs[joint][tit->oldstatus];

        tit = tit->next;
    } while(tit != NULL);

    // compute trait mutation contrib
    treesum = treedata->GetDiseaseSummary()->GetLongPoint();
    tit = treesum;
    jointend = stairlengths[0];
    joint = 0;
    do {
        while (jointend < tit->endtime)
        {
            ++joint;
            jointend += stairlengths[joint];
        }
        // JDEBUG--is the macro correct??
        answ +=  lparam[(tit->newstatus ? m_disstart+1 : m_disstart)] +
            stairlnfreqs[joint][tit->oldstatus]-stairlnfreqs[joint][tit->newstatus];

        tit = tit->next;
    } while(tit != NULL);

    // compute rec-rate contrib--Sum_j(recevents * log(r))
    // as noted earlier (see Likelihood::FillForces()), the existence
    // of a localpartforces vector is used to flag the presence/absence
    // of recombination

    // MDEBUG MREVIEW let's review this code--Mary doesn't get it.
    if (!m_localpartforces.empty())
    {
        treesum = treedata->GetRecSummary()->GetLongPoint();
        tit = treesum;
        jointend = stairlengths[0];
        joint = 0;
        do {
            while (jointend < tit->endtime)
            {
                ++joint;
                jointend += stairlengths[joint];
            }
            // we can use partnerpicks[0] because there is only one lpforce
            // in existence
            answ += lparam[m_r_here] +
                stairlnfreqs[joint][treesum->partnerpicks[0]];
            tit = tit->next;
        } while(tit != NULL);
    }

    return answ;
}

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

double StickSelectPL::lnPoint(const DoubleVec1d& param, const DoubleVec1d& lparam, const TreeSummary* treedata)
// there aren't separate "wait" and "point" terms with stick selection,
// all the actual calculations are handled in "point".
//
// The actual equation for the log tree posterior likelihood:
// Sum_Over_All_Stick_Intervals[
//    -log(sqrt(2*pi))-log(sqrt(Var(prev_freq)*int_length))-Exponent]
//
// Exponent =
//    (curr_freq-(prev_freq-Mean(prev_freq)*int_length))**2
//    ------------------------------------------------------------
//             2*Var(prev_freq)*int_length
//
// Mean(freq) = sel_coeff*freq*(1-freq)+mu_to*(1-freq)-mu_from*freq
// Var(freq) = freq*(1-freq) / theta_total
{
    assert(m_disstart != FLAGLONG);
    double answ(lnPointTreeStick(param,lparam,treedata));

    const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
    const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
    double tipthetaA(param[m_thetastart]);
    double tipprobA(stairfreqs[0][0]);
    double s(param[m_s_here]),toBigA(param[m_toBigA_here]),
        toSmallA(param[m_toSmallA_here]);
    assert(stairfreqs.size() == stairlengths.size());

    DoubleVec1d::size_type step;
    for(step = 1; step < stairfreqs.size(); ++step)
    {
        double variance = StickVar(stairfreqs[step-1][0],tipprobA,tipthetaA);
        answ += m_minuslnsqrt2pi - 0.5 * log(variance*stairlengths[step-1]);
        double numer = CommonFunctional(stairfreqs[step][0],
                                        stairfreqs[step-1][0],s,toBigA,toSmallA,stairlengths[step-1]);
        double denom(2.0*variance*stairlengths[step-1]);
        numer = -1.0 * numer * numer;
        answ += numer/denom;
    }

    return answ;
}

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

double StickSelectPL::DlnWait(const DoubleVec1d& param, const TreeSummary* treedata, const long& whichparam)
{
    long which;
    double answ(0.0);

    // compute theta contrib
    if (m_thetastart <= whichparam && whichparam < m_thetastart + m_nxparts)
    {
        const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
        const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
        const list<Interval>& treesum(treedata->GetCoalSummary()->GetLongWait());
        list<Interval>::const_iterator interval(treesum.begin());
        double tipthetaA(param[m_thetastart]);
        double tipprobA(stairfreqs[0][0]);
        double jointend(stairlengths[0]),intervalend(interval->endtime);
        double jointstart(0.0),intervalstart(0.0);
        long joint(0);
        which = whichparam-m_thetastart;
        while (interval != treesum.end())
        {
            // how long is the bit of tree where nothing changes?
            double lengthininterval = std::min(jointend, intervalend) -
                std::max(jointstart, intervalstart);

            answ += lengthininterval * tipprobA *
                (interval->xpartlines[which])*(interval->xpartlines[which]-1)
                // because there is no migration we can use xpart, otherwise
                // we'll need to add a loop in local partitions
                / (tipthetaA * tipthetaA * stairfreqs[joint][which]);

            if (jointend < intervalend) // joint ends first
            {
                ++joint;
                jointstart = jointend;
                jointend += stairlengths[joint];
            }
            else
                if (intervalend < jointend) // else interval ends first
                {
                    ++interval;
                    intervalstart = intervalend;
                    intervalend = interval->endtime;
                }
                else                      // they both ended at the exact same time!
                {
                    ++joint;
                    jointstart = jointend;
                    jointend += stairlengths[joint];
                    ++interval;
                    intervalstart = intervalend;
                    intervalend = interval->endtime;
                }
        }

        return answ;
    }

    // compute disease contrib
    if (m_disstart <= whichparam && whichparam < m_disstart + m_ndis)
    {
        assert(m_disstart != FLAGLONG);
        if (whichparam == m_disstart) which = 0;
        else which = 1;
        const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
        const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
        const list<Interval>& treesum(treedata->GetDiseaseSummary()->GetLongWait());
        list<Interval>::const_iterator interval(treesum.begin());
        double jointend(stairlengths[0]),intervalend(interval->endtime);
        double jointstart(0.0),intervalstart(0.0);
        long joint(0);

        while (interval != treesum.end())
        {
            // how long is the bit of tree where nothing changes?
            double lengthininterval = std::min(jointend, intervalend) -
                std::max(jointstart, intervalstart);

            answ -= lengthininterval * interval->partlines[m_disindex][which] *
                // the macro relies on only two disease states present
                // equation relies on nxparts == nparts!
                stairfreqs[joint][(which ? 0 : 1)]/stairfreqs[joint][which];

            if (jointend < intervalend) // joint ends first
            {
                ++joint;
                jointstart = jointend;
                jointend += stairlengths[joint];
            }
            else
                if (intervalend < jointend) // else interval ends first
                {
                    ++interval;
                    intervalstart = intervalend;
                    intervalend = interval->endtime;
                }
                else                      // they both ended at the exact same time!
                {
                    ++joint;
                    jointstart = jointend;
                    jointend += stairlengths[joint];
                    ++interval;
                    intervalstart = intervalend;
                    intervalend = interval->endtime;
                }
        }

        return answ;
    }

    // compute recombination contrib
    if (whichparam == m_r_here)
    {
        assert(!m_localpartforces.empty());
        const vector<double>& rWait = treedata->GetRecSummary()->GetShortWait();
        return -rWait[0];
    }

    // We compute all terms involving the stick in lnPoint and DlnPoint
    // (arbitrarily) and thus return 0 if asked to take dWait of those parameters
    return answ;
}

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

double StickSelectPL::DlnPoint(const DoubleVec1d& param, const TreeSummary* treedata, const long& whichparam)
{
    // NB Assumes only one localpartforce and no migration, no growth!
    long which;

    // compute recombination contrib (does not involve stick)
    if (whichparam == m_r_here)
    {
        const vector<double>& nrecs = treedata->GetRecSummary()->GetShortPoint();
        if(nrecs[0]==0)
        {
            return 0.0;
        }
        else
        {
            return SafeDivide(nrecs[0], *(param.begin() + m_r_here));
        }
    }

    // useful variables for all stick computations
    const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
    const DoubleVec1d& stairlengths(treedata->GetStickSummary().lengths);
    double s(param[m_s_here]),toBigA(param[m_toBigA_here]),
        toSmallA(param[m_toSmallA_here]);
    double tipthetaA(param[m_thetastart]);
    double tipprobA(stairfreqs[0][0]);
    assert(stairfreqs.size() == stairlengths.size());
    double answ(0.0);

    // compute theta contrib
    if (m_thetastart <= whichparam && whichparam < m_thetastart + m_nxparts)
    {
        which = whichparam-m_thetastart;
        const DoubleVec1d& ncoal = treedata->GetCoalSummary()->GetShortPoint();

        answ = -ncoal[which] / tipthetaA;

        // deal with the stick contribution to theta
        DoubleVec1d::size_type step;
        for(step = 1; step < stairfreqs.size(); ++step)
        {
            double numer = CommonFunctional(stairfreqs[step][0],
                                            stairfreqs[step-1][0],s,toBigA,toSmallA,stairlengths[step-1]);
            numer *= numer;
            answ -= numer / (4.0 * stairfreqs[step-1][0] * tipprobA *
                             (1.0 - stairfreqs[step-1][0]) * stairlengths[step-1]);
            answ += 1.0 / (2 * tipthetaA);
        }

        return answ;
    }

    // compute disease contrib
    if (m_disstart <= whichparam && whichparam < m_disstart + m_ndis)
    {
        // non-stick terms
        assert(m_disstart != FLAGLONG);
        if (whichparam == m_disstart) which = 0;
        else which = 1;

        answ += (treedata->GetDiseaseSummary()->GetShortPoint()[which]) / param[whichparam];

        // stick terms
        const DoubleVec2d& stairfreqs(treedata->GetStickSummary().freqs);
        DoubleVec1d::size_type step;
        // terms in toBigA
        if (whichparam == m_toBigA_here)
        {
            for(step = 1; step < stairfreqs.size(); ++step)
            {
                answ += CommonFunctional(stairfreqs[step][0],stairfreqs[step-1][0],s,
                                         toBigA,toSmallA,stairlengths[step-1]) * tipthetaA /
                    (2.0 * tipprobA * stairfreqs[step-1][0]);
            }
            return answ;
        }
        // terms in toSmallA
        if (whichparam == m_toSmallA_here)
        {
            for(step = 1; step < stairfreqs.size(); ++step)
            {
                answ -= CommonFunctional(stairfreqs[step][0],stairfreqs[step-1][0],s,
                                         toBigA,toSmallA,stairlengths[step-1]) * tipthetaA /
                    ((2.0 * tipprobA) * (1.0 - stairfreqs[step-1][0]));
            }
            return answ;
        }
        assert(false); // neither toBigA nor toSmallA?!
    }

    // compute s contrib
    if (whichparam == m_s_here)
    {
        DoubleVec1d::size_type step;
        for(step = 1; step < stairfreqs.size(); ++step)
        {
            answ += CommonFunctional(stairfreqs[step][0],stairfreqs[step-1][0],
                                     s,toBigA,toSmallA,stairlengths[step-1]) * tipthetaA /
                (2.0 * tipprobA);
        }
        return answ;
    }

    assert(false); // should have computed something!
    return answ;
}

//____________________________________________________________________________________