src/datalike/dlmodel.cpp

// $Id: dlmodel.cpp,v 1.137 2012/05/22 20:25:32 jyamato Exp $

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

  This software is distributed free of charge for non-commercial use
  and is copyrighted.  Of course, we do not guarantee that the software
  works, and are not responsible for any damage you may cause or have.
*/

#include <cassert>
#include <algorithm>
#include <cstring>
#include <numeric>
#include <iostream>

#include "local_build.h"

#include "calculators.h"
#include "datapack.h"
#include "datatype.h"
#include "defaults.h"                   // for defaults::threshhold in StepwiseModel::ctor
#include "dlcell.h"
#include "dlmodel.h"
#include "errhandling.h"
#include "funcMax.h"
#include "locus.h"
#include "mathx.h"
#include "registry.h"
#include "runreport.h"
#include "stringx.h"
#include "xml_strings.h"                // for ToXML()
#include "xmlsum_strings.h"             // For WriteAlpha()

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

// turns on detailed local variables for debugging
// JRM 3/10
//#define DEBUG_VARIABLES
#ifdef DEBUG_VARIABLES
    int printdbgvar = 0;
#endif

using namespace std;

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

DataModel::DataModel(long nmarkers,
                     long numCategories,
                     DoubleVec1d categoryRates,
                     DoubleVec1d categoryProbabilities,
                     double userAutoCorrelationValue,
                     bool doNormalize,
                     long numBins,
                     double relmurate)
    :
    m_ncategories(numCategories),
    m_catrates(categoryRates),
    m_catprobs(categoryProbabilities),
    m_acratio(1.0 / userAutoCorrelationValue),
    m_normalize(doNormalize),
    m_nbins(numBins),
    m_nmarkers(nmarkers),
    m_relmurate(relmurate)
{
    m_ucratio = 1.0 - m_acratio;
    ScaleCatProbabilities();
    ScaleCatRates();
    assert(DataModel::IsValidDataModel());

} // DataModel constructor

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

string DataModel::GetDataModelName() const
{
    // "true" argument gets long version of name
    return ToString(GetModelType(),true);
}

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

string DataModel::GetDataModelShortName() const
{
    // "false" argument gets short version of name
    return ToString(GetModelType(),false);
}

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

DoubleVec1d DataModel::ChooseRandomRates(long nsites) const
{
    DoubleVec1d rates;
    Random& rand = registry.GetRandom();
    double rate = m_catrates[0];
    //choose a random rate to start
    double r = rand.Float();
    for (long i = 0; i < m_ncategories; ++i)
    {
        if (r < m_catprobs[i]) rate = m_catrates[i];
        else r -= m_catprobs[i];
    }

    for (long site=0; site<nsites; site++)
    {
        if (rand.Float() <= m_acratio)
        {
            //choose a new rate.
            r = rand.Float();
            for (long i = 0; i < m_ncategories; ++i)
            {
                if (r < m_catprobs[i]) rate = m_catrates[i];
                else r -= m_catprobs[i];
            }
        }
        rates.push_back(rate);
    }
    return rates;

} // ChooseRandomRate

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

#if 0 // Vestigial but possibly future code
void DataModel::SetGamma(bool gam)
{
    usegamma = gam;
    if (usegamma)
    {
        // DEBUG debug warning WARNING
        // need to error check presence of gammashape and ncats > 1?
        SetCatRates(gamma_rates(gammashape,m_ncategories));
        DoubleVec1d newprobs(m_ncategories,1.0/m_ncategories);
        SetCatProbabilities(newprobs);
    }
}
#endif

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

bool DataModel::IsValidDataModel() const
{
    if (m_nbins < 1) return false;
    if (m_nmarkers < 1) return false;

    if (m_acratio < 0) return false;

    if (m_ncategories < 1) return false;
    size_t ncats = m_ncategories;
    if (m_catrates.size() != ncats) return false;
    if (m_catprobs.size() != ncats) return false;

    size_t i;
    double totalprob = 0.0;
    for (i = 0; i < ncats; ++i)
    {
        if (m_catrates[i] < 0.0) return false;
        if (m_catprobs[i] <= 0.0) return false;
        if (m_catprobs[i] > 1.0) return false;
        totalprob += m_catprobs[i];
    }

    if (fabs(1.0 - totalprob) > EPSILON) return false;

    return true;
} // DataModel::IsValidDataModel

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

StringVec1d DataModel::CreateDataModelReport() const
{

    StringVec1d report;
    string line;

    if (m_ncategories > 1)
    {
        line = ToString(m_ncategories) + " rate categories with correlated length " +
            ToString(1.0 / m_acratio);
        report.push_back(line);
        long cat;
        for (cat = 0; cat < m_ncategories; ++cat)
        {
            line = "Relative rate " + ToString(m_catrates[cat]) + "  Frequency " +
                ToString(m_catprobs[cat]);
            report.push_back(line);
        }
    }
    if (m_relmurate != 1)
    {
        line = "The relative marker mutation rate for this model's segment was ";
        line += ToString(m_relmurate) + ".";
        report.push_back(line);
    }

    return report;
} // CreateDataModelReport

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

StringVec1d DataModel::ToXML(size_t nspaces) const
{
    StringVec1d xmllines;
    string line = MakeIndent(
        MakeTagWithName(xmlstr::XML_TAG_MODEL,GetDataModelShortName()),
        nspaces);
    xmllines.push_back(line);

    nspaces += INDENT_DEPTH;
    string mytag(MakeTag(xmlstr::XML_TAG_NORMALIZE));
    line = MakeIndent(mytag,nspaces) + ToStringTF(ShouldNormalize()) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);
    line = MakeIndent(MakeTag(xmlstr::XML_TAG_CATEGORIES),nspaces);
    xmllines.push_back(line);

    nspaces += INDENT_DEPTH;
    mytag = MakeTag(xmlstr::XML_TAG_NUM_CATEGORIES);
    line = MakeIndent(mytag,nspaces) + ToString(GetNcategories()) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);
    mytag = MakeTag(xmlstr::XML_TAG_RATES);
    line = MakeIndent(mytag,nspaces) + ToString(GetCatRates(),6) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);
    mytag = MakeTag(xmlstr::XML_TAG_PROBABILITIES);
    line = MakeIndent(mytag,nspaces) + ToString(GetCatProbabilities()) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);
    mytag = MakeTag(xmlstr::XML_TAG_AUTOCORRELATION);
    line = MakeIndent(mytag,nspaces) + ToString(GetUserAcratio()) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);
    nspaces -= INDENT_DEPTH;

    line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_CATEGORIES),nspaces);
    xmllines.push_back(line);

    mytag = MakeTag(xmlstr::XML_TAG_RELATIVE_MURATE);
    line = MakeIndent(mytag,nspaces) + ToString(m_relmurate) +
        MakeCloseTag(mytag);
    xmllines.push_back(line);

    nspaces -= INDENT_DEPTH;
    line = MakeIndent(MakeCloseTag(xmlstr::XML_TAG_MODEL),nspaces);
    xmllines.push_back(line);

    return xmllines;
} // ToXML

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

void DataModel::ResetCatCells()
{
    m_catcells.assign(m_ncategories,1.0);
} // DataModel::ResetCatCells

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

void DataModel::ScaleCatProbabilities()
{
    long cat;
    double totalprob = 0.0;
    m_ncategories = m_catprobs.size();
    for(cat = 0; cat < m_ncategories; cat++)
    {
        totalprob += m_catprobs[cat];
    }
    if (fabs(1.0 - totalprob) > EPSILON)
    {
        for(cat = 0; cat < m_ncategories; cat++)
        {
            m_catprobs[cat] = m_catprobs[cat] / totalprob;
        }
    }
} // DataModel::ScaleCatProbabilities

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

void DataModel::ScaleCatRates()
{
    // We renormalize the category rates to a weighted mean of 1.0
    double meanrate = 0.0;
    size_t rate;

    if (m_catrates.size() != m_catprobs.size())
    {
        string msg = "DataModel::ScaleCatRates() was called before ";
        msg += "m_catrates and m_catprobs were properly initialized.";
        throw implementation_error(msg);
    }

    size_t numrates = m_catrates.size();

    for (rate = 0; rate < numrates; ++rate)
    {
        meanrate += m_catrates[rate] * m_catprobs[rate];
    }

    for (rate = 0; rate < numrates; ++rate)
    {
        m_catrates[rate] /= meanrate;
    }

} // DataModel::ScaleCatRates

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

void DataModel::TryToNormalizeAndThrow(long posn, model_type mtype)
{
    if (ShouldNormalize())
    {
        string errmsg("Encountered a subtree of zero likelihood at ");
        errmsg += "position " + ToString(posn) + ".";
        switch (mtype)
        {
            case F84:
            case GTR:
            case KAllele:
                assert(false);
                //Encountering a tree of zero data likelihood on a tree with DNA
                // is almost certainly a programming error, though occasionally
                // one gets a tree so large that mispairing tips can be
                // catastrophic.  Likewise, a K-allele model should allow mutations
                // of a single step to have reasonable likelihoods even if the
                // microsats are very disparate.
                //
                // However, a *user* is much more likely to have such a data set,
                // (and hopefully we debug the program thoroughly before release)
                // so we fall through here, throw the 'bad tree' error, and try
                // a different tree.  --Lucian
            case MixedKS:
                // The MixedKS model has been seen to have problems here when it
                // is optimizing.
            case Brownian:
            case Stepwise:
                throw zero_dl_error(errmsg);
        }
    }
    else
    {
        SetNormalize(true);
        throw datalikenorm_error("Datalikelihood normalization turned on.");
    }
}

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

NucModel::NucModel(long nmarkers,
                   long numCategories,
                   DoubleVec1d categoryRates,
                   DoubleVec1d categoryProbabilities,
                   double userAutoCorrelationValue,
                   bool doNormalize,
                   double relMuRate,
                   double freqA,
                   double freqC,
                   double freqG,
                   double freqT,
                   bool calcFreqsFromData,
                   double perBaseErrorRate)
    : DataModel(nmarkers,
                numCategories,
                categoryRates,
                categoryProbabilities,
                userAutoCorrelationValue,
                doNormalize,
                defaults::nucleotideBins,
                relMuRate),
      m_basefreqs(defaults::nucleotideBins),
      m_freqsfromdata(calcFreqsFromData),
      m_perBaseErrorRate(perBaseErrorRate)
{
    m_basefreqs[baseA] = freqA;
    m_basefreqs[baseC] = freqC;
    m_basefreqs[baseG] = freqG;
    m_basefreqs[baseT] = freqT;
    NormalizeBaseFrequencies();
    assert(NucModel::IsValidDataModel());
} // NucModel constructor

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

void NucModel::NormalizeBaseFrequencies()
{
    double sumOfFreqs = accumulate(m_basefreqs.begin(),m_basefreqs.end(),0.0);
    transform(m_basefreqs.begin(), m_basefreqs.end(),
              m_basefreqs.begin(),
              bind2nd(divides<double>(),sumOfFreqs));
}

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

bool NucModel::IsValidDataModel() const
{
    size_t index;
    double totalfreqs = 0.0;
    if(m_basefreqs.size() != 4) return false;
    for (index = 0; index < m_basefreqs.size(); index++)
    {
        double thisFreq = m_basefreqs[index];
        if(thisFreq <= 0.0) return false;
        if(thisFreq >= 1.0) return false;
        totalfreqs += thisFreq;
    }
    if (fabs(1.0 - totalfreqs) > EPSILON) return false;
    return DataModel::IsValidDataModel();
} // NucModel::ValidBaseFrequencies

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

vector<double> NucModel::DataToLikes(const string& datum, long) const
{
    return NucModel::StaticDataToLikes(datum, GetPerBaseErrorRate());
}

//------------------------------------------------------------------------------------
// This function implements the standard ambiguity codes for nucleotide data,
// returning a vector of four doubles indicating the likelihood for A, C, G, T in that order.

vector<double> NucModel::StaticDataToLikes(const string& datum, double perBaseErrorRate)
{
    // We assume this is a single nucleotide base, passed as a string only for generality
    assert(datum.size() == 1);

    vector<double> likes(BASES, 0.0);  // initialize to zero
    char nucleotide = datum[0];

    // resolve code
    switch(nucleotide)
    {
        case 'A':
            likes[baseA] = 1.0;
            break;

        case 'C':
            likes[baseC] = 1.0;
            break;

        case 'G':
            likes[baseG] = 1.0;
            break;

        case 'T':
        case 'U':
            likes[baseT] = 1.0;
        break;

        case 'M':
            likes[baseA] = 1.0;
            likes[baseC] = 1.0;
            break;

        case 'R':
            likes[baseA] = 1.0;
            likes[baseG] = 1.0;
            break;

        case 'W':
            likes[baseA] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'S':
            likes[baseC] = 1.0;
            likes[baseG] = 1.0;
            break;

        case 'Y':
            likes[baseC] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'K':
            likes[baseG] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'V':
            likes[baseA] = 1.0;
            likes[baseC] = 1.0;
            likes[baseG] = 1.0;
            break;

        case 'H':
            likes[baseA] = 1.0;
            likes[baseC] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'D':
            likes[baseA] = 1.0;
            likes[baseG] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'B':
            likes[baseC] = 1.0;
            likes[baseG] = 1.0;
            likes[baseT] = 1.0;
            break;

        case 'N':
        case 'O':
        case 'X':
        case '?':
        case '-':
            likes[baseA] = 1.0;
        likes[baseC] = 1.0;
        likes[baseG] = 1.0;
        likes[baseT] = 1.0;
        break;

        default:
            assert(false);    // how did an unknown nucleotide get past proofreading?
            likes[baseA] = 1.0;
            likes[baseC] = 1.0;
            likes[baseG] = 1.0;
            likes[baseT] = 1.0;
            break;
    }

    long   num_ones  = (long)(likes[baseA] + likes[baseC] + likes[baseG] + likes[baseT]);
    double new_0     = (double)num_ones * perBaseErrorRate / 3.0;
    double new_1     = 1.0 - (double)(4-num_ones) * perBaseErrorRate / 3.0;
    likes[baseA] = (likes[baseA] > 0.5) ? new_1 : new_0;
    likes[baseC] = (likes[baseC] > 0.5) ? new_1 : new_0;
    likes[baseG] = (likes[baseG] > 0.5) ? new_1 : new_0;
    likes[baseT] = (likes[baseT] > 0.5) ? new_1 : new_0;

    return(likes);

} // DataToLikes

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

StringVec1d NucModel::CreateDataModelReport() const
{
    StringVec1d report = DataModel::CreateDataModelReport();

    string line = "Base frequencies: " + ToString(m_basefreqs[baseA]) + ", ";
    line += ToString(m_basefreqs[baseC]) + ", ";
    line += ToString(m_basefreqs[baseG]) + ", ";
    line += ToString(m_basefreqs[baseT]);
    report.push_back(line);

    return report;

} // NucModel::CreateDataModelReport

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

StringVec1d NucModel::ToXML(size_t nspaces) const
{
    StringVec1d xmllines(DataModel::ToXML(nspaces));

    nspaces += INDENT_DEPTH;

    string line(MakeIndent(MakeTag(xmlstr::XML_TAG_BASE_FREQS),nspaces));
    if (m_freqsfromdata)
        line += " " + xmlstr::XML_TAG_CALCULATED + " ";
    else
        line += ToString(m_basefreqs,6);
    line += MakeCloseTag(xmlstr::XML_TAG_BASE_FREQS);
    StringVec1d::iterator endtag = --xmllines.end();
    xmllines.insert(endtag,line);

    string line2(MakeIndent(MakeTag(xmlstr::XML_TAG_PER_BASE_ERROR_RATE),nspaces));
    line2 += ToString(m_perBaseErrorRate);
    line2 += MakeCloseTag(xmlstr::XML_TAG_PER_BASE_ERROR_RATE);
    endtag = --xmllines.end();
    xmllines.insert(endtag,line2);

    nspaces -= INDENT_DEPTH;
    return xmllines;

} // NucModel::ToXML

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

DoubleVec1d NucModel::ChooseAncestralState(long marker)
{
    DoubleVec1d result(BASES, 0.0);
    double r = registry.GetRandom().Float();
    long i;
    for (i = 0; i < BASES; ++i)
    {
        if (r < m_basefreqs[i])
        {
            result[i] = 1.0;
            return result;
        }
        else
        {
            r -= m_basefreqs[i];
        }
    }

    // this code could be reached due to rounding errors
    result[0] = 1.0;
    return result;
} // ChooseAncestralState

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

string NucModel::CellToData(Cell_ptr cell, long marker) const
{
    LongVec1d ones = cell->GetOnes(marker);
    if (ones.size() == 1)
    {
        switch (ones[0])
        {
            case 0:
                return "A";
            case 1:
                return "C";
            case 2:
                return "G";
            case 3:
                return "T";
            default:
                throw data_error("Tried to convert " + ToString(ones[0]+1)
                                 + " to a nucleotide, but there should only be four bins.");
        }
    }
    throw implementation_error("Cannot convert nucleotide data from the internal format"
                               " that is not simply one of the four bases.");
}

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

F84Model::F84Model(
    long nmarkers,
    long numCategories,
    DoubleVec1d categoryRates,
    DoubleVec1d categoryProbabilities,
    double userAutoCorrelationValue,
    bool doNormalize,
    double relMuRate,
    double freqA,
    double freqC,
    double freqG,
    double freqT,
    double ttRatio,
    bool calculateFreqsFromData,
    double perBaseErrorRate)
    : NucModel(
        nmarkers,
        numCategories,
        categoryRates,
        categoryProbabilities,
        userAutoCorrelationValue,
        doNormalize,
        relMuRate,
        freqA,
        freqC,
        freqG,
        freqT,
        calculateFreqsFromData,
        perBaseErrorRate),
      m_ttratio(ttRatio),
      // set some potentially useful defaults
      computed1(true), computed2(true),
      freqar(FLAGDOUBLE), freqcy(FLAGDOUBLE),
      freqgr(FLAGDOUBLE), freqty(FLAGDOUBLE),
      basefreqarray(NULL),
      daughter1(NULL), daughter2(NULL), target(NULL)
{
    // The data-likelihood buffer vectors such as expA1 are default-constructed
    // to an empty state, which is fine.
    assert(F84Model::IsValidDataModel());
    //cout << "*****new F84 model created*****" << this << endl;
    Finalize();

} // F84Model::F84Model

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

void F84Model::EmptyBuffers()
{
    if (basefreqarray)
    {
        delete [] basefreqarray[0];
        delete [] basefreqarray;
    }
    if (daughter1)
    {
        delete [] daughter1[0];
        delete [] daughter1;
    }
    if (daughter2)
    {
        delete [] daughter2[0];
        delete [] daughter2;
    }
    if (target)
    {
        delete [] target[0];
        delete [] target;
    }

} // F84Model::EmptyBuffers

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

DataModel* F84Model::Clone() const
{
    DataModel* newmodel = new F84Model(*this);
    return newmodel;
}

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

void F84Model::AllocateBuffers()
{
    basefreqarray   = new double* [m_ncategories];
    daughter1 = new double* [m_ncategories];
    daughter2 = new double* [m_ncategories];
    target = new double* [m_ncategories];
    basefreqarray[0]   = new double [m_ncategories*BASES];
    daughter1[0] = new double [m_ncategories*BASES];
    daughter2[0] = new double [m_ncategories*BASES];
    target[0] = new double [m_ncategories*BASES];
    long cat;

    for (cat = 0; cat < m_ncategories; ++cat)
    {
        basefreqarray[cat]   = basefreqarray[0] + cat*BASES;
        daughter1[cat] = daughter1[0] + cat*BASES;
        daughter2[cat] = daughter2[0] + cat*BASES;
        target[cat] = target[0] + cat*BASES;
        basefreqarray[cat][baseA] = m_basefreqs[baseA];
        basefreqarray[cat][baseC] = m_basefreqs[baseC];
        basefreqarray[cat][baseG] = m_basefreqs[baseG];
        basefreqarray[cat][baseT] = m_basefreqs[baseT];
    }

} // F84Model::AllocateBuffers

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

F84Model::~F84Model()
{
    // if the array does not exist we suppose that the object is
    // not Finalized, and don't try any deallocations.
    EmptyBuffers();

} // F84Model::~F84Model

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

void F84Model::CopyMembers(const F84Model& src)
{
    m_ttratio   = src.m_ttratio;

    freqar    = src.freqar;
    freqcy    = src.freqcy;
    freqgr    = src.freqgr;
    freqty    = src.freqty;

    computed1 = src.computed1;
    computed2 = src.computed2;

    // don't copy the data likelihood buffers
}

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

bool F84Model::IsValidDataModel() const
{
    if (m_ttratio <= 0.5) return false;
    return NucModel::IsValidDataModel();

} // F84Model::IsValidDataModel

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

void F84Model::CopyBuffers(const F84Model& src)
{
    xcatrates = src.xcatrates;
    ycatrates = src.ycatrates;
    expA1     = src.expA1;
    expB1     = src.expB1;
    expC1     = src.expC1;
    expA2     = src.expA2;
    expB2     = src.expB2;
    expC2     = src.expC2;
    catlikes = src.catlikes;

    EmptyBuffers();
    AllocateBuffers();

    // Do not try to copy the buffers if they are unallocated
    if (src.daughter1)
        memcpy(daughter1[0],src.daughter1[0],m_ncategories*BASES*sizeof(double));
    if (src.daughter2)
        memcpy(daughter2[0],src.daughter2[0],m_ncategories*BASES*sizeof(double));
    if (src.target)
        memcpy(target[0],src.target[0],m_ncategories*BASES*sizeof(double));

} // F84Model::CopyBuffers

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

F84Model::F84Model(const F84Model& src)
    : NucModel(src), basefreqarray(NULL), daughter1(NULL), daughter2(NULL),
      target(NULL)
{
    //cout << "*****F84 model copied*****" << this << endl;

    CopyMembers(src);
    CopyBuffers(src);
}

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

F84Model& F84Model::operator=(const F84Model& src)
{
    NucModel::operator=(src);
    CopyMembers(src);
    CopyBuffers(src);
    return *this;
}

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

void F84Model::SetTTratio(double tr)
{
    if (tr <= 0.5)
    {
        data_error e("Transition/transversion ratio must be > 0.5");
        throw e;
    }
    m_ttratio = tr;
}

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

void F84Model::Finalize()
{
    double pur, pyr, ag, ct, m, n, x, y;
    long cat;

    double freqa(m_basefreqs[baseA]), freqc(m_basefreqs[baseC]),
        freqg(m_basefreqs[baseG]), freqt(m_basefreqs[baseT]);

    pur    = freqa + freqg;
    pyr    = freqc + freqt;
    freqar = freqa / pur;
    freqcy = freqc / pyr;
    freqgr = freqg / pur;
    freqty = freqt / pyr;

    ag = freqa * freqg;
    ct = freqc * freqt;
    m  = m_ttratio * pur * pyr - (ag + ct);
    n  = ag / pur + ct / pyr;

    // code stolen from DNAMLK here
    y = m / (m + n);
    // if globalfreqs has failed....
    if (y < 0) y = 0;
    else if (y > 1) y = 1;
    x = 1.0 - y;
    double fracchange = y * (2.0 * freqa * freqgr + 2.0 * freqc * freqty)
        + x * (1.0 - freqa * freqa - freqc * freqc - freqg * freqg - freqt * freqt);
    y /= - fracchange;
    x /= - fracchange;

    for (cat = 0; cat < m_ncategories; ++cat)
    {
        xcatrates.push_back(x*m_catrates[cat]);
        ycatrates.push_back(y*m_catrates[cat]);
    }

    // Allocate additional space for likelihood calculations.
    expA1.insert(expA1.begin(),m_ncategories,0.0);
    expA2.insert(expA2.begin(),m_ncategories,0.0);
    expB1.insert(expB1.begin(),m_ncategories,0.0);
    expB2.insert(expB2.begin(),m_ncategories,0.0);
    expC1.insert(expC1.begin(),m_ncategories,0.0);
    expC2.insert(expC2.begin(),m_ncategories,0.0);

    AllocateBuffers();

    double zero = 0.0;
    catlikes = CreateVec2d(m_nmarkers,m_ncategories,zero);

} // F84Model::Finalize

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

void F84Model::RescaleLength1(double length1)
{
    long cat;
    double n;

    if (length1 > FLAGDOUBLE)
    {
#ifdef DEBUG_VARIABLES
        if (printdbgvar > 0)
        {
            cout << "RescaleLength1 raw length1: " << length1 << endl;
            //cout << " m_relmurate: " << m_relmurate << endl;
        }
#endif
        length1 *= m_relmurate;
        computed1 = true;

        for(cat = 0; cat < m_ncategories; cat++)
        {
            n = exp(length1 * xcatrates[cat]);
            expA1[cat] = 1.0 - n;
            expB1[cat] = n * exp(length1 * ycatrates[cat]);
            expC1[cat] = n - expB1[cat];
#ifdef DEBUG_VARIABLES
            if (0)
            //if (printdbgvar > 0)
            {
                cout << " cat: " << cat << endl;
                cout << " n: " << n << endl;
                cout << " xcatrates: " << xcatrates[cat] << endl;
                cout << " ycatrates: " << ycatrates[cat] << endl;
                cout << " expA1: "<< expA1[cat] << endl;
                cout << " expB1: "<< expB1[cat] << endl;
                cout << " expC1: "<< expC1[cat] << endl;
            }
#endif
        }
    }
    else
    {
        computed1 = false;
    }
#ifdef DEBUG_VARIABLES
    if (printdbgvar > 0)
    {
        cout << " scaled length1: " << length1 << " computed1: "<< computed1 << endl << endl;
    }
#endif

} // F84Model::RescaleLength1

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

void F84Model::RescaleLength2(double length2)
{
    long cat;
    double n;

    if (length2 > FLAGDOUBLE)
    {
#ifdef DEBUG_VARIABLES
        if (printdbgvar > 0)
        {
            cout << "RescaleLength2 raw length2: " << length2 << endl;
            //cout << " m_relmurate: " << m_relmurate << endl;
        }
#endif
        length2 *= m_relmurate;
        computed2 = true;

        for(cat = 0; cat < m_ncategories; cat++)
        {
            n = exp(length2 * xcatrates[cat]);
            expA2[cat] = 1.0 - n;
            expB2[cat] = n * exp(length2 * ycatrates[cat]);
            expC2[cat] = n - expB2[cat];
#ifdef DEBUG_VARIABLES
            if (0)
            //if (printdbgvar > 0)
            {
                cout << " cat: " << cat << endl;
                cout << " n: " << n << endl;
                cout << " xcatrates: " << xcatrates[cat] << endl;
                cout << " ycatrates: " << ycatrates[cat] << endl;
                cout << " expA2: "<< expA2[cat] << endl;
                cout << " expB2: "<< expB2[cat] << endl;
                cout << " expC2: "<< expC2[cat] << endl;
            }
#endif
        }
    }
    else
    {
        computed2 = false;
    }
#ifdef DEBUG_VARIABLES
    if (printdbgvar > 0)
    {
        cout << " scaled length2: " << length2 << " computed2: "<< computed2 << endl << endl;
    }
#endif
}

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

void F84Model::RescaleLengths(double length1, double length2)
{
    RescaleLength1(length1);
    RescaleLength2(length2);

} // F84Model::RescaleLengths

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

double** F84Model::ComputeSiteDLs(double** siteDLs1, double** siteDLs2)
{
    double sumAll, sumPur, sumPyr;

    long cat;
#ifdef DEBUG_VARIABLES
    if (0)
    //if (printdbgvar > 0)
    {
        cout << "basefreqA: "<< m_basefreqs[baseA] << endl;
        cout << "basefreqC: "<< m_basefreqs[baseC] << endl;
        cout << "basefreqG: "<< m_basefreqs[baseG] << endl;
        cout << "basefreqT: "<< m_basefreqs[baseT] << endl;
        if (computed1)
            cout << "computed1: true"<< endl;
        else
            cout << "computed1: false"<< endl;
        if (computed2)
            cout << "computed2: true"<< endl;
        else
            cout << "computed2: false"<< endl;
    }
#endif

    if (computed1)
    {
        for (cat = 0; cat < m_ncategories; ++cat)
        {
            double* catDLs1 = siteDLs1[cat];
            sumAll  = expA1[cat] *
                (m_basefreqs[baseA]*catDLs1[baseA] +
                 m_basefreqs[baseC]*catDLs1[baseC] +
                 m_basefreqs[baseG]*catDLs1[baseG] +
                 m_basefreqs[baseT]*catDLs1[baseT]);

            sumPur  = freqar*catDLs1[baseA] + freqgr*catDLs1[baseG];
            sumPyr  = freqcy*catDLs1[baseC] + freqty*catDLs1[baseT];

#ifdef DEBUG_VARIABLES
            if (printdbgvar > 0)
            {
                //cout<< "cat: " << cat << endl;
                //cout<< "sumA1: " << sumAll << endl;
                //cout<< "freqar: " << freqar << endl;
                //cout<< "freqcy: " << freqcy << endl;
                //cout<< "freqgr: " << freqgr << endl;
                //cout<< "freqty: " << freqty << endl;
                //cout<< "sumPur1: " << sumPur << endl;
                //cout<< "sumPyr1: " << sumPyr << endl;
                cout<< "dls1: " << catDLs1[baseA] << " "<< catDLs1[baseC] << " "<< catDLs1[baseG] << " "<< catDLs1[baseT]  << endl;
                //cout<< "dls1A: " << catDLs1[baseA] << endl;
                //cout<< "dls1C: " << catDLs1[baseC] << endl;
                //cout<< "dls1G: " << catDLs1[baseG] << endl;
                //cout<< "dls1T: " << catDLs1[baseT] << endl;
                //cout<< "expA1: " << expA1[cat] <<endl;
                //cout<< "expB1: " << expB1[cat] <<endl;
                //cout<< "expC1: " << expC1[cat] <<endl<<endl;
            }
#endif
            double expC1cat = expC1[cat];
            daughter1[cat][baseA] = sumAll + expB1[cat]*catDLs1[baseA] +
                expC1cat*sumPur;
            daughter1[cat][baseC] = sumAll + expB1[cat]*catDLs1[baseC] +
                expC1cat*sumPyr;
            daughter1[cat][baseG] = sumAll + expB1[cat]*catDLs1[baseG] +
                expC1cat*sumPur;
            daughter1[cat][baseT] = sumAll + expB1[cat]*catDLs1[baseT] +
                expC1cat*sumPyr;
        }
    }
    else
    {
        memcpy(daughter1[0],basefreqarray[0],m_ncategories*BASES*sizeof(double));
    }

    if (computed2)
    {
        for (cat = 0; cat < m_ncategories; cat++)
        {
            double* catDLs2 = siteDLs2[cat];
            sumAll  = expA2[cat] *
                (m_basefreqs[baseA]*catDLs2[baseA] +
                 m_basefreqs[baseC]*catDLs2[baseC] +
                 m_basefreqs[baseG]*catDLs2[baseG] +
                 m_basefreqs[baseT]*catDLs2[baseT]);

            sumPur  = freqar*catDLs2[baseA] + freqgr*catDLs2[baseG];
            sumPyr  = freqcy*catDLs2[baseC] + freqty*catDLs2[baseT];

#ifdef DEBUG_VARIABLES
            if (printdbgvar > 0)
            {
                //cout<< "cat: " << cat << endl;
                //cout<< "sumA2: " << sumAll << endl;// JRM debug
                //cout<< "sumPur2: " << sumPur << endl;
                //cout<< "sumPyr2: " << sumPyr << endl;
                //cout<< "dls2A: " << catDLs2[baseA] << endl;
                //cout<< "dls2C: " << catDLs2[baseC] << endl;
                //cout<< "dls2G: " << catDLs2[baseG] << endl;
                //cout<< "dls2T: " << catDLs2[baseT] << endl;
                //cout<< "expA2: " << expA2[cat] <<endl;
                //cout<< "expB2: " << expB2[cat] <<endl;
                //cout<< "expC2: " << expC2[cat] <<endl<<endl;
                cout<< "dls2: " << catDLs2[baseA] << " "<< catDLs2[baseC] << " "<< catDLs2[baseG] << " "<< catDLs2[baseT]  << endl;
           }
#endif

            double expC2cat = expC2[cat];
            daughter2[cat][baseA] = sumAll + expB2[cat]*catDLs2[baseA] +
                expC2cat*sumPur;
            daughter2[cat][baseC] = sumAll + expB2[cat]*catDLs2[baseC] +
                expC2cat*sumPyr;
            daughter2[cat][baseG] = sumAll + expB2[cat]*catDLs2[baseG] +
                expC2cat*sumPur;
            daughter2[cat][baseT] = sumAll + expB2[cat]*catDLs2[baseT] +
                expC2cat*sumPyr;
        }
    }
    else
    {
        memcpy(daughter2[0],basefreqarray[0],m_ncategories*BASES*sizeof(double));
    }

    for (cat = 0; cat < m_ncategories; cat++)
    {
        target[cat][baseA] = daughter1[cat][baseA] *
            daughter2[cat][baseA];
        target[cat][baseC] = daughter1[cat][baseC] *
            daughter2[cat][baseC];
        target[cat][baseG] = daughter1[cat][baseG] *
            daughter2[cat][baseG];
        target[cat][baseT] = daughter1[cat][baseT] *
            daughter2[cat][baseT];
 #ifdef DEBUG_VARIABLES
        if (printdbgvar > 0)
        {
            //cout<< "cat: " << cat << endl;
            //cout<< "targetA: " << target[cat][baseA] <<endl;
            //cout<< "targetC: " << target[cat][baseC] <<endl;
            //cout<< "targetG: " << target[cat][baseG] <<endl;
            //cout<< "targetT: " << target[cat][baseT] <<endl<<endl;
                cout<< "target: " << target[cat][baseA] << " "<< target[cat][baseC] << " "<< target[cat][baseG] << " "<< target[cat][baseT]  << endl;
        }
#endif
  }

    return target;

} // F84Model::ComputeSiteDLs

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

double F84Model::ComputeSubtreeDL(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    double total=0.0, subtotal;
    double** marker;
    long cat;
    DoubleVec1d prior(m_ncategories);
    long firstposn = posn;
    int im = 0;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {
        subtotal = 0.0;

        for (cat = 0; cat < m_ncategories; cat++)
        {
            prior[cat] = m_basefreqs[baseA]*marker[cat][baseA] +
                m_basefreqs[baseC]*marker[cat][baseC] +
                m_basefreqs[baseG]*marker[cat][baseG] +
                m_basefreqs[baseT]*marker[cat][baseT];

            subtotal += m_catprobs[cat] * prior[cat];
#if 0
            cout << "marker " << im
                 << " cat " << cat
                 << " probs " << marker[cat][baseA]
                 << " " << marker[cat][baseC]
                 << " " << marker[cat][baseG]
                 << " " << marker[cat][baseT]
                 << " prior " << prior[cat]
                 << endl;
#endif
        }
        im++;

        if (!subtotal)
        {
            DataModel::TryToNormalizeAndThrow(posn, GetModelType());
        }

        if (ShouldNormalize())
        {
            total += (log(subtotal) + rootdls.GetNorms(posn));
        }
        else
        {
            total +=  log(subtotal);
        }

        // de-normalization not needed here, since we are only interested in
        // the ratio
        if (m_ncategories > 1)
        {
            for (cat = 0; cat < m_ncategories; cat++)
                catlikes[posn][cat] = prior[cat]/subtotal;
        }
        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);
    //    cout << "total: " << total << endl;
    return total;

} // F84Model::ComputeSubtreeDL

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

double F84Model::ComputeCatDL(long startmarker, long endmarker)
{
    double subtotal;
    long marker, cat;

    DoubleVec1d like = m_catcells;
    DoubleVec1d nulike(m_ncategories);

    for (marker = startmarker; marker != endmarker; ++marker)
    {
        subtotal = 0.0;
        for (cat = 0; cat < m_ncategories; cat++)
            subtotal += m_catprobs[cat] * like[cat];

        subtotal *= m_acratio;

        for (cat = 0; cat < m_ncategories; cat++)
            nulike[cat] = catlikes[marker][cat] *
                (subtotal + m_ucratio * like[cat]);

        // the following puts the nulike values into like.  It
        // also puts the like values into nulike, but we will not
        // be using values from nulike so we don't care.
        like.swap(nulike);
    }

    subtotal = 0.0;
    for (cat = 0; cat < m_ncategories; cat++)
        subtotal += m_catprobs[cat] * like[cat];

    // the following puts the like values into catcells for
    // long-term storage.  It also puts the catcells values into
    // like, but we don't care.
    m_catcells.swap(like);

    return log(subtotal);

} // F84Model::ComputeCatDL

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

StringVec1d F84Model::CreateDataModelReport() const
{
    StringVec1d report = NucModel::CreateDataModelReport();

    string line = "Transition/transversion ratio: " + ToString(m_ttratio);
    report.push_back(line);

    return report;

} // F84Model::CreateDataModelReport

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

StringVec1d F84Model::ToXML(size_t nspaces) const
{
    StringVec1d xmllines(NucModel::ToXML(nspaces));

    nspaces += INDENT_DEPTH;
    string line(MakeIndent(MakeTag(xmlstr::XML_TAG_TTRATIO),nspaces));
    line += ToString(GetTTratio());
    line += MakeCloseTag(xmlstr::XML_TAG_TTRATIO);
    nspaces -= INDENT_DEPTH;

    StringVec1d::iterator endtag = --xmllines.end();
    xmllines.insert(endtag,line);

    return xmllines;

} // F84Model::ToXML

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

// This routine simulates data for a single node under the F84
// model.  It assumes that the given branch lengths have already
// been rescaled for the rate category desired.
// Nomenclature follows _Inferring Phylogenies_ pp. 202-203,
//  printing date 2004 (earlier printings have different page numberings)

// The second argument is unused since nucleotides don't have differing
// number of states by marker position
DoubleVec1d F84Model::SimulateMarker(double branchlength, long whichmarker, const DoubleVec1d& state) const
{
    double freqa(m_basefreqs[baseA]), freqc(m_basefreqs[baseC]), freqg(m_basefreqs[baseG]), freqt(m_basefreqs[baseT]);
    double pur = freqa + freqg;
    double pyr = freqc + freqt;
    double beta = 1.0 / (2.0 * pur * pyr * (1.0 + m_ttratio));
    double alpha = (pur * pyr * m_ttratio - freqa * freqg - freqc * freqt) /
        (2.0 * (1.0 + m_ttratio) * (pyr * freqa * freqg + pur * freqc * freqt));

    DoubleVec1d cumprob;
    cumprob.push_back(freqa);
    cumprob.push_back(freqa + freqc);
    cumprob.push_back(freqa + freqc + freqg);

    double expB = exp(-beta * branchlength);

    DoubleVec1d answer(BASES, 0.0);

    double general = 1.0 - expB;
    double r = registry.GetRandom().Float();

    // We compute a chance of drawing from a pool of all four nucleotides
    // in proportion to their frequency; if that doesn't happen, we
    // compute a chance of drawing from a pool of only purines or only
    // pyrimidines; if that doesn't happen either there is no change
    // and we return the initial state.

    if (r < general)
    {
        // perform a draw from the general pool
        double r2 = registry.GetRandom().Float();
        long i;
        for (i = 0; i < BASES-1; ++i)
        {
            if (r2 < cumprob[i])
            {
                answer[i] = 1.0;
                return answer;
            }
        }
        // flowthrough if previous bases not picked
        answer[BASES-1] = 1.0;
        return answer;
    }
    else
    {
        r -= general;
        double transition = expB * (1.0 - exp(-alpha * branchlength));
        if (r < transition)
        {
            // perform a draw from the transition pool
            double r2 = registry.GetRandom().Float();
            if (state[baseA] == 1.0 || state[baseG] == 1.0) // purine
            {
                if (r2 < freqa / (freqa + freqg))
                {
                    answer[baseA] = 1.0;
                }
                else
                {
                    answer[baseG] = 1.0;
                }
            }
            else                        // pyrimidine
            {
                if (r2 < freqc / (freqc + freqt))
                {
                    answer[baseC] = 1.0;
                }
                else
                {
                    answer[baseT] = 1.0;
                }
            }
            return answer;
        }
    }

    // otherwise, no event happens
    return state;

} // SimulateMarke

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

GTRModel::GTRModel(long nmarkers,
                   long numCategories,
                   DoubleVec1d categoryRates,
                   DoubleVec1d categoryProbabilities,
                   double userAutoCorrelationValue,
                   bool doNormalize,
                   double relMuRate,
                   double freqA,
                   double freqC,
                   double freqG,
                   double freqT,
                   double freqAC,
                   double freqAG,
                   double freqAT,
                   double freqCG,
                   double freqCT,
                   double freqTG,
                   double perBaseErrorRate)
    : NucModel(nmarkers,
               numCategories,
               categoryRates,
               categoryProbabilities,
               userAutoCorrelationValue,
               doNormalize,
               relMuRate,
               freqA,
               freqC,
               freqG,
               freqT,
               false,          // calculate freqs from data
               perBaseErrorRate),
      AC(freqAC),
      AG(freqAG),
      AT(freqAT),
      CG(freqCG),
      CT(freqCT),
      TG(freqTG),
      // set some potentially useful defaults
      basefreqarray(NULL),
      daughter(NULL),target(NULL)
{
    // reset the names from the base values
    computed.assign(2,true);
    DoubleVec1d empty(BASES,0.0);
    scratch.assign(BASES,empty);
    Finalize();

    assert(GTRModel::IsValidDataModel());
} // GTRModel::GTRModel

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

GTRModel::GTRModel(const GTRModel& src)
    : NucModel(src), basefreqarray(NULL),daughter(NULL),target(NULL)
{
    CopyMembers(src); // doesn't copy target or daughter buffers!
} // GTRModel copy ctor

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

GTRModel::~GTRModel()
{
    EmptyBuffers();
} // GTRModel dtor

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

DataModel* GTRModel::Clone() const
{
    DataModel* newmodel = new GTRModel(*this);
    return newmodel;
}

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

GTRModel& GTRModel::operator=(const GTRModel& src)
{
    NucModel::operator=(src);
    CopyMembers(src);
    return *this;
}

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

void GTRModel::AllocateBuffers()
{
    DoubleVec1d zero(m_ncategories,0.0);
    catlikes.assign(m_nmarkers,zero);

    DoubleVec1d base1(BASES,0);
    DoubleVec2d base2(BASES,base1);
    DoubleVec3d empty(m_ncategories,base2);
    pchange.assign(2,empty);

    daughter         = new double** [2];
    target           = new double* [m_ncategories];
    basefreqarray    = new double* [m_ncategories];
    daughter[0]      = new double* [2*m_ncategories];
    daughter[0][0]   = new double [2*m_ncategories*BASES];
    target[0]        = new double [m_ncategories*BASES];
    basefreqarray[0] = new double [m_ncategories*BASES];
    daughter[1] = daughter[0] + m_ncategories;
    daughter[1][0] = daughter[0][0] + m_ncategories*BASES;
    long cat;
    for (cat = 0; cat < m_ncategories; ++cat)
    {
        basefreqarray[cat]   = basefreqarray[0] + cat*BASES;
        daughter[0][cat] = daughter[0][0] + cat*BASES;
        daughter[1][cat] = daughter[1][0] + cat*BASES;
        target[cat] = target[0] + cat*BASES;
        long base;
        for(base = baseA; base <= baseT; ++base)
            basefreqarray[cat][base] = m_basefreqs[base];
    }
} // GTRModel::AllocateBuffers

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

void GTRModel::EmptyBuffers()
{
    if (basefreqarray)
    {
        delete [] basefreqarray[0];
        delete [] basefreqarray;
    }
    if (daughter)
    {
        delete [] daughter[0][0];
        delete [] daughter[0];
        delete [] daughter;
    }
    if (target)
    {
        delete [] target[0];
        delete [] target;
    }

} // GTRModel::EmptyBuffers

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

void GTRModel::CopyMembers(const GTRModel& src)
{
    AG = src.AG;
    AC = src.AC;
    AT = src.AT;
    CG = src.CG;
    CT = src.CT;
    TG = src.TG;
    eigvals = src.eigvals;
    eigvecs1 = src.eigvecs1;
    eigvecs2 = src.eigvecs2;
    pchange = src.pchange;
    computed = src.computed;
    scratch = src.scratch;
    m_basefreqs = src.m_basefreqs;

    if (src.basefreqarray)
    {
        EmptyBuffers();
        AllocateBuffers();

        memcpy(basefreqarray[0],src.basefreqarray[0],
               m_ncategories*BASES*sizeof(double));
    }

    // we don't copy target or daughter!
} // GTRModel::CopyMembers

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

void GTRModel::GTRDotProduct(const DoubleVec2d& first, const DoubleVec2d& second, DoubleVec2d& answer)
// dot product of first and second put into PRE-EXISTING answer!
// second has been PRE-TRANSPOSED!!
{
    // should be square
    assert(first.size() == first[0].size());
    // and all the same size!
    assert(first.size() == second.size());
    assert(first.size() == answer.size());
    double initial = 0.0;

    long i, j, n=first.size();
    for (i = 0; i < n; ++i)
    {
        for (j = 0; j < n; ++j)
        {
            answer[i][j] = inner_product(first[i].begin(), first[i].end(), second[j].begin(), initial);
        }
    }
} // GTRModel::GTRDotProduct

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

void GTRModel::BarfOnBadGTRRates(const DoubleVec1d& rts) const
{
    if (rts.size() != 6)
    {
        data_error e("Incorrect number of GTR rates:  expected 6, found " + ToString(rts.size()));
        throw e;
    }

    size_t i;
    for (i = 0; i < rts.size(); ++i)
    {
        if (rts[i] <= 0.0)
        {
            data_error e("All rates for the GTR model must be greater than 0.");
            throw e;
        }
    }
} // GTRModel::BarfOnBadGTRRates

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

void GTRModel::SetRates(const DoubleVec1d& rts)
{
    BarfOnBadGTRRates(rts);

    AC = rts[0];
    AG = rts[1];
    AT = rts[2];
    CG = rts[3];
    CT = rts[4];
    TG = rts[5];
} // GTRModel::SetRates

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

DoubleVec1d GTRModel::GetRates() const
{
    DoubleVec1d rts;
    rts.push_back(AC);
    rts.push_back(AG);
    rts.push_back(AT);
    rts.push_back(CG);
    rts.push_back(CT);
    rts.push_back(TG);

    return rts;

} // GTRModel::GetRates

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

void GTRModel::Finalize()
{
    AllocateBuffers();

    // calculate eigvals and eigvecs1 & 2
    // assemble rate matrix; using scratch
    double* bf = basefreqarray[0];

    // rescale rates to a mean of 1 event per unit branchlength
    double scalefactor = 2.0 * (AC * bf[baseA] * bf[baseC] +
                                AG * bf[baseA] * bf[baseG] +
                                AT * bf[baseA] * bf[baseT] +
                                CG * bf[baseC] * bf[baseG] +
                                CT * bf[baseC] * bf[baseT] +
                                TG * bf[baseG] * bf[baseT]);
    double nAC = AC / scalefactor;
    double nAG = AG / scalefactor;
    double nAT = AT / scalefactor;
    double nCG = CG / scalefactor;
    double nCT = CT / scalefactor;
    double nTG = TG / scalefactor;

    scratch[0][0] = -(nAC*bf[baseC] + nAG*bf[baseG] + nAT*bf[baseT])/bf[baseA];
    scratch[0][1] = nAC;
    scratch[0][2] = nAG;
    scratch[0][3] = nAT;
    scratch[1][0] = nAC;
    scratch[1][1] = -(nAC*bf[baseA] + nCG*bf[baseG] + nCT*bf[baseT])/bf[baseC];
    scratch[1][2] = nCG;
    scratch[1][3] = nCT;
    scratch[2][0] = nAG;
    scratch[2][1] = nCG;
    scratch[2][2] = -(nAG*bf[baseA] + nCG*bf[baseC] + nTG*bf[baseT])/bf[baseG];
    scratch[2][3] = nTG;
    scratch[3][0] = nAT;
    scratch[3][1] = nCT;
    scratch[3][2] = nTG;
    scratch[3][3] = -(nAT*bf[baseA] + nCT*bf[baseC] + nTG*bf[baseG])/bf[baseT];

    // assemble square root of base frequencies
    DoubleVec1d zero(BASES,0.0);
    DoubleVec2d diag(BASES,zero);
    long i;
    for (i = 0; i < BASES; ++i)
    {
        diag[i][i] = sqrt(bf[i]);
    }

    GTRDotProduct(diag,Transpose(scratch),scratch);
    DoubleVec2d answ(BASES,zero);
    GTRDotProduct(scratch,diag,answ);

    EigenCalculator eig;
    pair<DoubleVec1d,DoubleVec2d> eigsys(eig.Eigen(answ));
    eigvecs1 = diag;  // needed to satisfy GTRDotProduct size reqs
    // calculate: diag . Transpose[eigsys.second], which would require
    // a call to the transpose of the transpose for GTRDotProduct purposes,
    // so just use the matrix straight.
    GTRDotProduct(diag,eigsys.second,eigvecs1);
    eigvals = eigsys.first;
    eigvecs2 = Transpose(Invert(eigvecs1));

} // GTRModel::Finalize

//------------------------------------------------------------------------------------
// calculate pchange [eigvecs1 . Exp[eigvals*length] . eigvecs2]

void GTRModel::RescaleLengths(double length1, double length2)
{
    computed.assign(2,true);

    // if a length is FLAGDOUBLE we are in a one-legged coalescence
    // and we avoid this whole computation
    if (CloseEnough(length1, FLAGDOUBLE)) computed[0] = false;
    if (CloseEnough(length2, FLAGDOUBLE)) computed[1] = false;

    DoubleVec1d zero(BASES,0.0);
    DoubleVec2d zeros(BASES,zero);
    DoubleVec3d zeroes(m_ncategories,zeros);
    DoubleVec4d diag(2,zeroes);

    long cat;
    double expon;

    length1 *= m_relmurate;
    length2 *= m_relmurate;

    for(cat = 0; cat < m_ncategories; ++cat)
    {
        long i;
        for (i = 0; i < BASES; ++i)
        {
            double scalar = m_catrates[cat]*eigvals[i];

            if (computed[0])
            {
                expon = length1 * scalar;
                if (expon >= EXPMIN)
                {
                    if (expon <= EXPMAX)
                        diag[0][cat][i][i] = exp(length1*scalar);
                    else
                        diag[0][cat][i][i] = EXP_OF_EXPMAX;
                }
                // else it remains zero
            }

            if (computed[1])
            {
                expon = length2 * scalar;
                if (expon >= EXPMIN)
                {
                    if (expon <= EXPMAX)
                        diag[1][cat][i][i] = exp(length2*scalar);
                    else
                        diag[1][cat][i][i] = EXP_OF_EXPMAX;
                }
                // else it remains zero
            }
        }
    }

    for(cat = 0; cat < m_ncategories; ++cat)
    {
        long br;
        for(br = 0; br < 2; ++br)
        {
            if (!computed[br]) continue;
            GTRDotProduct(eigvecs1,diag[br][cat],scratch);
            GTRDotProduct(scratch,eigvecs2,pchange[br][cat]);
        }
    }

} // GTRModel::RescaleLengths

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

double** GTRModel::ComputeSiteDLs(double** siteDL1, double** siteDL2)
{
    long cat, base;

    if (computed[0])
    {
        for(cat = 0; cat < m_ncategories; ++cat)
        {
            double* catDLs = siteDL1[cat];
            DoubleVec2d& prob = pchange[0][cat];
            for(base = baseA; base <= baseT; ++base)
            {
                daughter[0][cat][base] = prob[baseA][base] * catDLs[baseA] +
                    prob[baseC][base] * catDLs[baseC] +
                    prob[baseG][base] * catDLs[baseG] +
                    prob[baseT][base] * catDLs[baseT];
            }
        }
    }
    else
    {
        memcpy(daughter[0][0],basefreqarray[0],m_ncategories*BASES*sizeof(double));
    }

    if (computed[1])
    {
        for(cat = 0; cat < m_ncategories; ++cat)
        {
            double* catDLs = siteDL2[cat];
            DoubleVec2d& prob = pchange[1][cat];
            for(base = baseA; base <= baseT; ++base)
            {
                daughter[1][cat][base] = prob[baseA][base] * catDLs[baseA] +
                    prob[baseC][base] * catDLs[baseC] +
                    prob[baseG][base] * catDLs[baseG] +
                    prob[baseT][base] * catDLs[baseT];
            }
        }
    }
    else
    {
        memcpy(daughter[1][0],basefreqarray[0],m_ncategories*BASES*sizeof(double));
    }

    for (cat = 0; cat < m_ncategories; cat++)
        for(base = baseA; base <= baseT; ++base)
            target[cat][base] = daughter[0][cat][base] *
                daughter[1][cat][base];

    return target;

} // GTRModel::ComputeSiteDLs

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

double GTRModel::ComputeSubtreeDL(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    double total=0.0, subtotal;
    double** marker;
    long cat;
    DoubleVec1d prior(m_ncategories);
    long firstposn = posn;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {
        subtotal = 0.0;

        for (cat = 0; cat < m_ncategories; cat++)
        {
            prior[cat] = basefreqarray[cat][baseA]*marker[cat][baseA] +
                basefreqarray[cat][baseC]*marker[cat][baseC] +
                basefreqarray[cat][baseG]*marker[cat][baseG] +
                basefreqarray[cat][baseT]*marker[cat][baseT];

            subtotal += m_catprobs[cat] * prior[cat];
        }

        if (!subtotal)
        {
            DataModel::TryToNormalizeAndThrow(posn, GetModelType());
        }

        if (ShouldNormalize())
        {
            total += (log(subtotal) + rootdls.GetNorms(posn));
        }
        else
        {
            total += log(subtotal);
        }

        // de-normalization not needed here, since we are only interested in the ratio
        if (m_ncategories > 1)
        {
            for (cat = 0; cat < m_ncategories; cat++)
                catlikes[posn][cat] = prior[cat]/subtotal;
        }
        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);
    return total;

} // GTRModel::ComputeSubtreeDL

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

double GTRModel::ComputeCatDL(long startmarker, long endmarker)
{
    double subtotal;
    long marker, cat;

    DoubleVec1d like = m_catcells;
    DoubleVec1d nulike(m_ncategories);

    for (marker = startmarker; marker != endmarker; ++marker)
    {
        subtotal = 0.0;
        for (cat = 0; cat < m_ncategories; cat++)
            subtotal += m_catprobs[cat] * like[cat];

        subtotal *= m_acratio;

        for (cat = 0; cat < m_ncategories; cat++)
            nulike[cat] = catlikes[marker][cat] *
                (subtotal + m_ucratio * like[cat]);

        // the following puts the nulike values into like.  It
        // also puts the like values into nulike, but we will not
        // be using values from nulike so we don't care.
        like.swap(nulike);
    }

    subtotal = 0.0;
    for (cat = 0; cat < m_ncategories; cat++)
        subtotal += m_catprobs[cat] * like[cat];

    // the following puts the like values into catcells for
    // long-term storage.  It also puts the catcells values into
    // like, but we don't care.
    m_catcells.swap(like);

    return log(subtotal);

} // GTRModel::ComputeCatDL

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

StringVec1d GTRModel::CreateDataModelReport() const
{
    StringVec1d report = NucModel::CreateDataModelReport();

    string line = "Mutation parameters: ";
    report.push_back(line);
    line = "Between A and (C, G, T):  ";
    line += ToString(AC) + ", " + ToString(AG) + ", " + ToString(AT);
    report.push_back(line);
    line = "Between C and (G, T):  ";
    line += ToString(CG) + ", " + ToString(CT);
    report.push_back(line);
    line = "Between G and (T):  ";
    line += ToString(TG);
    report.push_back(line);

    return report;

} // GTRModel::CreateDataModelReport

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

StringVec1d GTRModel::ToXML(size_t nspaces) const
{
    StringVec1d xmllines(NucModel::ToXML(nspaces));

    nspaces += INDENT_DEPTH;
    string line(MakeIndent(MakeTag(xmlstr::XML_TAG_GTRRATES),nspaces));
    line += ToString(GetRates(),6);
    line += MakeCloseTag(xmlstr::XML_TAG_GTRRATES);
    nspaces -= INDENT_DEPTH;

    StringVec1d::iterator endtag = --xmllines.end();
    xmllines.insert(endtag,line);

    return xmllines;

} // GTRModel::ToXML

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

DoubleVec1d GTRModel::SimulateMarker(double branchlength, long whichmarker,
                                     const DoubleVec1d& state) const
{
    throw implementation_error("Cannot simulate data with the GTR model yet.");
} // SimulateMarker

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

bool GTRModel::IsValidDataModel() const
{
    if (AC <= 0) return false;
    if (AG <= 0) return false;
    if (AT <= 0) return false;
    if (CG <= 0) return false;
    if (CT <= 0) return false;
    if (TG <= 0) return false;
    return NucModel::IsValidDataModel();
} // GTRModel::IsValidDataModel

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

AlleleModel::AlleleModel(long nmarkers,
                         long numCategories,
                         DoubleVec1d categoryRates,
                         DoubleVec1d categoryProbabilities,
                         double userAutoCorrelationValue,
                         bool doNormalize,
                         long numBins,
                         double relMuRate)
    : DataModel(nmarkers,
                numCategories,
                categoryRates,
                categoryProbabilities,
                userAutoCorrelationValue,
                doNormalize,
                numBins,
                relMuRate),
      m_likes(CreateVec2d(m_nmarkers, m_ncategories, 0.0)),
      m_bincounts()
{
    // intentionally blank
} // AlleleModel constructor

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

DoubleVec1d AlleleModel::RescaleLength(double length)
{
    length *= m_relmurate;

    long cat;
    DoubleVec1d scaled(m_ncategories);
    for (cat = 0; cat < m_ncategories; ++cat)
    {
        scaled[cat] = length * m_catrates[cat];
    }
    return scaled;

} // AlleleModel::RescaleLengths

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

double AlleleModel::ComputeCatDL(long startmarker, long endmarker)
{
    double subtotal;
    long marker, cat;

    DoubleVec1d previous = m_catcells;
    DoubleVec1d current(m_ncategories, 0.0);

    for (marker = startmarker; marker != endmarker; ++marker)
    {
        subtotal = 0.0;
        for (cat = 0; cat < m_ncategories; cat++)
            subtotal += m_catprobs[cat] * previous[cat];

        subtotal *= m_acratio;

        for (cat = 0; cat < m_ncategories; cat++)
            current[cat] = m_likes[marker][cat] *
                (subtotal + m_ucratio * previous[cat]);

        // This line puts the current values in previous cheaply.  We
        // don't care that it puts the previous values in current,
        // because current will be overwritten anyway.
        previous.swap(current);
    }

    subtotal = 0.0;
    for (cat = 0; cat < m_ncategories; cat++)
        subtotal += m_catprobs[cat] * previous[cat];

    m_catcells.swap(previous);

    return log(subtotal);

} // AlleleModel::ComputeCatDL

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

DoubleVec1d AlleleModel::ChooseAncestralState(long marker)
{
    long size = m_bincounts[marker];
    DoubleVec1d result(size, 0.0);
    result[registry.GetRandom().Long(size)] = 1.0;
    return result;
} // ChooseAncestralState

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

StepwiseModel::StepwiseModel(long nmarkers,
                             const StringVec2d& uniqueAlleles,
                             long numCategories,
                             DoubleVec1d categoryRates,
                             DoubleVec1d categoryProbabilities,
                             double userAutoCorrelationValue,
                             bool doNormalize,
                             double relMuRate)
    :   AlleleModel(nmarkers,
                    numCategories,
                    categoryRates,
                    categoryProbabilities,
                    userAutoCorrelationValue,
                    doNormalize,
                    defaults::bins,
                    relMuRate),
        m_allowance(defaults::step_allowance)
{
    Initialize(uniqueAlleles);
} // StepwiseModel constructor

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

DataModel* StepwiseModel::Clone() const
{
    DataModel* newmodel = new StepwiseModel(*this);
    return newmodel;
}

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

StringVec1d StepwiseModel::CreateDataModelReport() const
{
    StringVec1d report = DataModel::CreateDataModelReport();
    string maxbins = "Maximum number of bins:  " + ToString(m_bincounts[0]);
    report.push_back(maxbins);
    return report;
} // StepwiseModel::CreateDataModelReport

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

vector<double> StepwiseModel::DataToLikes(const string& datum, long marker) const
{
    if (m_bincounts.empty())
    {
        string msg = "StepwiseModel::DataToLikes() was called before m_bincounts ";
        msg += "was initialized.";
        throw implementation_error(msg);
    }

    if (datum == "?")       // unknown data fills all bins with 1
    {
        vector<double> result(m_bincounts[marker], 1.0);
        return result;
    }
    else
    {
        vector<double> result(m_bincounts[marker], 0.0);
        long allele;
        FromString(datum,allele);
        allele -= m_offsets[marker];
        if (allele < 0 || allele >= m_bincounts[marker])
        {
            string msg = "StepwiseModel::DataToLikes() was called on an ";
            msg += "uninitialized (or incorrectly initialized) object.";
            throw implementation_error(msg);
        }
        result[allele] = 1.0;
        return result;
    }

} // DataToLikes

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

void StepwiseModel::Initialize(const StringVec2d& uniqueAlleles)
{
    if (static_cast<unsigned long>(m_nmarkers) != uniqueAlleles.size())
    {
        string msg = "StepwiseModel::Initialize() encountered m_nmarkers = ";
        msg += ToString(m_nmarkers) + " and uniqueAlleles.size() = ";
        msg += ToString(uniqueAlleles.size()) + "; these numbers should be equal.";
        throw implementation_error(msg);
    }
    long marker;

    // find the biggest and smallest allele for each marker
    // NB We assume that allele sizes are in number of repeats,
    // *not* base pair count or anything else, and that "missing
    // data" is coded as ?.

    for (marker = 0; marker < m_nmarkers; ++marker)
    {
        bool real_allele = false;  // did we ever see a non-? allele
        long smallone = MAXLONG, largeone = 0;
        for (size_t nallele = 0; nallele<uniqueAlleles[marker].size();
             nallele++)
        {
            // convert to number
            // do not count "unknown data" markers
            string allele = uniqueAlleles[marker][nallele];
            if (allele == "?")
            {
                assert(false); //need to catch this earlier
                continue;
            }
            real_allele = true;
            long tipval;
            FromString(allele, tipval);  // convert to long

            if (tipval < smallone) smallone = tipval;
            if (tipval > largeone) largeone = tipval;
        }

        // if no non-? were ever found, use arbitrary values
        if (!real_allele)
        {
            smallone = 10;
            largeone = 10;
        }

        long newoffset = max(smallone - m_allowance, 0L);
        m_offsets.push_back(newoffset);
        m_bincounts.push_back(largeone + m_allowance + 1 - newoffset);
        if (real_allele && largeone != smallone)
            // m_threshhold: +1 makes it work for (large-small) both odd/even
            m_threshhold.push_back((largeone - smallone + 1L)/2L);
        else
            m_threshhold.push_back(1L);

        // Pre-calculate table of steps
        CalculateSteps(marker);
    }

    AlleleModel::m_nbins = *max_element(m_bincounts.begin(),
                                        m_bincounts.end());
    fill(AlleleModel::m_bincounts.begin(),
         AlleleModel::m_bincounts.end(), AlleleModel::m_nbins);
} // Initialize

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's microsatellite likelihood nuview_micro routine from Migrate.

void StepwiseModel::ComputeSiteDLs (Cell_ptr child1, Cell_ptr child2,
                                    Cell_ptr thiscell, const DoubleVec1d& lengthOfBranchToChild1ScaledByRateCat,
                                    const DoubleVec1d& lengthOfBranchToChild2ScaledByRateCat, long marker)
{
    double **pSiteDLsForChild1 = child1->GetSiteDLs(marker);
    double **pSiteDLsForChild2 = child2->GetSiteDLs(marker);
    double normForChild1 = child1->GetNorms(marker);
    double normForChild2 = child2->GetNorms(marker);

    if (!pSiteDLsForChild1 && !pSiteDLsForChild2) // they can't both be NULL
    {
        string msg = "StepwiseModel::ComputeSiteDLs() found no existing ";
        msg += "site data-likelihoods for either child.";
        throw implementation_error(msg);
    }

    // in case of a one-legged coalescence, copy values and return
    if (!pSiteDLsForChild1)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild2);
        thiscell->SetNorms(normForChild2, marker);
        return;
    }

    if (!pSiteDLsForChild2)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild1);
        thiscell->SetNorms(normForChild1, marker);
        return;
    }

    double **jointProbChild1Child2 = new double*[m_ncategories * sizeof(double *)];
    jointProbChild1Child2[0] = new double[m_ncategories * m_bincounts[marker] * sizeof(double)];
    for (long cat = 0; cat < m_ncategories; ++cat)
        jointProbChild1Child2[cat] = jointProbChild1Child2[0] + cat * m_bincounts[marker];

    long bmax = m_bincounts[marker];
    long threshold = m_threshhold[marker];
    double maxJointProbChild1Child2 = -DBL_MAX;
    double mutationProbChild1, mutationProbChild2;

    // Compute the bin contents for each possible microsat allele.
    // "b" is the "bin number," i.e., the scaled value of the starting allele
    // ("scaled" meaning "shifted so that the smallest allele gets set to zero").
    // For each starting allele "b", we sweep through a range of alleles "a" to which
    // allele "b" can mutate, from b - threshold to b + threshold.
    for (long b = 0; b < bmax; b++)
    {
        for (long cat = 0; cat < m_ncategories; ++cat)
        {
            mutationProbChild1 = mutationProbChild2 = 0.0;
            for (long a = max(0L, b - threshold); a <= min(b + threshold, bmax-1); a++)
            {
                // Note:  The probability of mutating n "steps" downward
                // (e.g., from a microsat allele of 23 repeats to an allele of 20 repeats)
                // equals the prob. of mutating n steps upward
                // (e.g., from 23 repeats to 26 repeats).
                long netNumPositiveSteps = labs(b - a);
                if (pSiteDLsForChild1[cat][a] > 0)
                {
                    mutationProbChild1 += Probability(lengthOfBranchToChild1ScaledByRateCat[cat],
                                                      netNumPositiveSteps, marker) * pSiteDLsForChild1[cat][a];
                }
                if (pSiteDLsForChild2[cat][a] > 0)
                {
                    mutationProbChild2 += Probability(lengthOfBranchToChild2ScaledByRateCat[cat],
                                                      netNumPositiveSteps, marker) * pSiteDLsForChild2[cat][a];
                }
            }
            jointProbChild1Child2[cat][b] = mutationProbChild1*mutationProbChild2;

            if (jointProbChild1Child2[cat][b] > maxJointProbChild1Child2)
            {
                maxJointProbChild1Child2 = jointProbChild1Child2[cat][b];
            }
        }
    }

    // normalize to further protect against overflow, if requested
    if (ShouldNormalize())
    {
        if (0.0 == maxJointProbChild1Child2)
        {
            thiscell->SetNorms(-DBL_MAX,marker);
        }
        else
        {
            for (long b = 0; b < bmax; b++)
            {
                for (long cat = 0; cat < m_ncategories; ++cat)
                {
                    jointProbChild1Child2[cat][b] /= maxJointProbChild1Child2;
                }
            }
            thiscell->SetNorms(log(maxJointProbChild1Child2) +
                               normForChild1 + normForChild2, marker);
        }
    }

    thiscell->SetSiteDLs(marker, jointProbChild1Child2);

    delete[] jointProbChild1Child2[0];
    delete[] jointProbChild1Child2;
} // ComputeSiteDLs

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

double StepwiseModel::ComputeSubtreeDLs(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    double total=0.0, subtotal;
    double** marker;
    long cat, bin;
    long firstposn = posn;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {
        subtotal = 0.0;
        DoubleVec1d buffer(m_ncategories, 0.0);

        for (cat = 0; cat < m_ncategories; ++cat)
        {
            for (bin = 0; bin < m_bincounts[posn]; ++bin)
            {
                // NB:  We assume a flat allele frequency prior here.
                //cout << "cat " << cat << ", bin " << bin << ", prob " << marker[cat][bin] << endl;
                buffer[cat] += marker[cat][bin];
            }
            subtotal += m_catprobs[cat] * buffer[cat];
        }

        if (!subtotal)
        {
            DataModel::TryToNormalizeAndThrow(posn, GetModelType());
        }

        total += (log(subtotal) + rootdls.GetNorms(posn));

        assert (total != 0);  // that would be *too* likely

        if (m_ncategories > 1)
        {
            for (cat = 0; cat < m_ncategories; cat++)
                m_likes[posn][cat] = buffer[cat]/subtotal;
        }
        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);
    return total;

} // StepwiseModel::ComputeSubtreeDLs

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

// This method returns the probability of changing an allele by "diff" net steps
// in time t.  For example, if the data type is microsatellite, we could compute
// the probability that the sequence ACACAC (3 repeats) mutates to ACACACACAC
// (5 repeats, for a "diff" of 2) along a branch whose length is scaled to be t.
// "diff" received by this method is nonnegative, representing a net increase
// in the size of the allele; the probability of decreasing by "diff" net steps
// is defined to be the same as increasing by "diff" net steps.
// All nonvanishingly unlikely net steps are considered; e.g., simply taking 2
// steps along the branch, or taking 2 steps forward plus k steps backward
// plus k steps forward, where k = 0, 1, 2, 3, ..., infinity.
// The formula is Prob(i net steps in time t) =
//    exp(-t)*sum_over_k((t/2)^(i+2k) / ((i+k)!k!)),
// where the sum runs from k = 0 to k = infinity.
// This formula is eq. 15.26 in Joe's "Inferring Phylogenies" book (p. 242).
// The code looks different because the factorials are obtained via a precomputed
// lookup table, and exp/log are used to counteract underflow.
// Contributions to the sum drop off rapidly as k increases.
// (Side note:  An equivalent and more compact version of this formula is
// Prob(i,t) = exp(-t)*BesselI(i,t), where BesselI(n,x) is the
// modified Bessel function of integer order n evaluated at real argument x.)

double StepwiseModel::Probability(double t, long diff, long marker) const
{
    long threshold = m_threshhold[marker]; // max. number of positive steps
    if (diff > threshold)
        return 0.0; // approximately infinitely unlikely to mutate that much in time t

    double sum(0.0), oldsum(0.0);
    const DoubleVec2d& PrecomputedTerms = m_steps[marker];
    double log_tOver2 = log(0.5 * t);

    for (long k = 0; k <= threshold; k++) // num steps = diff + k <= threshold
    {
        sum += exp(-t + log_tOver2*(diff + 2.0*k) - PrecomputedTerms[diff][k]);

        // quit if the contributions have become trivial
        if (fabs (oldsum - sum) < DBL_EPSILON) break;
        oldsum = sum;
    }

    return sum;
}

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's calculate_steps routine from MIGRATE.
// This routine precomputes values needed by Probability, for speed.

void StepwiseModel::CalculateSteps(long marker)
{
    long k, diff;
    DoubleVec1d tempvec;
    DoubleVec2d steps;
    long threshhold = m_threshhold[marker];

    for (diff = 0; diff <= threshhold; diff++)
    {
        tempvec.clear();
        for (k = 0; k <= threshhold; k++)
        {
            tempvec.push_back(logfac (diff + k) + logfac (k));
        }
        steps.push_back(tempvec);
    }

    m_steps.push_back(steps); // Note:  This is okay from a speed standpoint,
    // but not so good from a design standpoint,
    // because this method assumes it's being called
    // within a loop over markers.
} // calculate_steps

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

DoubleVec1d StepwiseModel::SimulateMarker(double branchlength, long whichmarker, const DoubleVec1d& state) const
{
    throw implementation_error("Cannot simulate data with the stepwise model yet.");
} // SimulateMarker

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

string StepwiseModel::CellToData(Cell_ptr cell, long marker) const
{
    LongVec1d ones = cell->GetOnes(marker);
    assert(static_cast<size_t>(marker) < m_offsets.size());
    if (ones.size() == 1)
    {
        return ToString(ones[0] + m_offsets[marker]);
    }
    throw implementation_error
        ("Cannot convert stepwise data from the internal format for data not simply a single number.");
}

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

BrownianModel::BrownianModel(long nmarkers,
                             long numCategories,
                             DoubleVec1d categoryRates,
                             DoubleVec1d categoryProbabilities,
                             double userAutoCorrelationValue,
                             bool doNormalize,
                             double relMuRate)
    : AlleleModel(nmarkers,
                  numCategories,
                  categoryRates,
                  categoryProbabilities,
                  userAutoCorrelationValue,
                  doNormalize,
                  defaults::brownianBins, // 3, for mean, variance, and cumulative total
                  relMuRate)
{
    // intentionally blank
} // BrownianModel constructor

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

DataModel* BrownianModel::Clone() const
{
    DataModel* newmodel = new BrownianModel(*this);
    return newmodel;
}

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

void BrownianModel::SetNormalize(bool norm)
{
    if(norm)
    {
        data_error e("Normalization cannot be set for Brownian Model");
        throw e;
    }
} // BrownianModel::SetNormalize

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

vector<double> BrownianModel::DataToLikes(const string& datum, long) const
{
    vector<double> result(m_nbins,0.0);

    if (datum == "?")                   // unknown data
    {
        result[0] = 5.0;                // this is an arbitrary value
        result[1] = DBL_BIG;
        result[2] = 0.0;
    }
    else
    {
        FromString(datum,result[0]);
        result[1] = 0.0;
        result[2] = 0.0;
    }

    return result;

} // DataToLikes

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's nuview_brownian() from Migrate-1.2.4
// by Jon Yamato 2002/05/06

// N[Log[1/Sqrt[2 Pi]], 30]
#define LOG2PIHALF -0.918938533204672741780329736406

void BrownianModel::ComputeSiteDLs(Cell_ptr child1, Cell_ptr child2,
                                   Cell_ptr thiscell, const DoubleVec1d& vv1, const DoubleVec1d& vv2,
                                   long marker)
{
    double mean1, mean2, xx1, xx2, c12, v1, v2, vtot, f1, f2;
    double **c1dls = child1->GetSiteDLs(marker),
        **c2dls = child2->GetSiteDLs(marker);

    if (!c1dls && !c2dls)
    {
        string msg = "BrownianModel::ComputeSiteDLs() failed to find ";
        msg += "data likelihoods for either child.";
        throw implementation_error(msg);
    }

    if (!c1dls)
    {
        thiscell->SetSiteDLs(marker,c2dls);
        return;
    }

    if (!c2dls)
    {
        thiscell->SetSiteDLs(marker,c1dls);
        return;
    }

    // temporary space needed for interface with dlcell::SetSiteDLs()
    double **mydls = new double*[m_ncategories*sizeof(double*)];
    mydls[0] = new double[m_ncategories*m_nbins*sizeof(double)];
    long cat;
    for(cat = 1; cat < m_ncategories; ++cat)
        mydls[cat] = mydls[0] + cat*m_nbins;

    for (cat = 0; cat < m_ncategories; ++cat)
    {

        mean1 = c1dls[cat][0];
        xx1 = c1dls[cat][2];
        v1 = vv1[cat] + c1dls[cat][1];

        mean2 = c2dls[cat][0];
        xx2 = c2dls[cat][2];
        v2 = vv2[cat] + c2dls[cat][1];

        vtot = v1 + v2;

        // the weights are set reciprocally so that the value coming from
        // the shorter branch is given more weight.
        if (vtot > 0.0) f1 = v2/vtot;
        else f1 = 0.5;
        f2 = 1.0 - f1;

        mydls[cat][0] = f1*mean1 + f2*mean2;

        mydls[cat][1] = v1*f1;

        mydls[cat][2] = xx1 + xx2;

        c12 = (mean1-mean2)*(mean1-mean2) / vtot;
        mydls[cat][2] += min(0.0,-0.5 * (log(vtot)+c12) + LOG2PIHALF);
    }

    thiscell->SetSiteDLs(marker,mydls);

    delete [] mydls[0];
    delete [] mydls;

} // BrownianModel::ComputeSiteDLs

#undef LOG2PIHALF

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

double BrownianModel::ComputeSubtreeDLs(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    // NB:  Brownian likelihoods are stored as logs!

    double total=0.0, subtotal;
    double** marker;
    long cat;
    long firstposn = posn;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {

        if (m_ncategories > 1)
        {
            // in order to add up likelihoods over categories, we
            // must un-log them.  We normalize them first to avoid
            // underflow.

            DoubleVec1d buffer(m_ncategories, 0.0);
            subtotal = 0.0;

            double biggest = NEGMAX;

            for (cat = 0; cat < m_ncategories; ++cat)
            {
                if (marker[cat][2] > biggest) biggest = marker[cat][2];
            }

            for (cat = 0; cat < m_ncategories; ++cat)
            {
                buffer[cat] += exp(marker[cat][2] - biggest);
                subtotal += m_catprobs[cat] * buffer[cat];
            }

            for (cat = 0; cat < m_ncategories; cat++)
                m_likes[posn][cat] = buffer[cat]/subtotal;

            if (!subtotal)
            {
                // we shouldn't be here!  Normalization shouldn't happen
                // for Brownian.
                datalike_error ex("invalid subtree found in Brownian model");
                throw(ex);
            }

            total += (log(subtotal) + biggest);

        }
        else
        {
            // If there is only one category, there is no normalization
            // and we MUST NOT un-log the likelihood or it will underflow.
            total += marker[0][2];
        }

        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);

    return total;

} // BrownianModel::ComputeSubtreeDLs

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

StringVec1d BrownianModel::CreateDataModelReport() const
{
    StringVec1d report = DataModel::CreateDataModelReport();
    string rptline("(The brownian approximation for microsatellite ");
    rptline += "evolution has no extra parameters.)";
    report.push_back(rptline);

    return report;

} // BrownianModel::CreateDataModelReport

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

DoubleVec1d BrownianModel::SimulateMarker(double branchlength, long whichmarker, const DoubleVec1d& state) const
{
    throw implementation_error("Cannot simulate data with the Brownian model yet.");
} // SimulateMarker

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

string BrownianModel::CellToData(Cell_ptr cell, long marker) const
{
    double** dls = cell->GetSiteDLs(marker);
    if (dls != NULL)
    {
        return ToString(dls[0][0]);
    }
    return "0";
}

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

KAlleleModel::KAlleleModel(
    long nmarkers,
    const StringVec2d& uniqueAlleles,
    long numCategories,
    DoubleVec1d categoryRates,
    DoubleVec1d categoryProbabilities,
    double userAutoCorrelationValue,
    bool doNormalize,
    double relMuRate
    )
    : AlleleModel(
        nmarkers,
        numCategories,
        categoryRates,
        categoryProbabilities,
        userAutoCorrelationValue,
        doNormalize,
        defaults::bins,
        relMuRate
        )
{
    Initialize(uniqueAlleles);
} // KAlleleModel constructor

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

void KAlleleModel::ComputeSiteDLs(Cell_ptr child1, Cell_ptr child2,
                                  Cell_ptr thiscell, const DoubleVec1d& lengthOfBranchToChild1ScaledByRateCat,
                                  const DoubleVec1d& lengthOfBranchToChild2ScaledByRateCat, long marker)
{
    double **pSiteDLsForChild1 = child1->GetSiteDLs(marker);
    double **pSiteDLsForChild2 = child2->GetSiteDLs(marker);
    double normForChild1 = child1->GetNorms(marker);
    double normForChild2 = child2->GetNorms(marker);

    if (!pSiteDLsForChild1 && !pSiteDLsForChild2)
    {
        string msg = "KAlleleModel::ComputeSiteDLs() found no existing ";
        msg += "site data-likelihoods for either child.";
        throw implementation_error(msg);
    }

    // in case of a one-legged coalescence, copy values and return
    if (!pSiteDLsForChild1)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild2);
        thiscell->SetNorms(normForChild2, marker);
        return;
    }

    if (!pSiteDLsForChild2)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild1);
        thiscell->SetNorms(normForChild1, marker);
        return;
    }

    long smax = m_bincounts[marker];

    // allocate temporary working space;
    // OPTIMIZE should be done more cheaply!
    double **jointProbChild1Child2 = new double*[m_ncategories * sizeof(double*)];
    jointProbChild1Child2[0] = new double[m_ncategories * m_bincounts[marker] * sizeof(double)];
    for (long cat = 0; cat < m_ncategories; ++cat)
        jointProbChild1Child2[cat] = jointProbChild1Child2[0] + cat * m_bincounts[marker];

    // compute the bin contents for each possible allele size
    long s, ss;
    double prob1, prob2, a1, a2, b1, b2, temp1, temp2, sum1, sum2;
    double maxJointProbChild1Child2(-DBL_MAX);

    for (long cat = 0; cat < m_ncategories; ++cat)
    {
        prob1 = probMathFunc ( lengthOfBranchToChild1ScaledByRateCat[cat], (smax-1.0)/smax);
        prob2 = probMathFunc ( lengthOfBranchToChild2ScaledByRateCat[cat], (smax-1.0)/smax);

        a1 = 1.0 - prob1;
        a2 = 1.0 - prob2;
        b1 = prob1 / ( smax - 1.0); // smax was set to at least 2 in Initialize()
        b2 = prob2 / ( smax - 1.0);
        for (s = 0; s < smax; s++)
        {
            sum1 = 0; sum2 = 0;
            for(ss = 0; ss < smax; ss++)
            {
                if (s == ss)
                {
                    temp1 = a1; temp2 = a2;
                }
                else
                {
                    temp1 = b1; temp2 = b2;
                }
                sum1 += pSiteDLsForChild1[cat][ss] * temp1;
                sum2 += pSiteDLsForChild2[cat][ss] * temp2;
            }

            jointProbChild1Child2[cat][s] = sum1 * sum2;
            if (jointProbChild1Child2[cat][s] > maxJointProbChild1Child2)
                maxJointProbChild1Child2 = jointProbChild1Child2[cat][s];  // overflow protection
        }
    }

    // normalize to further protect against overflow, if requested
    if (ShouldNormalize())
    {
        if (maxJointProbChild1Child2 == 0.0)
        {
            thiscell->SetNorms(-DBL_MAX,marker);
        }
        else
        {
            for (s = 0; s < smax; s++)
            {
                for (long cat = 0; cat < m_ncategories; ++cat)
                {
                    jointProbChild1Child2[cat][s] /= maxJointProbChild1Child2;
                }
            }
            thiscell->SetNorms(log(maxJointProbChild1Child2) + normForChild1 + normForChild2, marker);
        }
    }

    thiscell->SetSiteDLs(marker, jointProbChild1Child2);

    delete[] jointProbChild1Child2[0];
    delete[] jointProbChild1Child2;
} // KAlleleModel::ComputeSiteDLs

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

DataModel* KAlleleModel::Clone() const
{
    DataModel* newmodel = new KAlleleModel(*this);
    return newmodel;
} // KAllelModel::Clone

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

double KAlleleModel::probMathFunc( double ut, double coef)
{
    double mathValue = coef*(1.0-exp(-1.0/coef*ut));
    if (systemSpecificIsnan(mathValue)) // BUGBUG Presumably some other action should be taken?
        cout<<"coef:"<<coef<<endl
            <<"time:"<<ut<<endl;
    return mathValue;
} // KAllelModel::probMathFunc

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

StringVec1d KAlleleModel::CreateDataModelReport() const
{
    StringVec1d report = DataModel::CreateDataModelReport();
    string maxbins = "Maximum number of bins:  " + ToString(m_bincounts[0]);
    report.push_back(maxbins);
    return report;

} // KAllelModel::CreateDataModelReport

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

void KAlleleModel::Initialize(const StringVec2d& uniqueAlleles)
{
    if (static_cast<unsigned long>(m_nmarkers) != uniqueAlleles.size())
    {
        string msg = "KAlleleModel::Initialize() encountered m_nmarkers = ";
        msg += ToString(m_nmarkers) + " and uniqueAlleles.size() = ";
        msg += ToString(uniqueAlleles.size()) + "; these numbers should be equal.";
        throw implementation_error(msg);
    }

    for (long marker = 0; marker < m_nmarkers; ++marker)
    {
        map <string,long> tmpMap;
        m_allelemaps.push_back(tmpMap);
        size_t nallele = 0;
        for (; nallele<uniqueAlleles[marker].size(); nallele++)
        {
            string allelename = uniqueAlleles[marker][nallele];
            if (allelename == "?")
            {
                assert(false); //we need to catch this earlier.
                continue;
            }
            assert(m_allelemaps[marker].find(allelename)==m_allelemaps[marker].end());
            m_allelemaps[marker].insert(make_pair(allelename, nallele));
        }

        if(nallele <= 1) nallele=2; // must be at least this much
        m_bincounts.push_back(nallele);
    }

    AlleleModel::m_nbins = *max_element(m_bincounts.begin(),
                                        m_bincounts.end());
    fill(AlleleModel::m_bincounts.begin(),
         AlleleModel::m_bincounts.end(), AlleleModel::m_nbins);
} // KAlleleModel::Initialize

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

vector<double> KAlleleModel::DataToLikes(const string& datum, long marker) const
{
    if (m_bincounts.empty())
    {
        string msg = "KAlleleModel::DataToLikes() was called on an uninitialized ";
        msg += "object.";
        throw implementation_error(msg);
    }

    if (datum == "?")       // unknown data fills all bins with 1
    {
        vector<double> result(m_bincounts[marker], 1.0);
        return result;
    }
    else
    {
        vector<double> result(m_bincounts[marker], 0.0);
        map<string, long>::const_iterator allele = m_allelemaps[marker].find(datum);
        assert(allele != m_allelemaps[marker].end()); //Uncounted allele
        result[allele->second] = 1.0;
        return result;
    }

} // DataToLikes

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

double KAlleleModel::ComputeSubtreeDLs(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    double total=0.0, subtotal;
    double** marker;
    long cat, bin;
    long firstposn = posn;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {
        subtotal = 0.0;
        DoubleVec1d buffer(m_ncategories, 0.0);

        for (cat = 0; cat < m_ncategories; ++cat)
        {
            for (bin = 0; bin < m_bincounts[posn]; ++bin)
            {
                // NB:  We assume a flat allele frequency prior here.
                buffer[cat] += marker[cat][bin]/m_nbins;
                //LS DEBUG MAPPING:  we eventually want unique allele freqs here instead
                // of dividing by m_nbins.
            }

            subtotal += m_catprobs[cat] * buffer[cat];
            //cout<<m_catprobs[cat]<<endl<<buffer[cat]<<endl;
            //assert(subtotal != 0);  // that would not be likely enough
            if (!subtotal)
            {
                DataModel::TryToNormalizeAndThrow(posn, GetModelType());
            }
        }

        total += (log(subtotal) + rootdls.GetNorms(posn));

        assert (total != 0);  // that would be *too* likely

        if (m_ncategories > 1)
        {
            for (cat = 0; cat < m_ncategories; cat++)
                m_likes[posn][cat] = buffer[cat]/subtotal;
        }
        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);

    if(systemSpecificIsnan(total))
    {
        cout<<endl<<"ComputeSubtreeDLs";
        cout<<endl;
    }

    return total;

}; // KAlleleModel::ComputeSubtreeDLs

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

StringVec1d KAlleleModel::ToXML(size_t nspaces) const
{
    StringVec1d xmllines(DataModel::ToXML(nspaces));

    // Yes, this model has no fields of its own.

    return xmllines;

} // ToXML

//------------------------------------------------------------------------------------
// This routine simulates data for a single node under the KAllele
// model.  It assumes that the given branch lengths have already
// been rescaled for the rate category desired.

DoubleVec1d KAlleleModel::SimulateMarker(double branchlength, long whichmarker, const DoubleVec1d& state) const
{
    double rnd = registry.GetRandom().Float();
    long nstates(m_bincounts[whichmarker]);
    long oldstate = 0;
    for (long staten=0; staten<static_cast<long>(state.size()); staten++)
    {
        if (state[staten] == 1.0) oldstate = staten;
    }

    // if something happens
    double newratio = nstates/(nstates-1);
    if (rnd < (1/newratio) * (1.0 - exp(-newratio * branchlength)))
    {
        long chosenstate = registry.GetRandom().Long(nstates-1);
        if (chosenstate >= oldstate) chosenstate++;
        DoubleVec1d answer(state.size(),0.0);
        answer[chosenstate] = 1.0;
        return answer;
    }

    // return nothing happens
    return state;

} // SimulateMarker

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

string KAlleleModel::CellToData(Cell_ptr cell, long marker) const
{
    LongVec1d ones = cell->GetOnes(marker);
    if (ones.size() == 1)
    {
        for (map<string, long>::const_iterator key=m_allelemaps[marker].begin();
             key != m_allelemaps[marker].end(); key++)
        {
            if ((*key).second == ones[0])
            {
                return (*key).first;
            }
        }
        throw data_error("Cannot find any any alleles for bin " + ToString(ones[0]) + ".");
    }
    throw implementation_error
        ("Cannot convert K-Allele data from the internal format if the data is not a single allele.");

}

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

MixedKSModel::MixedKSModel(
    long nmarkers,
    const StringVec2d& uniqueAlleles,
    long numCategories,
    DoubleVec1d categoryRates,
    DoubleVec1d categoryProbabilities,
    double userAutoCorrelationValue,
    bool doNormalize,
    double relMuRate,
    double alphaVal,
    bool doOptimization
    )
    : AlleleModel(
        nmarkers,
        numCategories,
        categoryRates,
        categoryProbabilities,
        userAutoCorrelationValue,
        doNormalize,
        defaults::bins,
        relMuRate
        ),
      m_allowance(defaults::mixedks_allowance),
      m_isOpt(doOptimization),
      m_alphaRpt()
{
    SetAlpha(alphaVal); // instead of initializing, because there's
                        // logic to check value of alpha in SetAlpha
    m_origAlpha = m_alpha;
    Initialize(uniqueAlleles);
}

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

DataModel* MixedKSModel::Clone() const
{
    DataModel* newmodel = new MixedKSModel(*this);
    return newmodel;
}

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

vector<double> MixedKSModel::DataToLikes(const string& datum, long marker) const
{
    if (m_bincounts.empty())
    {
        string msg = "MixedKSModel::DataToLikes() was called on an uninitialized ";
        msg += "object.";
        throw implementation_error(msg);
    }

    if (datum == "?")                   // unknown data fills all bins with 1
    {
        vector<double> result(m_bincounts[marker], 1.0);
        return result;
    }
    else if (StringType(datum)!=2)
    {
        throw data_error ("Your data must consist of numbers, not letters or punctuation."
                          "  The Mixed KS model is inappropriate for DNA or RNA data.");
    }
    else
    {
        vector<double> result(m_bincounts[marker], 0.0);
        long allele;
        FromString(datum,allele);
        allele -= m_offsets[marker];
        if (allele < 0 || allele >= m_bincounts[marker])
        {
            string msg = "MixedKSModel::DataToLikes() was called on an ";
            msg += "uninitialized (or incorrectly initialized) object.";
            throw implementation_error(msg);
        }
        result[allele] = 1.0;
        return result;
    }

} // DataToLikes

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

void MixedKSModel::Initialize(const StringVec2d& uniqueAlleles)
{
    if (static_cast<unsigned long>(m_nmarkers) != uniqueAlleles.size())
    {
        string msg = "MixedKSModel::Initialize() encountered m_nmarkers = ";
        msg += ToString(m_nmarkers) + " and uniqueAlleles.size() = ";
        msg += ToString(uniqueAlleles.size()) + "; these numbers should be equal.";
        throw implementation_error(msg);
    }

    // find the biggest and smallest allele for each marker
    // NB We assume that allele sizes are in number of repeats,
    // *not* base pair count or anything else, and that "missing
    // data" is coded as ?.

    long tipval;

    for (long marker = 0; marker < m_nmarkers; ++marker)
    {
        bool real_allele = false;  // did we ever see a non-? allele
        long smallone = MAXLONG, largeone = 0;

        for (size_t nallele = 0; nallele<uniqueAlleles[marker].size();
             nallele++)
        {
            // convert to number if possible;
            // do not count "unknown data" markers
            string allele = uniqueAlleles[marker][nallele];
            if (allele == "?")
            {
                assert(false); //catch this earlier
                continue;
            }
            real_allele = true;
            FromString(allele, tipval);  // convert to long

            if (tipval < smallone) smallone = tipval;
            if (tipval > largeone) largeone = tipval;
        }

        // if no non-? were ever found, use arbitrary values
        if (!real_allele)
        {
            smallone = 10;
            largeone = 10;
        }

        long newoffset = max(smallone - m_allowance, 0L);
        m_offsets.push_back(newoffset);
        m_bincounts.push_back(largeone + m_allowance + 1 - newoffset);
        if (real_allele && largeone != smallone)
            // m_threshhold: +1 makes it work for (large-small) both odd/even
            m_threshhold.push_back((largeone - smallone + 1L)/2L);
        else
            m_threshhold.push_back(1L);

        // Pre-calculate table of steps
        CalculateSteps(marker);
    }

    AlleleModel::m_nbins = *max_element(m_bincounts.begin(),
                                        m_bincounts.end());
    fill(AlleleModel::m_bincounts.begin(),
         AlleleModel::m_bincounts.end(), AlleleModel::m_nbins);
} // Initialize

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

void MixedKSModel::SetAlpha(double alphaVal)
{
    assert(alphaVal >= 0.0);
    assert(alphaVal <= 1.0);
    m_alpha = alphaVal;
    m_beta = 1 - m_alpha;
}

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

void MixedKSModel::SetAlpha(double alphaVal, long rep, long chain)
{
    if (m_alphaRpt.size() == static_cast<size_t>(rep))
    {
        ResetAlpha();
    }
    assert(m_alphaRpt[rep].size() == static_cast<size_t>(chain));
    m_alphaRpt[rep].push_back(alphaVal);
    SetAlpha(alphaVal);
}

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's microsatellite likelihood nuview_micro routine from Migrate.

void MixedKSModel::ComputeSiteDLs (Cell_ptr child1, Cell_ptr child2,
                                   Cell_ptr thiscell, const DoubleVec1d& lengthOfBranchToChild1ScaledByRateCat,
                                   const DoubleVec1d& lengthOfBranchToChild2ScaledByRateCat, long marker)
{
    double **pSiteDLsForChild1 = child1->GetSiteDLs(marker);
    double **pSiteDLsForChild2 = child2->GetSiteDLs(marker);
    double normForChild1 = child1->GetNorms(marker);
    double normForChild2 = child2->GetNorms(marker);

    if (!pSiteDLsForChild1 && !pSiteDLsForChild2) // they can't both be NULL
    {
        string msg = "MixedKSModel::ComputeSiteDLs() found no existing ";
        msg += "site data-likelihoods for either child.";
        throw implementation_error(msg);
    }

    // in case of a one-legged coalescence, copy values and return
    if (!pSiteDLsForChild1)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild2);
        thiscell->SetNorms(normForChild2, marker);
        return;
    }

    if (!pSiteDLsForChild2)
    {
        thiscell->SetSiteDLs(marker, pSiteDLsForChild1);
        thiscell->SetNorms(normForChild1, marker);
        return;
    }

    long smax = m_bincounts[marker];
    long threshold = m_threshhold[marker];

    // allocate temporary working space;
    // OPTIMIZE should be done more cheaply!
    double **jointProbChild1Child2 = new double*[m_ncategories * sizeof(double *)];
    jointProbChild1Child2[0] = new double[m_ncategories * m_bincounts[marker] * sizeof(double)];
    for (long cat = 0; cat < m_ncategories; ++cat)
        jointProbChild1Child2[cat] = jointProbChild1Child2[0] + cat * m_bincounts[marker];
    double maxJointProbChild1Child2 = -DBL_MAX;
    double mutationProbChild1, mutationProbChild2;

    long s, a, diff;

    //KAllele part
    // compute the bin contents for each possible allele size
    long ss;
    double prob1, prob2, a1, a2, b1, b2, temp1, temp2, sum1, sum2;

    for (long cat = 0; cat < m_ncategories; ++cat)
    {
        prob1 = probMathFunc (lengthOfBranchToChild1ScaledByRateCat[cat], (smax-1.0)/smax);
        prob2 = probMathFunc (lengthOfBranchToChild2ScaledByRateCat[cat], (smax-1.0)/smax);

        a1 = 1.0 - prob1;
        a2 = 1.0 - prob2;
        b1 = prob1 / ( smax - 1.0); //need to check when smax = 1
        b2 = prob2 / ( smax - 1.0);
        for (s = 0; s < smax; s++)
        {
            sum1 = 0; sum2 = 0;
            for(ss = 0; ss < smax; ss++)
            {
                if (s == ss)
                {
                    temp1 = a1; temp2 = a2;
                }
                else
                {
                    temp1 = b1; temp2 = b2;
                }
                sum1 += pSiteDLsForChild1[cat][ss] * temp1;
                sum2 += pSiteDLsForChild2[cat][ss] * temp2;
            }

            jointProbChild1Child2[cat][s] = sum1 * sum2 * m_alpha;  //alpha
            if (jointProbChild1Child2[cat][s] > maxJointProbChild1Child2)
                maxJointProbChild1Child2 = jointProbChild1Child2[cat][s];  // overflow protection
        }
    }
    //LS NOTE:  The following code firing was due to a bug in the past;
    // you should probably look carefully to make sure that's not true
    // again.
    if(maxJointProbChild1Child2 < 0)
    {
        string msg = "Warning:  maximum probability value of "
            + ToString(maxJointProbChild1Child2)
            + " in MixedKS.";
        registry.GetRunReport().ReportDebug(msg);
    }

    //StepWise part
    // compute the bin contents for each possible allele size
    for (s = 0; s < smax; s++)
    {
        for (long cat = 0; cat < m_ncategories; ++cat)
        {
            mutationProbChild1 = mutationProbChild2 = 0.0;
            for (a = max (0L, s - threshold); a <= s + threshold && a < smax; a++)
            {
                diff = labs(s - a);
                if (pSiteDLsForChild1[cat][a] > 0)
                {
                    mutationProbChild1 += Probability(lengthOfBranchToChild1ScaledByRateCat[cat],
                                                      diff, marker) * pSiteDLsForChild1[cat][a];
                }
                if (pSiteDLsForChild2[cat][a] > 0)
                {
                    mutationProbChild2 += Probability(lengthOfBranchToChild2ScaledByRateCat[cat],
                                                      diff, marker) * pSiteDLsForChild2[cat][a];
                }
            }
            jointProbChild1Child2[cat][s] += mutationProbChild1 * mutationProbChild2 * m_beta;
            if (jointProbChild1Child2[cat][s] > maxJointProbChild1Child2)
                maxJointProbChild1Child2 = jointProbChild1Child2[cat][s];  // overflow protection
        }
    }

    // normalize to further protect against overflow, if requested
    if (ShouldNormalize())
    {
        if (maxJointProbChild1Child2 == 0.0)
        {
            thiscell->SetNorms(-DBL_MAX,marker);
        }
        else
        {
            for (s = 0; s < smax; s++)
            {
                for (long cat = 0; cat < m_ncategories; ++cat)
                {
                    jointProbChild1Child2[cat][s] /= maxJointProbChild1Child2;
                }
            }
            thiscell->SetNorms(log(maxJointProbChild1Child2) + normForChild1 + normForChild2,
                               marker);
        }
    }

    thiscell->SetSiteDLs(marker, jointProbChild1Child2);

    delete[] jointProbChild1Child2[0];
    delete[] jointProbChild1Child2;
} // ComputeSiteDLs

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

double MixedKSModel::ComputeSubtreeDLs(Cell& rootdls, double** startmarker, double** endmarker, long posn)
{
    double total=0.0, subtotal;
    double** marker;
    long cat, bin;
    long firstposn = posn;

    for (marker = startmarker; marker != endmarker;
         marker = rootdls.GetNextMarker(marker))
    {
        subtotal = 0.0;
        DoubleVec1d buffer(m_ncategories, 0.0);

        for (cat = 0; cat < m_ncategories; ++cat)
        {
            for (bin = 0; bin < m_bincounts[posn]; ++bin)
            {
                // NB:  We assume a flat allele frequency prior here.
                buffer[cat] += marker[cat][bin];
            }

            subtotal += m_catprobs[cat] * buffer[cat];
        }

        if (!subtotal)
        {
            DataModel::TryToNormalizeAndThrow(posn, GetModelType());
        }

        total += (log(subtotal) + rootdls.GetNorms(posn));

        assert (total != 0);  // that would be *too* likely

        if (m_ncategories > 1)
        {
            for (cat = 0; cat < m_ncategories; cat++)
                m_likes[posn][cat] = buffer[cat]/subtotal;
        }
        ++posn;
    }

    if (m_ncategories > 1) total += ComputeCatDL(firstposn, posn);
    return total;

} // MixedKSModel::ComputeSubtreeDLs

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

StringVec1d MixedKSModel::ToXML(size_t nspaces) const
{
    StringVec1d xmllines(AlleleModel::ToXML(nspaces));

    nspaces += INDENT_DEPTH;

    string line(MakeIndent(MakeTag(xmlstr::XML_TAG_ISOPT),nspaces));
    line += ToStringTF(m_isOpt);
    line += MakeCloseTag(xmlstr::XML_TAG_ISOPT);
    StringVec1d::iterator endtag = --xmllines.end();
    xmllines.insert(endtag,line);

    string line2(MakeIndent(MakeTag(xmlstr::XML_TAG_ALPHA),nspaces));
    line2 += ToString(GetAlpha());
    line2 += MakeCloseTag(xmlstr::XML_TAG_ALPHA);
    endtag = --xmllines.end();
    xmllines.insert(endtag,line2);

    nspaces -= INDENT_DEPTH;
    return xmllines;

} // ToXML

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's prob_micro routine from MIGRATE.
// This routine computes the probability of a change of "diff" steps in time "t".
//
//   Mary Kuhner 2002/01/02

double MixedKSModel::Probability (double t, long diff, long marker) const
{
    long threshold = m_threshhold[marker];
    if (diff > threshold)
        return 0.0; // approximately infinitely unlikely to mutate that much in time t

    const DoubleVec2d& PrecomputedTerms = m_steps[marker];
    double sum(0.0), oldsum(0.0), log_tOver2(log(0.5*t));

    for (long k = 0; k <= threshold; k++) // num steps = diff + k <= threshold
    {
        sum += exp(-t + log_tOver2*(diff + 2.0*k) - PrecomputedTerms[diff][k]);

        // quit if the contributions have become trivial
        if (fabs (oldsum - sum) < DBL_EPSILON) break;
        oldsum = sum;
    }
    return sum;
} // Probability

//------------------------------------------------------------------------------------
// Adaptation of Peter Beerli's calculate_steps routine from MIGRATE.
// This routine precomputes values needed by Probability, for speed.

void MixedKSModel::CalculateSteps (long marker)
{
    long k, diff;
    DoubleVec1d tempvec;
    DoubleVec2d steps;
    long threshhold = m_threshhold[marker];

    for (diff = 0; diff <= threshhold; diff++)
    {
        tempvec.clear();
        for (k = 0; k <= threshhold; k++)
        {
            tempvec.push_back(logfac (diff + k) + logfac (k));
        }
        steps.push_back(tempvec);
    }

    m_steps.push_back(steps);
    // Note:  This is okay from a speed standpoint,
    // but not so good from a design standpoint,
    // because this method assumes it's being called
    // within a loop over markers.
} // CalculateSteps

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

double MixedKSModel::probMathFunc( double ut, double coef)
{
    double mathValue = coef*(1.0-exp(-1.0/coef*ut));
    if (systemSpecificIsnan(mathValue))
        cout << "coef:" << coef << endl
             << "time:" << ut << endl;
    return mathValue;
} // MixedKSModel::probMathFunc

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

StringVec1d MixedKSModel::CreateDataModelReport() const
{
    StringVec1d report = DataModel::CreateDataModelReport();

    string line;
    if (m_isOpt)
    {
        line = "Final ";
    }
    line += "Multistep:single-step ratio: " + ToString(m_alpha);
    report.push_back(line);

    string maxbins = "Maximum number of bins:  " + ToString(m_bincounts[0]);
    report.push_back(maxbins);
    if (m_isOpt)
    {
        report.push_back(string("Multistep:single-step ratio used for each chain (optimized from the previous):"));
        //report alpha for each chain
        for(size_t rep=0; rep<m_alphaRpt.size(); rep++)
        {
            for (size_t chain=0; chain<m_alphaRpt[rep].size(); chain++)
            {
                line = "";
                if (m_alphaRpt.size() > 1)
                {
                    line = "Replicate " + indexToKey(rep) + ", ";
                }
                line += "Chain " + indexToKey(chain) + " :"
                    + ToString(m_alphaRpt[rep][chain]);
                report.push_back( line );
            }
        }
    }

    return report;

} // MixedKSModel::CreateDataModelReport

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

bool MixedKSModel::OptimizeDataModel(Tree* tree, const Locus& locus)
{
    //only do optimization when turn on 'm_isOpt', otherwize do nothing
    if(!m_isOpt) return false;

    FuncMax fm(*this, tree, locus);

    FMState fms(100,      //max_iters
                0.01,     //increment
                0.0001,   //threshold
                0,        //limits
                1,
                0.5       //initX
        );

    fm.setState(fms);
    fm.run();
    //Note:  The FuncMax object has a pointer to us, and uses it to set our alpha.
    //report final alpha for each chain
    m_alphaRpt[m_alphaRpt.size()-1].push_back(m_alpha);
    return true;
}

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

void MixedKSModel::ResetAlpha()
{
    SetAlpha(m_origAlpha);
    DoubleVec1d newalphas;
    newalphas.push_back(m_origAlpha);
    m_alphaRpt.push_back(newalphas);
}

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

void MixedKSModel::WriteAlpha(ofstream& sumout, long loc, long rep, long chain)
{
    if (!m_isOpt) return;
    assert(static_cast<size_t>(rep)<m_alphaRpt.size());
    if (static_cast<size_t>(chain+1)>=m_alphaRpt[rep].size()) return;
    //LS NOTE:  We  write out the alpha that we calculated for this chain
    // instead of the alpha we used for this chain, because of timing
    // issues--we don't want to lose the information if sumfile writing stops
    // unexpectedly.
    sumout << "\t"
           << xmlsum::ALPHA_START1 << " "
           << xmlsum::ALPHA_START2 << " " << loc << " "
           << xmlsum::ALPHA_START3 << " " << m_alphaRpt[rep][chain+1]
           << " " << xmlsum::ALPHA_END << endl;

}

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

DoubleVec1d MixedKSModel::SimulateMarker(double branchlength, long whichmarker, const DoubleVec1d& state) const
{
    throw implementation_error("Cannot simulate data with the MixedKS model yet.");
} // SimulateMarker

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

string MixedKSModel::CellToData(Cell_ptr cell, long marker) const
{
    LongVec1d ones = cell->GetOnes(marker);
    assert(static_cast<size_t>(marker) < m_offsets.size());
    if (ones.size() == 1)
    {
        return ToString(ones[0] + m_offsets[marker]);
    }
    throw implementation_error
        ("Cannot convert stepwise data from the internal format for data not simply a single number.");
}

//____________________________________________________________________________________