C++ SiteContainer類代碼示例

void RecursiveLikelihoodTree::initLikelihoods(const SiteContainer& sites, const SubstitutionProcess& process)
throw (Exception)
{
  if (sites.getNumberOfSequences() == 1)
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. Only 1 sequence in data set.");
  if (sites.getNumberOfSequences() == 0)
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. No sequence in data set.");
  if (!process.isCompatibleWith(sites))
    throw Exception("RecursiveLikelihoodTree::initLikelihoods. Data and model are not compatible.");
  alphabet_ = sites.getAlphabet();
  nbStates_ = process.getNumberOfStates();
  nbSites_  = sites.getNumberOfSites();
  unique_ptr<SitePatterns> patterns;

  if (usePatterns_)
  {
    patterns.reset(initLikelihoodsWithPatterns_(process.getTree().getRootNode(), sites, process));
    shrunkData_.reset(patterns->getSites());
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_ = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();

    setPatterns(patternLinks_);
  }
  else
  {
    patterns.reset(new SitePatterns(&sites));
    shrunkData_.reset(patterns->getSites());
    rootWeights_      = patterns->getWeights();
    rootPatternLinks_ = patterns->getIndices();
    nbDistinctSites_  = shrunkData_->getNumberOfSites();
    initLikelihoodsWithoutPatterns_(process.getTree().getRootNode(), *shrunkData_, process);
  }
}

void SiteContainerTools::removeGapSites(SiteContainer& sites, double maxFreqGaps)
{
  for (size_t i = sites.getNumberOfSites(); i > 0; i--) {
    map<int, double> freq;
    SiteTools::getFrequencies(sites.getSite(i - 1), freq);
    if (freq[-1] > maxFreqGaps)
      sites.deleteSite(i - 1);
  }
}

void SiteContainerTools::changeUnresolvedCharactersToGaps(SiteContainer& sites)
{
  // NB: use iterators for a better algorithm?
  int gapCode = sites.getAlphabet()->getGapCharacterCode();
  for (unsigned int i = 0; i < sites.getNumberOfSites(); i++)
  {
    for (unsigned int j = 0; j < sites.getNumberOfSequences(); j++)
    {
      int* element = &sites(j, i);
      if (sites.getAlphabet()->isUnresolved(*element))
        *element = gapCode;
    }
  }
}

void SiteContainerTools::changeGapsToUnknownCharacters(SiteContainer& sites)
{
  // NB: use iterators for a better algorithm?
  int unknownCode = sites.getAlphabet()->getUnknownCharacterCode();
  for (unsigned int i = 0; i < sites.getNumberOfSites(); i++)
  {
    for (unsigned int j = 0; j < sites.getNumberOfSequences(); j++)
    {
      int* element = &sites(j, i);
      if (sites.getAlphabet()->isGap(*element))
        *element = unknownCode;
    }
  }
}

void RewardMappingTools::writeToStream(
  const ProbabilisticRewardMapping& rewards,
  const SiteContainer& sites,
  ostream& out)
throw (IOException)
{
  if (!out)
    throw IOException("RewardMappingTools::writeToFile. Can't write to stream.");
  out << "Branches";
  out << "\tMean";
  for (size_t i = 0; i < rewards.getNumberOfSites(); i++)
  {
    out << "\tSite" << sites.getSite(i).getPosition();
  }
  out << endl;

  for (size_t j = 0; j < rewards.getNumberOfBranches(); j++)
  {
    out << rewards.getNode(j)->getId() << "\t" << rewards.getNode(j)->getDistanceToFather();
    for (size_t i = 0; i < rewards.getNumberOfSites(); i++)
    {
      out << "\t" << rewards(j, i);
    }
    out << endl;
  }
}

/*
 * Inheriting from SubstitutionProcess
 */
bool SubstitutionProcessCollectionMember::isCompatibleWith(const SiteContainer& data) const
{
  if (modelToNodes_.size() > 0)
    return data.getAlphabet()->getAlphabetType() == pSubProColl_->getModel(modelToNodes_.begin()->first).getAlphabet()->getAlphabetType();
  else
    return true;
}

SiteContainer* SequenceSimulationTools::simulateSites(const SiteSimulator& simulator, const vector<double>& rates)
{
  size_t numberOfSites = rates.size();
  vector<const Site*> vs(numberOfSites);
  for (size_t i = 0; i < numberOfSites; i++)
  {
    Site* s = simulator.simulateSite(rates[i]);
    s->setPosition(static_cast<int>(i));
    vs[i] = s;
  }
  SiteContainer* sites = new VectorSiteContainer(vs, simulator.getAlphabet());
  sites->setSequencesNames(simulator.getSequencesNames(), false);
  // Freeing memory:
  for (size_t i = 0; i < numberOfSites; i++)
  {
    delete vs[i];
  }

  return sites;
}

void DCSE::appendAlignmentFromStream(istream& input, SiteContainer& sc) const throw (Exception)
{
  // Checking the existence of specified file
  if (!input) { throw IOException ("DCSE::read : fail to open file"); }

  // Initialization
  const Alphabet * alpha = sc.getAlphabet();
  string line, name, sequence = "";

  line = FileTools::getNextLine(input); // Copy current line in temporary string
  //StringTokenizer st(line);
  //st.nextToken();
  //First line ignored for now!
  //int n1 = TextTools::toInt(st.nextToken());
  //int n2 = TextTools::toInt(st.nextToken());
  //int nbSites = n2 - n1
  //cout << nbSpecies << " species and " << nbSites << " sites." << endl;

  // Main loop : for all file lines
  while (!input.eof())
  {
    line = FileTools::getNextLine(input); // Copy current line in temporary string
    if(line == "") break;
    string::size_type endOfSeq = line.find("     ");
    if(endOfSeq == line.npos) break;
    sequence = string(line.begin(), line.begin() + static_cast<ptrdiff_t>(endOfSeq));
    sequence = TextTools::removeWhiteSpaces(sequence);
    sequence = TextTools::removeChar(sequence, '{');
    sequence = TextTools::removeChar(sequence, '}');
    sequence = TextTools::removeChar(sequence, '[');
    sequence = TextTools::removeChar(sequence, ']');
    sequence = TextTools::removeChar(sequence, '(');
    sequence = TextTools::removeChar(sequence, ')');
    sequence = TextTools::removeChar(sequence, '^');
    name     = string(line.begin() + static_cast<ptrdiff_t>(endOfSeq + 1), line.end()),
    name     = TextTools::removeFirstWhiteSpaces(name);
    if(name.find("Helix numbering") == name.npos
    && name.find("mask") == name.npos)
      sc.addSequence(BasicSequence(name, sequence, alpha), true);
  }
}

void SiteContainerTools::removeGapOrUnresolvedOnlySites(SiteContainer& sites)
{
  size_t n = sites.getNumberOfSites();
  size_t i = n;
  while (i > 1)
  {
    ApplicationTools::displayGauge(n - i + 1, n);
    const Site* site = &sites.getSite(i - 1);
    if (SiteTools::isGapOnly(*site))
    {
      size_t end = i;
      while (SiteTools::isGapOrUnresolvedOnly(*site) && i > 1)
      {
        --i;
        site = &sites.getSite(i - 1);
      }
      sites.deleteSites(i, end - i);
    }
    else
    {
      --i;
    }
  }
  ApplicationTools::displayGauge(n, n);
  const Site* site = &sites.getSite(0);
  if (SiteTools::isGapOrUnresolvedOnly(*site))
    sites.deleteSite(0);
}

void SiteContainerTools::merge(SiteContainer& seqCont1, const SiteContainer& seqCont2, bool leavePositionAsIs)
throw (AlphabetMismatchException, Exception)
{
  if (seqCont1.getAlphabet()->getAlphabetType() != seqCont2.getAlphabet()->getAlphabetType())
    throw AlphabetMismatchException("SiteContainerTools::merge.", seqCont1.getAlphabet(), seqCont2.getAlphabet());


  vector<string> seqNames1 = seqCont1.getSequencesNames();
  vector<string> seqNames2 = seqCont2.getSequencesNames();
  const SiteContainer* seqCont2bis = 0;
  bool del = false;
  if (seqNames1 == seqNames2)
  {
    seqCont2bis = &seqCont2;
  }
  else
  {
    // We shall reorder sequences first:
    SiteContainer* seqCont2ter = new VectorSiteContainer(seqCont2.getAlphabet());
    SequenceContainerTools::getSelectedSequences(seqCont2, seqNames1, *seqCont2ter);
    seqCont2bis = seqCont2ter;
    del = true;
  }

  if (leavePositionAsIs)
  {
    for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); i++)
    {
      seqCont1.addSite(seqCont2bis->getSite(i), false);
    }
  }
  else
  {
    int offset = static_cast<int>(seqCont1.getNumberOfSites());
    for (size_t i = 0; i < seqCont2bis->getNumberOfSites(); i++)
    {
      seqCont1.addSite(seqCont2bis->getSite(i), offset + seqCont2bis->getSite(i).getPosition(), false);
    }
  }

  if (del)
    delete seqCont2bis;
}

AbstractTreeParsimonyScore::AbstractTreeParsimonyScore(
  const Tree& tree,
  const SiteContainer& data,
  const StateMap* statesMap,
  bool verbose)
throw (Exception) :
  tree_(new TreeTemplate<Node>(tree)),
  data_(0),
  alphabet_(data.getAlphabet()),
  statesMap_(statesMap),
  nbStates_(statesMap->getNumberOfModelStates())
{
  init_(data, verbose);
}

SiteContainer* SequenceSimulationTools::simulateSites(const SiteSimulator& simulator, const vector<double>& rates, const vector<size_t>& states)
throw (Exception)
{
  size_t numberOfSites = rates.size();
  if (states.size() != numberOfSites)
    throw Exception("SequenceSimulationTools::simulateSites., 'rates' and 'states' must have the same length.");
  vector<const Site*> vs(numberOfSites);
  for (size_t i = 0; i < numberOfSites; i++)
  {
    Site* s = simulator.simulateSite(states[i], rates[i]);
    s->setPosition(static_cast<int>(i));
    vs[i] = s;
  }
  SiteContainer* sites = new VectorSiteContainer(vs, simulator.getAlphabet());
  sites->setSequencesNames(simulator.getSequencesNames(), false);
  // Freeing memory:
  for (size_t i = 0; i < numberOfSites; i++)
  {
    delete vs[i];
  }

  return sites;
}

std::map<size_t, size_t> SiteContainerTools::translateSequence(const SiteContainer& sequences, size_t i1, size_t i2)
{
  const Sequence* seq1 = &sequences.getSequence(i1);
  const Sequence* seq2 = &sequences.getSequence(i2);
  map<size_t, size_t> tln;
  size_t count1 = 0; // Sequence 1 counter
  size_t count2 = 0; // Sequence 2 counter
  int state1;
  int state2;
  for (size_t i = 0; i <  sequences.getNumberOfSites(); i++)
  {
    state1 = (*seq1)[i];
    if (state1 != -1)
      count1++;
    state2 = (*seq2)[i];
    if (state2 != -1)
      count2++;
    if (state1 != -1)
    {
      tln[count1] = (state2 == -1 ? 0 : count2);
    }
  }
  return tln;
}

AbstractTreeParsimonyScore::AbstractTreeParsimonyScore(
  const Tree& tree,
  const SiteContainer& data,
  bool verbose,
  bool includeGaps)
throw (Exception) :
  tree_(new TreeTemplate<Node>(tree)),
  data_(0),
  alphabet_(data.getAlphabet()),
  statesMap_(0),
  nbStates_(0)
{
  statesMap_ = new CanonicalStateMap(alphabet_, includeGaps);
  nbStates_  = statesMap_->getNumberOfModelStates();
  init_(data, verbose);
}

VectorSiteContainer::VectorSiteContainer(const SiteContainer& sc) :
  AbstractSequenceContainer(sc),
  sites_(0),
  names_(sc.getSequencesNames()),
  comments_(sc.getNumberOfSequences()),
  sequences_(sc.getNumberOfSequences())
{
  // Now try to add each site:
  for (size_t i = 0; i < sc.getNumberOfSites(); i++)
  {
    addSite(sc.getSite(i), false); // We assume that positions are correct.
  }
  // Seq comments:
  for (size_t i = 0; i < sc.getNumberOfSequences(); i++)
  {
    comments_[i] = new Comments(sc.getComments(i));
  }
}