src/datalike/locus.cpp

// $Id: locus.cpp,v 1.81 2012/01/05 01:28:42 ewalkup Exp $

/*
  Copyright 2003  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 <functional>                   // for Locus::CountNNucleotides()
#include <fstream>

#include "local_build.h"

#include "constants.h"
#include "dlcalc.h"
#include "dlmodel.h"
#include "force.h"                      // for TipData::GetBranchPartitions()
#include "individual.h"                 // for use of Individuals in setting up Locus objects
#include "locus.h"
#include "mathx.h"                      // for IsEven
#include "rangex.h"                     // LS DEBUG SIM
#include "registry.h"
#include "runreport.h"
#include "stringx.h"
#include "xml_strings.h"

using namespace std;

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

Locus::Locus(long int ind, bool movable, string name)
    : m_index(ind),
      m_name(name),
      m_nmarkers(FLAGLONG),
      m_nsites(FLAGLONG),
      m_regionalmapposition(0),
      m_globalmapposition(0),
      m_offset(0),
      m_movable(movable),
      m_type(mloc_data),
      m_defaultlocations(true),
      m_positions(),
      m_pDatatype(),
      m_pDatamodel(),
      m_tipdata(),
      m_protoCell(),
      m_pDLCalculator(),
      m_tipcells(),
      m_allowedrange(),
      m_variablerange(),
      m_unknownrange(),
      m_map(),
      m_simulate(false),
      m_truesite(FLAGLONG),
      m_variability(),
      m_phenotypes(name)
{
    if (movable)
    {
        //this is a different sort of locus:
        SetNmarkers(1); //This requires alleles of 1 marker
        SetNsites(1);
        SetPositions();
    }
} // Locus constructor

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

string Locus::GetName() const
{
    if (!m_name.empty()) return m_name;
    string tempname("#");
    tempname += ToString(GetIndex()+1);
    return tempname;
} // GetName

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

long int Locus::GetNsites() const
{
    // if it has never been set, we assume it's m_nmarkers
    if (m_nsites == FLAGLONG) return m_nmarkers;
    return m_nsites;

} // GetNsites

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

long int Locus::GetOffset() const
{
    return m_offset;

} // GetOffset

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

LongVec1d Locus::GetUserMarkerLocations() const
{
    LongVec1d userpos(m_positions);
    std::transform(userpos.begin(),userpos.end(),userpos.begin(),
                   std::bind2nd(std::plus<long int>(),m_offset-m_regionalmapposition));
    return userpos;
} // GetUserMarkerLocations

//------------------------------------------------------------------------------------
// Once everything is ready, make this Locus into a fully functional one containing likelihood cells for its tips.

void Locus::Setup(const IndVec& individuals)
{
    m_tipcells.clear();

    m_protoCell = m_pDatatype->CreateDLCell(*this);
    unsigned long int tip;
    for (tip = 0; tip < m_tipdata.size(); ++tip)
    {
        if (m_tipdata[tip].m_nodata)
        {
            //The data is stored in the haplotypes, not this tip.
            vector<LocusCell> cellsbymarkers;
            //Find which individual codes for this tip.
            for (unsigned long int ind = 0; ind < individuals.size(); ind++)
            {
                if (m_tipdata[tip].individual == individuals[ind].GetId())
                {
                    for (long int marker = 0; marker < GetNmarkers(); marker++)
                    {
                        vector<LocusCell> cells = individuals[ind].GetLocusCellsFor(GetName(), marker);
                        assert(static_cast<long int>(cells.size()) > m_tipdata[tip].m_hap);
                        cellsbymarkers.push_back(cells[m_tipdata[tip].m_hap]);
                    }
                    continue;
                }
            }
            LocusCell onecell(cellsbymarkers);
            m_tipcells.push_back(onecell);
        }
        else
        {
            m_tipcells.push_back(m_pDatatype->CreateInitializedDLCell(*this, m_tipdata[tip].data));
        }
    }
    m_pDLCalculator = DLCalc_ptr(m_pDatatype->CreateDLCalculator(*this));

} // Setup

//------------------------------------------------------------------------------------
// took away clone of src since it should be uniquely generated
// for each locus by Registry::InstallDataModels

void Locus::SetDataModelOnce(DataModel_ptr src)
{
    if (src.get() != NULL) m_pDatamodel = src;
    else m_pDatamodel.reset();

} // SetDataModel

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

void Locus::SetAnalysisType(mloc_type type)
{
    m_type = type;
    switch(type)
    {
        case mloc_mapjump:
        case mloc_mapfloat:
            m_movable = true;
            break;
        case mloc_data:
            m_movable = false;
            break;
        case mloc_partition:
            assert(false);
            throw implementation_error("You shouldn't be able to set the analysis type to 'partition' yet.");
            break;
            //LS DEBUG MAPPING:  We need to throw it away here or something.
    }
}

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

void Locus::SetAllowedRange(rangeset rs, long int regoffset)
{
    //rs comes in here on the global scale.
    rangeset offsetrs;
    for (rangesetiter rpair = rs.begin(); rpair != rs.end(); rpair++)
    {
        long int low = rpair->first - regoffset;
        long int high = rpair->second - regoffset;
        offsetrs.insert(make_pair(low, high));
    }

    m_allowedrange = offsetrs;
}

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

void Locus::SetEmptyTipData(vector<TipData> td)
{
    for (unsigned long int tip = 0; tip < td.size(); tip++)
    {
        td[tip].data.clear();
        td[tip].m_nodata = true;
    }
    m_tipdata = td;
}

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

void Locus::SetVariableRange(rangeset rs)
{
    m_variablerange = rs;
}

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

void Locus::SetDataType(DataType_ptr src)
{
    m_pDatatype = src; // shallow copy!
} // SetDataType

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

void Locus::SetNmarkers(long int n)
{
    if (m_nmarkers == FLAGLONG)
    {
        m_nmarkers = n;
    }
    else
    {
        if (m_nmarkers != n)
        {
            data_error e("Inconsistent number of markers");
            throw e;
        }
    }
} // SetNmarkers

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

void Locus::SetGlobalMapposition(long int site)
{
#if 0  // Potentially DEAD CODE (bobgian, Feb 2010)
    if (site == 0)
    {
        throw data_error("Assuming the biologist's convention of the nonexistence of site zero,"
                         " we assume that the position left of site 1 is site -1."
                         "  As such, you may not set the map position of any segment to '0'.");
    }
#endif
    m_globalmapposition = site;
} // Setglobalmapposition

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

void Locus::SetRegionalMapposition(long int site)
{
    transform(m_positions.begin(),m_positions.end(),m_positions.begin(),
              std::bind2nd(std::minus<long int>(),m_regionalmapposition));

    if(IsMoving())
    {
        long int movement = site - m_regionalmapposition;
        SetGlobalMapposition(m_globalmapposition + movement);
    }

    m_regionalmapposition = site;
    std::transform(m_positions.begin(),m_positions.end(),m_positions.begin(),
                   std::bind2nd(std::plus<long int>(),m_regionalmapposition));

} // Setregionalmapposition

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

void Locus::SetOffset(long int val)
{
    m_offset = val;
} // SetOffset

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

void Locus::SetPositions()
{
    m_defaultlocations = true;
    m_positions.clear();
    m_positions.reserve(m_nmarkers);  // for speed
    long int i;
    for (i = 0; i < m_nmarkers; ++i)
    {
        m_positions.push_back(i + m_regionalmapposition);
    }

} // SetPositions

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

void Locus::SetPositions(const LongVec1d& pos)
{
    m_defaultlocations = false;
    m_positions = pos;
    std::transform(m_positions.begin(),m_positions.end(),m_positions.begin(),
                   std::bind2nd(std::plus<long int>(),m_regionalmapposition));

} // SetPositions

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

LongVec1d Locus::CalcNVariableMarkers() const
{
    LongVec1d nvarmarkers;

    const DataPack& dpack = registry.GetDataPack();
    long int xpart, nxparts = dpack.GetNCrossPartitions();

    for (xpart = 0; xpart < nxparts; ++xpart)
    {
        nvarmarkers.push_back(CalcNVariableMarkers(dpack.GetTipId(xpart)));
    }

    return nvarmarkers;
} // CalcNVariableMarkers()

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

long int Locus::CalcNVariableMarkers(tipidtype xpart) const
{
    long int nvarmarkers = 0;

    const StringVec2d data = GetCrossPartitionGeneticData(xpart);
    if (!data.empty()) nvarmarkers = GetDataTypePtr()->CalcNVarMarkers(data);

    return nvarmarkers;

} // CalcNVariableMarkers(tipidtype xpart)

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

vector<TipData> Locus::GetPopulationTipData(const string& popname) const
{
    vector<TipData> popdata;
    vector<TipData>::const_iterator tip;
    for(tip = m_tipdata.begin(); tip != m_tipdata.end(); ++tip)
        if (tip->IsInPopulation(popname)) popdata.push_back(*tip);

    return popdata;

} // GetPopulationTipData

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

StringVec2d Locus::GetCrossPartitionGeneticData(tipidtype xpart) const
{
    StringVec2d data;
    vector<TipData>::const_iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
        if (tip->IsInCrossPartition(xpart)) data.push_back(tip->data);

    return data;

} // GetCrossPartitionGeneticData

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

StringVec3d Locus::GetPartitionGeneticData(force_type partname) const
{
    assert(partname == force_MIG || partname == force_DISEASE || partname == force_DIVMIG);

    StringVec3d data(registry.GetDataPack().GetNPartitionsByForceType(partname));
    vector<TipData>::const_iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
        data[tip->GetPartition(partname)].push_back(tip->data);

    return data;

} // GetPartitionGeneticData

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

StringVec1d Locus::GetMarkerDataWithLabels(const IndVec& individuals) const
{
    long int width = std::max(static_cast<long int>(GetName().size()), GetNmarkers()+5);
    StringVec1d labeleddata(1,MakeCentered(GetName(),width));
    vector<TipData>::const_iterator tip = GetTipData().begin();
    if (tip->m_nodata)
    {
        //We need to iterate over the individuals instead
        //LS DEBUG:  this is a pretty fragile check for this situation.
        for (long int marker = 0; marker < m_nmarkers; marker++)
        {
            for (unsigned long int ind = 0; ind < individuals.size(); ind++)
            {
                string label = MakeJustified(individuals[ind].GetName(), -9);
                string data = individuals[ind].GetMarkerDataFor(GetName(), marker);
                labeleddata.push_back(label + " " + data);
            }
            if (m_nmarkers > 1)
            {
                labeleddata.push_back("");
            }
        }
    }
    else
    {
        for( ; tip != GetTipData().end(); ++tip)
        {
            string label = MakeJustified(tip->label,-9);
            string data = tip->GetFormattedData(m_pDatatype->GetDelimiter());
            labeleddata.push_back(label+ " " + data);
        }
    }
    return labeleddata;

} // GetMarkerDataWithLabels

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

DoubleVec1d Locus::CountNNucleotides() const
{
    DoubleVec1d count(BASES,0L);

    if (!m_pDatatype->IsNucleotideData()) return count;

    vector<TipData>::const_iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
    {
        StringVec1d data = tip->data;

        StringVec1d::const_iterator sit;
        for (sit = data.begin(); sit != data.end(); ++sit)
        {
            // we can't use locus' inherent datamodel here because this code may
            // be called in the menu, where the locus may not have a datamodel yet!
            DoubleVec1d site = NucModel::StaticDataToLikes(*sit,GetPerBaseErrorRate());
            double zero(0.0);
            double total(accumulate(site.begin(),site.end(),zero));
            transform(site.begin(),site.end(),site.begin(),
                      bind2nd(divides<double>(),total));

            assert(site.size() == count.size());

            transform(count.begin(),count.end(),site.begin(),
                      count.begin(),plus<double>());
        }
    }
    return count;
} // CountNNucleotides

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

bool Locus::MultiSampleIndividuals() const
{
    set<long int> individuals; //a unique list
    for (unsigned long int ind = 0; ind < m_tipdata.size(); ind++)
    {
        long int individual = m_tipdata[ind].individual;
        if (individuals.find(individual) != individuals.end()) return true;
        individuals.insert(individual);
    }
    return false;
}

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

long int Locus::GetNTips(tipidtype xpart) const
{
    long int count = 0;
    vector<TipData>::const_iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
        if (tip->IsInCrossPartition(xpart)) ++count;

    return count;
} // GetNTips(tipidtype xpart)

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

double Locus::CalculateDataLikelihood(Tree& tree, bool moving) const
{
    return m_pDLCalculator->Calculate(tree, *this, moving);
} // CalculateDataLikelihood

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

void Locus::AddUniqueNamesTo(set<string>& popnames) const
{
    vector<TipData>::const_iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
        popnames.insert(tip->m_popname);  // only new names added because popnames is a std::set

} // AddUniqueNamesTo

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

void Locus::SetNewMappositionIfMoving()
{
    assert (m_simulate);
    if (IsMoving())
    {
        //Pick a site to actually live
        long int nsites = CountSites(m_allowedrange);
        do {
            m_truesite = registry.GetRandom().Long(nsites);
        } while (m_allowedrange != AddPairToRange(make_pair(m_truesite, m_truesite+1), m_allowedrange));
        SetRegionalMapposition(m_truesite);
    }
}

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

void Locus::SimulateData(Tree& tree, long int nsites)
{
    assert(m_simulate);
    ClearVariability();
    long int ntries = 0;
    while (IsNighInvariant() && ++ntries<5000)
    {
        m_pDLCalculator->SimulateData(tree, *this);
    }
    if (ntries >= 5000)
    {
        registry.GetRunReport().ReportNormal("Gave up trying to simulate non-invariant data for segment "
                                             + GetName() + ".");
    }
} // SimulateData

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

void Locus::CopyDataFrom(Locus& original, Tree& tree)
{
    m_pDLCalculator->CopyDataFrom(*this, original, tree);
}

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

void Locus::MakePhenotypesFor(IndVec& individuals)
{
    for (size_t ind = 0; ind<individuals.size(); ind++)
    {
        for (long int marker = 0; marker<m_nmarkers; marker++)
        {
            StringVec1d alleles = individuals[ind].GetAllelesFromDLs(m_index, marker, IsMoving(), m_pDatamodel);
            //LS DEBUG
            //cout << "Original haplotypes:  " << ToString(alleles) << endl;
            if (m_phenotypes.AnyDefinedPhenotypes())
            {
                individuals[ind].SetHaplotypes(GetName(), marker, m_phenotypes.ChooseHaplotypes(alleles));
            }
            else
            {
                Haplotypes haps(individuals[ind].GetHaplotypesFor(GetName(),marker), true);
                haps.AddHaplotype(alleles, 1.0);
                individuals[ind].SetHaplotypes(GetName(), marker, haps);
                //LS DEBUG
                //cout << "Final haplotypes:  " << ToString(haps.GetAlleles()) << endl;
                // individuals[ind].PrintHaplotypesFor(GetName(), marker);
            }
        }
    }
}

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

void Locus::SaveState(DoubleVec1d& state, long int marker, string label)
{
    //First, copy over our own dlcell:
    for (size_t tip = 0; tip < m_tipdata.size(); tip++)
    {
        if (m_tipdata[tip].label == label)
        {
            m_tipcells[tip].SetAllCategoriesTo(state, marker);
        }
    }

    //Now add it to the 'variability' list.
    map<long int, DoubleVec1d>::iterator freqs = m_variability.find(marker);
    if (freqs == m_variability.end())
    {
        m_variability.insert(make_pair(marker, state));
        return;
    }
    //Not new, so add it to the old
    std::transform(  state.begin(),
                     state.end(),
                     freqs->second.begin(),
                     freqs->second.begin(),
                     std::plus<double>());
}

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

void Locus::RandomizeHalf(Tree& tree, bool swath)
{
    m_unknownrange.clear();
    //The idea here is that we take a random quarter-sized swath out of the first
    // half, and and another quarter-sized swath out of the second half.

    if (swath)
    {
        //Pick a swath out of the first half, and another swath out of the
        // second half.
        long int quarter = static_cast<long int>(m_nsites/4);
        long int left = registry.GetRandom().Long(quarter);
        m_unknownrange.insert(make_pair(left, left+quarter+1));
        left = registry.GetRandom().Long(quarter);
        m_unknownrange.insert(make_pair(left+(2*quarter), left + (3*quarter) + 1));
    }
    else
    {
        //The 'every other site' version:
        for (long int left = 0; left < m_nsites; left = left+2)
        {
            m_unknownrange.insert(make_pair(left, left+1));
        }
    }
    m_pDLCalculator->Randomize(*this, m_unknownrange, tree);

    //Now we need to clear the variable sites out of the unknown range
    m_variablerange = RemoveRangeFromRange(m_unknownrange, m_variablerange);
}

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

bool Locus::SiteInLocus(long int site) const
{
    // is the given site in this locus?
    return (m_regionalmapposition <= site &&
            site < m_regionalmapposition + m_nsites);
} // SiteInLocus

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

long int Locus::SiteToMarker(long int site) const
{
    assert(SiteInLocus(site));
    for (long int i = 0; i < m_nmarkers; i++)
    {
        if (m_positions[i] == site)
        {
            return i;
        }
    }
    return FLAGLONG;
}

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

pair<long int, long int> Locus::GetSiteSpan() const
{
    return make_pair(m_regionalmapposition, m_regionalmapposition + m_nsites);
} // GetSiteSpan

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

pair<long int, long int> Locus::GetGlobalScaleSiteSpan() const
{
    return make_pair(m_globalmapposition, m_globalmapposition + m_nsites);
} // GetGlobalScaleSiteSpan

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

bool Locus::IsMovable() const
{
    //We might want to change this if split into phase 1/phase 2
    return m_movable;
}

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

bool Locus::IsMoving() const
{
    switch (GetAnalysisType())
    {
        case mloc_mapjump:
        case mloc_mapfloat:
            return true;
        case mloc_data:
            return false;
        case mloc_partition:
            assert(false);
            throw implementation_error("Shouldn't be checking this for a partition segment.");
    }
    throw implementation_error("Uncaught analysis type for segment.");
}

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

double Locus::GetPerBaseErrorRate() const
// EWFIX.CODESMELL -- horrible, horrible code smell
{
    const DataModel_ptr p = GetDataModel();
    if(p == NULL) return defaults::per_base_error_rate;

    DataModel * q = p->Clone();
    NucModel * nm = dynamic_cast<NucModel*>(q);
    if(nm != NULL) return nm->GetPerBaseErrorRate();

    return defaults::per_base_error_rate;
}

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

void Locus::RemovePartitionFromTipDatas(force_type forcename)
{
    vector<TipData>::iterator tip = GetTipData().begin();
    for( ; tip != GetTipData().end(); ++tip)
        tip->RemovePartition(forcename);

} // RemovePartitionFromTipDatas

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

bool Locus::IsValidLocus(string& errorString) const
{
    unsigned long int nmark = GetNmarkers();  // to avoid comparison warning
    if (nmark == 0)
    {
        errorString = "No data in segment " + GetName();
        return false;
    }

    if (nmark != GetMarkerLocations().size())
    {
        errorString = "The number of markers doesn't match their positions in segment " + GetName();
        return false;
    }

    if (m_nsites < static_cast<long int>(nmark))
    {
        errorString = "Number of markers exceeds sequence length in segment " + GetName();
        return false;
    }

    vector<long int> sortedpositions = GetMarkerLocations();
    std::sort(sortedpositions.begin(), sortedpositions.end());
    if (sortedpositions != GetMarkerLocations())
    {
        errorString = "Positions out of order in segment " + GetName();
        return false;
    }

#if 0 // Dead code
    // There should no longer be a need to validate the data
    // model since they are constructed to be valid and the
    // participating members don't change.
    if (m_pDatamodel.get() != NULL && !m_pDatamodel->IsValidDataModel())
    {
        errorString = "Invalid datamodel in segment " + m_name;
        return false;
    }
#endif // Dead code

    // We needn't validate datatype as it has no state.
    return true;

} // IsValidLocus

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

set<pair<double, long int> > Locus::MakeOrderedSites(long int regoffset) const
{
    assert(m_map.size() > 0);

    set<pair<double, long int> > orderedsites;
    for (unsigned long int site = 0; site < m_map.size(); site++)
    {
        long int place = site + regoffset;
        orderedsites.insert(make_pair(m_map[site], place));
    }
    return orderedsites;
}

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

rangeset Locus::GetBestSites(set<pair<double, long int> > orderedsites) const
{
    rangeset bestsites;

    for (set<pair<double, long int> >::reverse_iterator sitepair=orderedsites.rbegin();
         sitepair != orderedsites.rend(); sitepair++)
    {
        rangepair thissite = make_pair(sitepair->second, sitepair->second+1);
        if (sitepair->first == orderedsites.rbegin()->first)
        {
            bestsites = AddPairToRange(thissite, bestsites);
        }
    }

    return bestsites;
}

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

rangeset Locus::GetTopSites(set<pair<double, long int> > orderedsites, double percLimit) const
{
    rangeset topsites;

    double total = 0;
    for (set<pair<double, long int> >::reverse_iterator sitepair=orderedsites.rbegin();
         sitepair != orderedsites.rend(); sitepair++)
    {
        rangepair thissite = make_pair(sitepair->second, sitepair->second+1);
        if (total < percLimit)
        {
            topsites = AddPairToRange(thissite, topsites);
        }
        total += sitepair->first;
        //We add afterwards so that the site that pushes us over the edge still gets
        // included in the appropriate range.
    }

    return topsites;
}

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

StringVec1d Locus::ReportMappingInfo(long int regoffset, bool isshort) const
{
    set<pair<double, long int> > orderedsites = MakeOrderedSites(regoffset);


    rangeset bestsites = GetBestSites(orderedsites);
    rangeset topfivepercent = GetTopSites(orderedsites,0.05);
    rangeset topfiftypercent = GetTopSites(orderedsites,0.5);
    rangeset topninetyfivepercent = GetTopSites(orderedsites,0.95);

    // EWFIX -- 2010-09-02 -- don't remove until checked
#if 0
    double total = 0;
    for (set<pair<double, long int> >::reverse_iterator sitepair=orderedsites.rbegin();
         sitepair != orderedsites.rend(); sitepair++)
    {
        rangepair thissite = make_pair(sitepair->second, sitepair->second+1);
        if (sitepair->first == orderedsites.rbegin()->first)
        {
            bestsites = AddPairToRange(thissite, bestsites);
        }
        if (total < .05)
        {
            topfivepercent = AddPairToRange(thissite, topfivepercent);
        }
        if (total < .5)
        {
            topfiftypercent = AddPairToRange(thissite, topfiftypercent);
        }
        if (total < .95)
        {
            topninetyfivepercent = AddPairToRange(thissite, topninetyfivepercent);
        }
        total += sitepair->first;
        //We add afterwards so that the site that pushes us over the edge still gets
        // included in the appropriate range.
    }
#endif

    StringVec1d report;
    string msg = "Most likely site(s) for " + GetName() + ":  "
        + ToString(bestsites) + ".  Relative data likelihood = "
        + ToString(orderedsites.rbegin()->first);
    report.push_back(msg);
    if (isshort)
    {
        return report;
    }
    msg = "The top 5% of all sites in this region:  " + ToString(topfivepercent);
    report.push_back(msg);
    msg = "The top 50% of all sites in this region:  " + ToString(topfiftypercent);
    report.push_back(msg);
    msg = "The top 95% of all sites in this region:  " + ToString(topninetyfivepercent);
    report.push_back(msg);
    msg = "You have a total of " + ToString(CountSites(topninetyfivepercent)) +
        " sites in your 95% range.";
    report.push_back(msg);
    if (m_truesite != FLAGLONG)
    {
        rangepair truesite = make_pair(m_truesite, m_truesite+1);
        msg = "The true site (" + ToString(truesite) + ") ";
        if (topninetyfivepercent == AddPairToRange(truesite, topninetyfivepercent))
        {
            msg += "was";
        }
        else
        {
            msg += "was not";
        }
        msg += " included in the top 95%.";
        report.push_back(msg);

#if 0  // LS DEBUG SIM  Hard-coded information output
        rangeset unknownrange = registry.GetDataPack().GetRegion(0).GetLocus(0).GetUnknownRange();
        msg = "This site was in the ";
        if (unknownrange == AddPairToRange(truesite, unknownrange))
        {
            msg += "known";
        }
        else
        {
            msg += "unknown";
        }
        msg += " part of the segment.";
        report.push_back(msg);
#endif
    }
    if (m_variability.size() > 0 && IsMoving())
    {
        report.push_back("Data variability:");
        for (map<long int, DoubleVec1d>::const_iterator marker = m_variability.begin();
             marker != m_variability.end(); marker++)
        {
            for (unsigned long int allele = 0; allele < marker->second.size(); allele++)
            {
                msg = "Position " + ToString(marker->first + 1) + ", "
                    + "Allele " + ToString(allele+1) + ":  "
                    + ToString(marker->second[allele]);
                report.push_back(msg);
            }
        }
    }
    return report;
}

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

bool Locus::IsDuplicateTipName(const string& newname) const
{
    vector<TipData>::const_iterator tip;
    for(tip = m_tipdata.begin(); tip != m_tipdata.end(); ++tip)
        if (tip->label == newname) return true;

    return false;
}

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

StringVec1d Locus::MakeTraitXML(long int nspaces, long int regoffset) const
{
    StringVec1d xmllines;
    string line = MakeIndent(MakeTag(xmlstr::XML_TAG_TRAIT), nspaces);
    xmllines.push_back(line);
    nspaces += INDENT_DEPTH;

    line = MakeTag(xmlstr::XML_TAG_NAME) + " " + GetName() + " " +
        MakeCloseTag(xmlstr::XML_TAG_NAME);
    xmllines.push_back(MakeIndent(line, nspaces));

    line = MakeTag(xmlstr::XML_TAG_ANALYSIS) + " "
        + ToXMLString(GetAnalysisType())
        + " " + MakeCloseTag(xmlstr::XML_TAG_ANALYSIS);
    xmllines.push_back(MakeIndent(line, nspaces));

    line = MakeTag(xmlstr::XML_TAG_POSSIBLE_LOCATIONS);
    xmllines.push_back(MakeIndent(line, nspaces));

    nspaces += INDENT_DEPTH;
    for (rangeset::iterator range = m_allowedrange.begin();
         range != m_allowedrange.end(); range++)
    {
        long int start = (*range).first + regoffset;
        long int end   = (*range).second + regoffset - 1;

        line = MakeTag(xmlstr::XML_TAG_RANGE);
        xmllines.push_back(MakeIndent(line, nspaces));

        nspaces += INDENT_DEPTH;
        line = MakeTag(xmlstr::XML_TAG_START) + " " + ToString(start)
            + " " + MakeCloseTag(xmlstr::XML_TAG_START);
        xmllines.push_back(MakeIndent(line, nspaces));
        line = MakeTag(xmlstr::XML_TAG_END) + " " + ToString(end)
            + " " + MakeCloseTag(xmlstr::XML_TAG_END);
        xmllines.push_back(MakeIndent(line, nspaces));
        nspaces -= INDENT_DEPTH;

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

    StringVec1d dlmodelxml(GetDataModel()->ToXML(nspaces));
    xmllines.insert(xmllines.end(),dlmodelxml.begin(),dlmodelxml.end());

    StringVec1d phenotypexml(m_phenotypes.GetPhenotypesXML(nspaces));
    xmllines.insert(xmllines.end(),phenotypexml.begin(),phenotypexml.end());

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

    return xmllines;
}

//------------------------------------------------------------------------------------
//The general rule for IsNighInvariant is that if any marker is *not* 'nigh
// invariant' (false), the locus as a whole is also not.  If all markers are
// indeed 'nigh invariant' (true), so is the locus as a whole.

bool Locus::IsNighInvariant() const
{
    if (m_variability.size() == 0)
    {
        return true;
    }
    for (map<long int, DoubleVec1d>::const_iterator marker = m_variability.begin();
         marker != m_variability.end(); marker++)
    {
        if (!(IsNighInvariant(marker->first))) return false;
    }
    return true;
}

//------------------------------------------------------------------------------------
// The rule for IsNighInvariant (well, in its current 11-28-05 form) is that
//  if all of the data at this marker is one particular allele but two, the
//  data is 'nigh invariant' (true).  If there are at least three alleles which
//  are members of the non-majority allele, the data is not 'nigh invariant'
//  (false).

bool Locus::IsNighInvariant(long int marker) const
{
    unsigned long int ntips = m_tipcells.size();

    map<long int, DoubleVec1d>::const_iterator freqs = m_variability.find(marker);
    if (freqs == m_variability.end()) return false;
    double maxfreq = 0;
    for (unsigned long int allele = 0; allele < freqs->second.size(); allele++)
    {
        maxfreq = max(maxfreq, freqs->second[allele]);
    }
    long int mindiff = std::min(3L, static_cast<long int>(ntips) - 1);
    if (maxfreq <= ntips - mindiff)
    {
        return false;
    }
    return true;
}

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

bool Locus::IsCompletelyInvariant(long int marker) const
{
    unsigned long int ntips = m_tipcells.size();

    map<long int, DoubleVec1d>::const_iterator freqs = m_variability.find(marker);
    if (freqs == m_variability.end()) return false;
    for (unsigned long int allele = 0; allele < freqs->second.size(); allele++)
    {
        if (0 < freqs->second[allele] && freqs->second[allele] < ntips)
        {
            return false;
        }
    }
    return true;
}

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

long int Locus::ChooseVariableSiteFrom(rangeset rset)
{
    LongVec1d variables;
    for (rangeset::iterator rit = rset.begin(); rit != rset.end(); rit++)
    {
        for (long int site = rit->first; site<rit->second; site++)
        {
            if (SiteInLocus(site))
            {
                if (!IsNighInvariant(SiteToMarker(site)))
                {
                    variables.push_back(site);
                }
            }
        }
    }
    if (variables.size() == 0)
    {
        //No variable sites!
        return FLAGLONG;
    }
    return variables[registry.GetRandom().Long(variables.size())];
}

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

rangeset Locus::GetVariableRange()
{
    if (m_variablerange.size() == 0)
    {
        m_variablerange = CalculateVariableRange();
    }
    return m_variablerange;
}

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

rangeset Locus::CalculateVariableRange() const
{
    rangeset variables;
    for (long int site = 0; site < m_nsites; site++)
    {
        if (!IsNighInvariant(SiteToMarker(site)))
        {
            variables = AddPairToRange(make_pair(site, site+1), variables);
        }
    }
    return variables;

}

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

rangeset Locus::CalculateCompleteVariableRange() const
{
    rangeset variables;
    for (long int site = 0; site < m_nsites; site++)
    {
        if (!IsCompletelyInvariant(SiteToMarker(site)))
        {
            variables = AddPairToRange(make_pair(site, site+1), variables);
        }
    }
    return variables;
}

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

rangeset Locus::GetVariableAndUnknownRange() const
{
    return Union(m_variablerange, m_unknownrange);
}

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

long int Locus::GetVariabilityOfUnknownRange() const
{
    long int numvariable = 0;
    for (rangeset::iterator range = m_unknownrange.begin(); range != m_unknownrange.end(); range++)
    {
        for (long int site = range->first; site < range->second; site++)
        {
            if (!IsNighInvariant(SiteToMarker(site)))
            {
                numvariable++;
            }
        }
    }
    return numvariable;
}

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

long int Locus::GetVariabilityOfKnownRange() const
{
    return CountSites(m_variablerange);
}

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

rangeset Locus::GetKnownRange() const
{
    rangeset retset;
    retset.insert(make_pair(0, m_nsites));
    retset = RemoveRangeFromRange(m_unknownrange, retset);
    return retset;
}

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

StringVec1d Locus::CreateDataModelReport(string regionname) const
{
    StringVec1d out = m_pDatamodel->CreateDataModelReport();
    StringVec1d head;
    head.push_back("Parameters of a " + m_pDatamodel->GetDataModelName()
                   + " model for the " + GetName() + " segment of the "
                   + regionname + " region");
    head.insert(head.end(), out.begin(), out.end());
    return head;
}

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

void Locus::WriteMapping(string regname, long int regoffset) const
{

    ofstream mapfile;
    mapfile.precision(10);
    UserParameters& userparams = registry.GetUserParameters();
    string fname = userparams.GetMapFilePrefix() + "_" + GetName() + ".txt";
    mapfile.open(fname.c_str(), ios::out );
    userparams.AddMapFileName(fname);

    mapfile << "Mapping results for " + GetName() + " from the region \"" + regname + "\".\n";

    switch (GetAnalysisType())
    {
        case mloc_mapjump:
            mapfile << "The analysis for this trait was performed by allowing the location of "
                "the trait marker to move from place to place as trees were created.\n";
            break;
        case mloc_mapfloat:
            mapfile << "This analysis for this trait was performed by collecting trees, then calculating the "
                "data likelihood of the trait marker at all allowed sites on those trees, and then averaging.\n";
            break;
        case mloc_data:
        case mloc_partition:
            assert(false); //These loci should not be in the moving locus vector.
            return;
    }

    StringVec1d mapinfo = ReportMappingInfo(regoffset);
    for (size_t i = 0; i < mapinfo.size(); i++)
    {
        mapfile << mapinfo[i] + "\n";
    }
    mapfile << "\nSite\tData likelihood\tRaw likelihood\n";

    DoubleVec1d results = GetMappingInfo();
    DoubleVec1d rawresults = GetRawMappingInfo();
    StringVec1d resultsStrings(results.size(), "-");
    StringVec1d rawresultsStrings(rawresults.size(), "-");
    rangeset validrange = GetAllowedRange();
    for (rangeset::iterator range=validrange.begin();
         range != validrange.end(); range++)
    {
        for (long int site = range->first; site<range->second; site++)
        {
            resultsStrings[site] = Pretty(results[site], 10);
            rawresultsStrings[site] = Pretty(rawresults[site], 10);
        }
    }
    for (size_t site = 0; site < results.size(); site += 1)
    {
        long int place = ToNoZeroesIfNeeded(static_cast<long int>(site) + regoffset);
        mapfile << ToString(place) + "\t" + resultsStrings[site] + "\t" + rawresultsStrings[site] + "\n";
    }
}

//------------------------------------------------------------------------------------
// Debugging function.

void Locus::PrintVariability()
{
    for (map<long int, DoubleVec1d>::const_iterator marker = m_variability.begin();
         marker != m_variability.end();
         marker++)
    {
        string msg = "Position " + ToString(marker->first + 1) + ": ";
        for (unsigned long int allele = 0; allele < marker->second.size(); allele++)
        {
            msg += ToString(marker->second[allele]) + ", ";
        }
        cout << msg << endl;
    }
    cout << "Overall, the variable sites are:  " << ToString(GetVariableRange()) << endl;
}

//------------------------------------------------------------------------------------
// Debugging function.

void Locus::PrintOnesAndZeroesForVariableSites()
{
    rangeset variability = GetVariableRange();
    long int newl = 249;
    cout << "Variability: ";
    for (long int site = 0; site < m_nsites; site++)
    {
        if (variability == AddPairToRange(make_pair(site, site + 1), variability))
        {
            cout << "1 ";
        }
        else
        {
            cout << "0 ";
        }
        if (site == newl)
        {
            cout << endl << "Variability: ";
            newl += 250;
        }
    }
    cout << endl;
}

//____________________________________________________________________________________