PolyMoSim.cpp

/***************************************************************************************************
*  The PolyMoSim project is distributed under the following license:
*  
*  Copyright (c) 2006-2025, Christoph Mayer, Leibniz Institute for the Analysis of Biodiversity Change,
*  Bonn, Germany
*  All rights reserved.
*  
*  Redistribution and use in source and binary forms, with or without
*  modification, are permitted provided that the following conditions are met:
*  1. Redistributions of source code (complete or in parts) must retain
*     the above copyright notice, this list of conditions and the following disclaimer.
*  2. Redistributions in binary form must reproduce the above copyright
*     notice, this list of conditions and the following disclaimer in the
*     documentation and/or other materials provided with the distribution.
*  3. All advertising materials mentioning features or any use of this software
*     e.g. in publications must display the following acknowledgement:
*     This product includes software developed by Christoph Mayer, Forschungsmuseum
*     Alexander Koenig, Bonn, Germany.
*  4. Neither the name of the organization nor the
*     names of its contributors may be used to endorse or promote products
*     derived from this software without specific prior written permission.
*  
*  THIS SOFTWARE IS PROVIDED BY CHRISTOPH MAYER ''AS IS'' AND ANY
*  EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
*  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHTHOLDER OR ITS ORGANISATION BE LIABLE FOR ANY
*  DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
*  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
*  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
*  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
*  (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
*  SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*  
*  IMPORTANT (needs to be included, if code is redistributed):
*  Please not that this license is not compatible with the GNU Public License (GPL)
*  due to paragraph 3 in the copyright. It is not allowed under any
*  circumstances to use the code of this software in projects distributed under the GPL.
*  Furthermore, it is not allowed to redistribute the code in projects which are
*  distributed under a license which is incompatible with one of the 4 paragraphs above.
*  
*  This project makes use of code coming from other projects. What follows is a complete
*  list of files which make use of external code. Please refer to the copyright within
*  these files.
*  
*  Files in tclap foler:         Copyright (c) 2003 Michael E. Smoot
*                                See copyright in tclap/COPYRIGHT file for details.	
*  discrete_gamma.c:             Copyright 1993-2004 by Ziheng Yang.
*                                See copyright in this file for details.
*  CRandom.h:                    Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura
*                                See copyright in this file for details.
***************************************************************************************************/

#include "PolyMoSim.h"
#include "BasicTree.h"
#include "BasicNode.h"
#include "mymodel.h"
#include "model_admin.h"
#include "tree_admin.h"
// #include "sequence_gen.h"
#include "CRandom.h"
#include <sstream>
#include <ctime>
#include "symbol-cartesian-product.h"
#include "CSequences3.1.h"
#include <iomanip>

using namespace std;

// enum output_format{nexus, fasta, phylip, log_format};

inline void add_or_count_pattern(std::map<faststring, unsigned> &m, faststring &x)
{
  std::map<faststring,unsigned>::iterator it;
  
  it = m.find(x);
  if (it == m.end() )
  {
    m[x] = 1;
  }
  else
  {
    ++it->second;
  }
}


void welcome(FILE *of, const char *s="")
{
  myPrint(of, s, welcome_str);
}

void welcome(ofstream &of, const char *s="")
{
  of << s << welcome_str;
}

void report_seed(FILE *of, const char *s="")
{
  if (global_verbosity >= 1)
  {
    fprintf(of, "%sRandom number generator: %s\n",   s, randomNumberGeneratorNames[global_rg]);
    fprintf(of, "%sSeed of random generator: %u\n\n", s, global_seed_random_generator);
  }
}

void print_results_header(ostream& os, unsigned repetitions, unsigned ntaxa,
                          tree_admin &tchef, enum_outputformat format)
{
  // Write comments to output file:
  if (format == outputformat_nexus) {
    os << "#Nexus" << endl
    << "[============================================================"  << endl
    << "| Generated by PolyMoSim, version " << VERSION << endl 
    << "|" << endl
    << "| Model file:              " << global_model_file << endl
    << "| Tree file:               " << global_tree_file << endl
    << "| Replicates:              " << repetitions << endl
    << "| Taxa:                    " << ntaxa       << endl
    << "| Random number seed:      " << global_seed_random_generator << endl
    << "| Random number generator: " << randomNumberGeneratorNames[global_rg] << endl
    << "|" << endl;
    tchef.print(os, 1, "| ");
    os << "\\============================================================]"
    << endl << endl;
  }
  else if (format == outputformat_fasta)  {
  }
  else if (format == outputformat_phylip || format == outputformat_phylip_no_spaces) {
  }
  else if (format == outputformat_logfile)
  {
    os << endl
    << "Paramters used for simulation:" << endl 
    <<  endl
    << "Model file:              " << global_model_file << endl
    << "Tree file:               " << global_tree_file << endl
    << "Replicates:              " << repetitions << endl
    << "Taxa:                    " << ntaxa       << endl
    << "Random number seed:      " << global_seed_random_generator << endl
    << "Random number generator: " << randomNumberGeneratorNames[global_rg] << endl
    << endl;
    tchef.print(os, 1, "");
    os << endl << endl;
  }
}

void print_results_preanalysis(ostream& os, unsigned rep, unsigned ntaxa, unsigned nchar, const faststring& dt, enum_outputformat format)
{
  UNUSED(rep);
  
  char c;
  if (!global_preanalysis_file.empty()) {
    ifstream is(global_preanalysis_file.c_str());
    while ( is.get(c) )
    {
      os.put(c);
    }
    is.close();
  }
  
  if (format == outputformat_nexus) {
    os << "Begin data;" << endl
    << "   Dimensions ntax=" << ntaxa << " nchar=" << nchar << ";" << endl
    << "   Format missing=? gap=- datatype=" << dt << ";" << endl
    << "   Matrix" << endl;
  }
  else if (format == outputformat_fasta) {
  }
  else if (format == outputformat_phylip || format == outputformat_phylip_no_spaces ) {
    os << "\t" << ntaxa << "\t" << nchar << endl;
  }
}

void print_results_postanalysis(ostream& os, unsigned rep, unsigned ntaxa,
                                unsigned nchar, const faststring& dt, enum_outputformat format)
{
  UNUSED(rep);
  UNUSED(ntaxa);
  UNUSED(nchar);
  UNUSED(dt);
  
  
  char c;
  if (format == outputformat_nexus) {
    os << ";"          << endl;
    os << "end;"  << endl << endl;
  }
  else if (format == outputformat_fasta) {
  }
  else if (format == outputformat_phylip || format == outputformat_phylip_no_spaces) {
  }
  if ( !global_postanalysis_file.empty() ) {
    ifstream is(global_postanalysis_file.c_str() );
    while ( is.get(c) )
    {
      os.put(c);
    }
    is.close();
  }
}

void print_results_simdata(ostream& os, vector<BasicTree::map_of_OTUs> &OTU_data, unsigned rep, unsigned ntaxa,
                           unsigned nchar, enum_outputformat format, int datatype)
{
  //  UNUSED(rep);
  UNUSED(ntaxa);
  UNUSED(nchar);
  
  unsigned N_OTUs       = (unsigned)OTU_data[0].size(); // Number of OTUs.
  unsigned N_partitions = (unsigned)OTU_data.size();    // Number of partitions.
  
  if (1)
  {
    BasicTree::map_of_OTUs::iterator  it, it_end;
    it     = OTU_data[0].begin();
    it_end = OTU_data[0].end();
    
    vector<unsigned> len_vec(N_OTUs, 0);
    
    // Check integrety:
    unsigned i_OTU=0;
    while (it != it_end)                 // for all OTUs
    {
      for (unsigned i=0; i < N_partitions; ++i)     // for all partitions
        len_vec[i_OTU] += OTU_data[i][it->first]->size();
      ++it;
      ++i_OTU;
    }
    
    bool equal_len = true;
    for (unsigned j=1; j < len_vec.size(); ++j)
    {
      if (len_vec[j-1] != len_vec[j])
        equal_len = false;
    }
    if (!equal_len)
    {
      cerr << "WARNING: Sequnces have unequal total lengths:\n";
      for (unsigned j=0; j < len_vec.size(); ++j)
      {
        cerr << len_vec[j] << ',';
      }
      cerr << endl;
      cerr << "INFO: Length of simulated sequences: " << len_vec[0] << endl;
      
      unsigned i;
      it     = OTU_data[0].begin();
      it_end = OTU_data[0].end();
      
      faststring fileout = "Problematic-alignment_"+faststring(rep)+".fas";
      ofstream os_al(fileout.c_str());
      
      while (it != it_end) {
        os_al << ">" << *(it->first) << endl;
        for (i=0; i<OTU_data.size(); ++i)
          os_al << *OTU_data[i][it->first] << " ";
        os_al << endl;
        ++it;
      }
      os_al.close();
    }
  }
  
  unsigned i;
  BasicTree::map_of_OTUs::iterator  it, it_end;
  it     = OTU_data[0].begin();
  it_end = OTU_data[0].end();
  
  if (format == outputformat_nexus) {                        //
    while (it != it_end) {                      // for all taxa
      os << *(it->first) << " ";              //
      for (i=0; i < OTU_data.size(); ++i)       // for all partitions
        os << *OTU_data[i][it->first] << " "; //
      os << endl;
      ++it;
    }
  }
  else if (format == outputformat_fasta) {
    while (it != it_end) {
      os << ">" << *(it->first) << endl;
      for (i=0; i<OTU_data.size(); ++i)
        os << *OTU_data[i][it->first] << " ";
      os << endl;
      ++it;
    }
  } else if (format == outputformat_phylip) {
    while (it != it_end) {
      os << *(it->first) << " ";
      for (i=0; i<OTU_data.size(); ++i)
        os << *OTU_data[i][it->first] << " ";
      os << endl;
      ++it;
    }
  } else if (format == outputformat_phylip_no_spaces) {
    while (it != it_end) {
      os << *(it->first) << " ";
      for (i=0; i<OTU_data.size(); ++i)
        os << *OTU_data[i][it->first];
      os << endl;
      ++it;
    }
  } else if (format == outputformat_pattern_absolute      || format == outputformat_pattern_relative      ||
             format == outputformat_pattern_absolute_fill || format == outputformat_pattern_relative_fill )
  {
    //pointer to fastring sequences with alignment:
    unsigned total_sites = 0;
    
    faststring pattern( (size_t)N_OTUs, (char)' ');
    
    char **aln = new char*[N_OTUs]; // array of pointers to OTU sequences
    
    map<faststring, unsigned> pattern_counter;
    
    if (format == outputformat_pattern_absolute_fill || format == outputformat_pattern_relative_fill)
    {
      vector<faststring> all_patterns;
      cartesian_product  cp;
      
      if (datatype == basic_model::DNA)
        cp.operator()(N_OTUs, "ACGT", all_patterns);
      else if (datatype == basic_model::Protein)
        cp.operator()(N_OTUs, "ARNDCQEGHILKMFPSTWYV", all_patterns);
      else
      {
        cerr << "Unknown data type passed to function void print_results_simdata(..). Please report this bug.\n";
        exit(-33);
      }
      for (unsigned i=0; i<all_patterns.size(); ++i)
      {
        pattern_counter[all_patterns[i]] = 0;
      }
      cerr << "Size of all-pattern-vector: " << all_patterns.size() << '\n';
    }
    
    // Count patterns over all partitions
    // For all partitions:
    for (unsigned k=0; k<N_partitions; ++k)
    {
      unsigned num_sites;
      // Initialise aln:
      it     = OTU_data[k].begin();
      it_end = OTU_data[k].end();
      
      num_sites = (unsigned)it->second->size();
      total_sites += num_sites;
      
      int z=0;
      while (it != it_end) {
        aln[z] = it->second->begin();
        ++it;
        ++z;
      }
      
      // For all sites:
      for (unsigned i=0; i < num_sites; ++i)
      {
        for (unsigned j=0; j < N_OTUs; ++j)
        {
          pattern[j] = aln[j][i];
        }
        add_or_count_pattern(pattern_counter, pattern);
      }
    }
    cerr << "Size of map pattern_counter: " << pattern_counter.size() << '\n';
    
    // Print data:
    map<faststring, unsigned>::iterator it_pat, it_pat_end;
    it_pat     = pattern_counter.begin();
    it_pat_end = pattern_counter.end();
    
    if (format == outputformat_pattern_absolute || format == outputformat_pattern_absolute_fill)
    {
      //      os << "### Site pattern frequencies: Repetition " << rep << endl;
      while (it_pat != it_pat_end)
      {
        if (it_pat->first.size() != N_OTUs)
        {
          cerr << "WARNING: Pattern count contains pattern with string len: " << (it_pat->first).size() << '\n';
        }
        os << it_pat->first << '\t' << it_pat->second << endl;
        ++it_pat;
      }
    }
    else
    {
      while (it_pat != it_pat_end)
      {
        if ((it_pat->first).size() != N_OTUs)
        {
          cerr << "WARNING: Pattern count contains pattern with string len: " << (it_pat->first).size() << '\n';
        }
        //      os << "### Site pattern frequencies: Repetition " << rep << endl;
        os << it_pat->first << '\t' << setprecision(8) << (double)it_pat->second/(double)total_sites << endl;
        ++it_pat;
      }
    }
    
  } //END if( format == outputformat_pattern_absolute || format == outputformat_pattern_relative ||
  //        format == outputformat_pattern_absolute_fill || format == outputformat_pattern_relative_fill )
}



int main(int argc, char** argv)
{
  time_t         start_time;
  time_t         end_time;
  
  unsigned       i, rep;
  
  ofstream       os_tree;
  ofstream       os_newm;
  
  ofstream       os_ancestral_seq;
  ofstream       os_siterates_data;
  ofstream       os_siterates_histogram;
  
  ofstream       os;
  ifstream       is;
  istringstream  iss;
  //  FILE       *logfile;
  ofstream       logfile;
  ostream        *seq_os = &cout;
  ofstream       outfile;
  
  //  model_admin modelchef;
  //  BasicTree* currTree;
  
  // Simulation data: (Could be combined in a class)
  tree_admin         tree_master;
  model_admin        model_master;
  
  unsigned           alignmentLength;
  unsigned           numPartitions;
  unsigned           repetitions;
  
  const faststring   *currTreeString;
  BasicTree    *currTree;
  const basic_model  *startModelCurrTree;
  unsigned           currPartitionSize;
  faststring         *currStartsequence;
  //  faststring         dataTypeString;
  
  double             (*random_gamma)(double, double);
  double             (*random_lf_co)();
  void               (*seedrandom)(unsigned long);
  
  vector<BasicTree*>    tree_vec;
  vector<const faststring*> startSeq_vec;
  
  bool               reinit_siterates = false;
  
  
  
  
  start_time = time(NULL);
  
  welcome(stderr);
  if ( read_and_init_parameters(argc, argv) <0 ) 
  {
    cerr << "Abording due to errors." << endl;
    exit(-1);
  }
  report_seed(stderr);
  
  if (global_logging)
  {
    logfile.open(global_log_file.c_str());
    //    logfile = fopen(global_log_file.c_str(), "w");
    welcome(logfile);
    print_analysis_parameters(logfile, "");
  }
  
  if (!global_output_filename.empty() )
  {
    outfile.open(global_output_filename.c_str() );
    seq_os = &outfile;
  }
  
  random_gamma   =  &(next_random_gamma_MT19937);
  random_lf_co   =  &(next_random_lf_co_MT19937);
  seedrandom     =  &(srandom_MT19937);
  
  seedrandom(global_seed_random_generator);
  
  if (global_verbosity >= 2) // More-progress
  {
    cerr << "Reading model file." << endl;
  }
  
  try {
    model_master.create(global_model_file.c_str());
  }
  catch(basic_model::readerror x)
  {
    cerr << "Error while reading the model file: " << global_model_file << endl;
    cerr << "Line:   " << x.getLine() << endl;
    cerr << "Reason: " << x.getUnknwonKeyword() << endl;
    exit(1);
  }
  if (global_verbosity >= 2)
  {
    cerr << "Model file has been read successfully.\n" << endl;
  }
  
  if (global_verbosity >= 2)
  {
    cerr << "Reading tree file." << endl;
  }
  try {
    tree_master.create(global_tree_file.c_str());
  }
  catch (tree_admin::readerror x)
  {
    cerr << "Error while reading the tree file: " << global_tree_file << endl;
    cerr << "Line:   " << x.getLine() << endl;
    cerr << "Reason: " << x.getUnknwonKeyword() << endl;
    exit(1);
  }
  // For all trees, add the partition size
  alignmentLength = tree_master.getAlignmentLength();
  numPartitions   = tree_master.getNumTrees();
  
  if (numPartitions == 0)
  {
    cerr << "Error: no input trees found." << endl;
    exit(0);
  }
  
  if (global_verbosity >= 2)
  {
    cerr << "Tree file has been read successfully.\n" << endl;
  }
  
  if (global_verbosity >= 1)
  {
    cerr << "Initializing models." << endl;
  }
  
  model_master.init_models(random_gamma, random_lf_co);
  
  if (global_verbosity >= 1)
  {
    cerr << "Finished initializing models.\n" << endl;
  }
  
  if (global_verbosity >= 2)
  {
    cerr << endl;
    cerr << "Models that have been read from the file:" << endl;
    model_master.print(cerr);
    cerr << "Trees that have been read from the file:\n" << endl;
    tree_master.print(cerr);
  }
  
  if (global_logging)
  {
    myPrint(logfile, "Models that have been read from the model file: ", global_model_file.c_str(), "\n");
    model_master.print(logfile);
    myPrint(logfile, "Trees that have been read from the file:", global_tree_file.c_str(), "\n");
    tree_master.print(logfile, 1, "");
    myPrint(logfile, "\n");
  }
  
  if (!global_ancestral_sequence_file.empty())
  {
    os_ancestral_seq.open(global_ancestral_sequence_file.c_str());
  }
  
  
  // Tree parsing and linking of models to tree could be done once before the loop over all repetitions.
  // The loop over all repetitions should then begin with computing and setting the start sequence.
  // We dont need a new tree_vec for each replicate.
  
  // Currently, the startsequence is copied - which could be done more efficient.
  
  // loesung: tree_struct in tree_admin soll einen Zeiger auf den Baum speichern.
  
  if (!global_siterateshist_file.empty() )
  {
    os_siterates_histogram.open(global_siterateshist_file.c_str());
  }
  
  repetitions = global_num_repetitions;
  for (rep=0; rep < repetitions; ++rep) // Beginning of outer loop - all repetitions
  {
    cerr << "Simulating rep: " << rep+1 << endl;
    
    if (global_logging)
    {
      logfile << "Simulating rep: " << rep+1 << endl;
    }
    
    for (i=0; i < numPartitions; ++i) // Beginning of inner loop - all partitions
    {
      if (global_logging)
      {
        logfile << "Simulating partition: " << i+1 << endl;
      }
      if (global_verbosity >= 4 )
      {
        cerr << "Simulating partition: " << i+1 << endl;
      }
      // Currently, this loop contains code that could be moved to the tree_master:
      //    The following code manages the conversion of the tree string to the
      //    tree data structure - which does not have to be done here.
      //    Finally, the pointers to the tree structure are given to the tree_admin.
      
      
      currTreeString           = &tree_master.getTreeString(i);
      currPartitionSize        = tree_master.getPartitionSize(i);
      try
      {
        startModelCurrTree       = model_master.get_model(tree_master.getDefaultModelName(i));
      }
      catch (model_admin::indexerror x)
      {
        cerr << x.getReason_c_str() << endl;
        cerr << "Abording due to errors." << endl;
        exit(0);
      }
      
      iss.clear();
      faststring tmp_str_ = *currTreeString;
      iss.str(tmp_str_.c_str());
      if (global_verbosity >= 100)
      {
        cerr << "iss.str():                   " << iss.str() << endl;
      }
      
      if (global_logging)
      {
        logfile << "Tree read from tree file and stored string that will be parsed: " << iss.str() << endl;
      }
      
      currTree          = new BasicTree;
      currStartsequence = new faststring;
      
      startModelCurrTree->get_random_sequence(*currStartsequence, currPartitionSize);
      
      if (global_logging)
      {
        myPrint(logfile, "Start sequence: ", currStartsequence->c_str(), "\n");
      }
      if (global_verbosity >= 200)
      {
        cerr << "Start sequence: " << endl << *currStartsequence << endl;
      }
      
      if (!global_ancestral_sequence_file.empty())
      {
        // report: rep and i (partition number)
        faststring rep_str(rep);
        //	faststring part_num(i);
        if (i==0)
          os_ancestral_seq << ">Ancestral-sequence:_rep_" << rep_str << endl;
        os_ancestral_seq << *currStartsequence;
        if (i == numPartitions - 1)
          os_ancestral_seq << endl;
      }
      
      // Pass the pointer of the tree to the tree master -- could be moved to tree_master
      tree_master.set_itsBasicTree(i, currTree);
      tree_vec.push_back(currTree);
      startSeq_vec.push_back( currStartsequence );
      
      // Parse trees and link models -- could be moved to tree_master
      //    Further advantages: errors while parsing the trees are detected
      //    before the first simulations might have been done.
      try
      {
        currTree->set_read_model_status(true);
        currTree->read_tree(iss, tree_master.getScalingFactor(i));
        if (global_verbosity >= 4)
        {
          cerr << "Tree for this partition as in memory (before linking models): " << endl;
          currTree->output(cerr, 1);
          cerr << endl;
        }
        
        if (global_verbosity >= 3)
        {
          cerr << "Linking nodes of tree to models."  << endl;
        }
        currTree->link_models_to_tree(&model_master, startModelCurrTree);
        if (global_verbosity >= 3)
        {
          cerr << "Finished linking nodes of tree to models." << endl;
        }
        
        if (global_verbosity >= 4)
        {
          cerr << "Tree for this partition (after linking models):           ";
          currTree->output(cerr, 1);
          cerr << endl;
          cerr << "root model name of tree:  " << tree_master.getDefaultModelName(i) << endl << endl;
        }
      }
      catch (model_admin::indexerror x)
      {
        cerr << x.getReason() << endl;
        exit(0);
      }
      
      currTree->get_treeroot()->set_sequence( *currStartsequence );
      
      // In future versions we can save a small amount of time if we do not have to reinit the siterates
      // every time.
      // In reinit mode: We do not link inherited models again. We do not reallocate memory.
      // But be careful: For different partitions we always need to reinit siterates, since
      //                 sequence lengths are different.
      // For now, we always set reinit to false, so that new partitions sizes are always recognized!!!
      
      model_master.reset_siterates();
      
      reinit_siterates = false;
      model_master.init_siterates( currPartitionSize, reinit_siterates);
      reinit_siterates = true;  // After first initialisation we will always reinit siterates
      
      //       if (global_verbosity >= 100)
      //       {
      // 	//      model_master.print_relative_site_rates(cerr);
      // 	model_master.print_site_rates_histogramm_data(cerr);
      //       }
      
      // Print siterates histogram
      if (!global_siterateshist_file.empty() )
      {
        faststring description = "rep " + faststring(rep) + " partition " + faststring(i+1);
        model_master.print_site_rates_histogramm_data(os_siterates_histogram, description);
      }
      
      // Print detailed siterates data
      if (!global_siteratesdata_file.empty())
      {
        os_siterates_data.open(global_siteratesdata_file.c_str());
        model_master.print_relative_site_rates(os_siterates_data);
        os_siterates_data.close();
      }
      if (global_verbosity >= 4)
      {
        cerr << "Tree/partition number i before evolve: " << i << endl;
      }
      currTree->evolve_tree();
    } // End of inner loop - all partitions
    
    if (!tree_master.equal_takon_sets() )
    {
      std::cerr << "Not all trees have the same terminal taxa. Simulation has been aborded.\n"  << std::endl;
      exit(0);
    }
    
    // The map_of_OTUs is defined and filled in BasicTree.(h|cpp). It is a map of node names to sequences of the OTUs.
    // Here we create a vector of these maps, for each partition of the simulated data.
    
    
    // Collect results:
    vector<BasicTree::map_of_OTUs> OTU_data(numPartitions);
    for (i=0; i<numPartitions; ++i)
      tree_vec[i]->get_map_of_OTUs(OTU_data[i]);
    
    // Check for equal set of OTUs
    //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    
    // Print results:
    if (rep == 0)
    {
      print_results_header(*seq_os, repetitions, (unsigned)OTU_data[0].size(), tree_master, global_outputformat);
      if (global_logging)
        print_results_header(logfile, repetitions, (unsigned)OTU_data[0].size(), tree_master, outputformat_logfile);
    }
    
    // TODO: Check: For mixed models, the dataTypeString contains all models. But this is determined and used for all partitions, which might be wrong.
    faststring dataTypeString = tree_master.get_dataTypeString();
    //    int        dataType       = tree_vec[i]->get_root_model()->get_datatype();
    int dataType;
    
    if (dataTypeString == "DNA")
      dataType = 0;
    else
      dataType = 1;
    print_results_preanalysis(*seq_os, rep, (unsigned)OTU_data[0].size(), alignmentLength, dataTypeString, global_outputformat);
    print_results_simdata(*seq_os, OTU_data, rep, (unsigned)OTU_data[0].size(), alignmentLength, global_outputformat, dataType);
    print_results_postanalysis(*seq_os, rep, (unsigned)OTU_data[0].size(), alignmentLength, dataTypeString, global_outputformat);
    
    // Free memory not needed any more:
    unsigned i, N = (unsigned)tree_vec.size();
    for (i=0; i<N; ++i) {
      delete tree_vec[i];
      delete startSeq_vec[i];
    }
    OTU_data.clear();
    tree_vec.clear();
    startSeq_vec.clear();
  } // end outer loop - all repetitions
  
  if (!global_siterateshist_file.empty() )
  {
    os_siterates_histogram.close();
  }
  
  
  if (!global_ancestral_sequence_file.empty())
  {
    os_ancestral_seq.close();
  }
  
  end_time = time(NULL);
  
  cerr << "Simulation finished successfully after " << end_time - start_time
  << " seconds." << endl;
  
  if (global_logging)
  {
    logfile  << "Simulation finished successfully after " << end_time - start_time
    << " seconds." << endl;
    
    logfile.close();
  }
  
#if defined(mingw32_HOST_OS) || defined(__MINGW32__) || defined(WIN32)  
  system("pause");
#endif
  
  if (!global_output_filename.empty() )
  {
    outfile.close();
  }
  
  return 0;
}

// Needs to be revised after changes in main and other parts of the program
// TODO:
int run_internal_simulation(char   data_type,             // 'n' for nuc or 'p' for protein
                            faststring modelname_param,
                            faststring modeltype_param,
                            vector<double> *rrates_param, // Parameters are supplied as in an upper triangular matrix.
                            vector<double> *base_param,
                            double         shape_param,
                            double         pinv_param,
                            unsigned       ncat_param,
                            double         *tstv_param,
                            unsigned       seq_len,
                            unsigned       myreps,
                            faststring         tree_str,
                            unsigned       sim_seed,
                            void (*call_back_analysis)(const CSequences3 * p_seqs)
                            )
{
  UNUSED(data_type);
  UNUSED(modelname_param);
  UNUSED(modeltype_param);
  UNUSED(rrates_param);
  UNUSED(base_param);
  UNUSED(shape_param);
  UNUSED(pinv_param);
  UNUSED(ncat_param);
  UNUSED(tstv_param);
  UNUSED(myreps);
  UNUSED(tree_str);
  UNUSED(call_back_analysis);
  
  unsigned       i, rep;
  
  ofstream       os_tree;
  ofstream       os_newm;
  
  istringstream  iss;
  ofstream       logfile;
  //  ostream        *seq_os = &cout;
  //  ofstream       outfile;
  
  
  //  model_admin modelchef;
  //  BasicTree* currTree;
  
  // Simulation data: (Could be combined in a class)
  tree_admin         tree_master;
  model_admin        model_master;
  
  unsigned           alignmentLength=seq_len;
  unsigned           numPartitions=1;
  unsigned           repetitions;
  
  //  const faststring       *currTreeString;
  BasicTree    *currTree;
  const basic_model  *startModelCurrTree;
  unsigned           currPartitionSize;
  faststring             *currStartsequence;
  faststring             dataTypeString;
  
  double             (*random_gamma)(double, double);
  double             (*random_lf_co)();
  void               (*seedrandom)(unsigned long);
  
  vector<BasicTree*>    tree_vec;
  vector<const faststring*> startSeq_vec;
  
  global_seed_random_generator = sim_seed;
  
  // welcome(stderr);
  // read_and_init_parameters(argc, argv);
  report_seed(stderr);
  
  //if (global_logging)
  //{
  //  logfile.open(global_log_file.c_str() );
  //  welcome(logfile);
  //  print_search_parameters(logfile, "");
  //}
  
  random_gamma   =  &(next_random_gamma_MT19937);
  random_lf_co   =  &(next_random_lf_co_MT19937);
  seedrandom     =  &(srandom_MT19937);
  
  seedrandom(global_seed_random_generator);
  
  try {
    // model_master.append_(global_model_file.c_str());
  }
  catch(basic_model::setmodelerror x)
  {
    cerr << "Error while setting the substitution model."  << endl;
    exit(1);
  }
  try {
    tree_master.create(global_tree_file.c_str());
  }
  catch (tree_admin::readerror x)
  {
    cerr << "Error while reading the tree file: " << global_tree_file << endl;
    cerr << "Line:   " << x.getLine() << endl;
    cerr << "Reason: " << x.getUnknwonKeyword() << endl;
    exit(1);
  }
  
  model_master.init_models(random_gamma, random_lf_co);
  
  if (global_logging)
  {
    logfile << "Models read in from model file: " << global_model_file << endl;
    model_master.print(logfile);
    logfile << "Trees read in from model file: " << global_model_file << endl;
    tree_master.print(logfile);
  }
  
  DEBUGCODE( model_master.print() );
  DEBUGCODE( tree_master.print()  );
  
  // For all trees, add the partition size
  alignmentLength = tree_master.getAlignmentLength();
  numPartitions   = tree_master.getNumTrees();
  
  if (numPartitions == 0)
  {
    cerr << "Error: no input trees found." << endl;
    exit(0);
  }
  
  // Tree parsing and linking of models to tree could be done once before the loop over all repetitions.
  // The loop over all repetitions should then begin with computing and setting the start sequence.
  // We dont need a new tree_vec for each replicate.
  
  // Currently, the startsequence is copied - which could be done more efficient.
  
  // Loesung: tree_struct in tree_admin soll einen Zeiger auf den Baum speichern.
  
  repetitions = global_num_repetitions;
  for (rep=0; rep < repetitions; ++rep)
  {
    for (i=0; i < numPartitions; ++i)
    {
      // Currently, this loop contains code that could be moved to the tree_master:
      //    The following code manages the conversion of the tree string to the
      //    tree data structure - which does not have to be done here.
      //    Finally, the pointers to the tree structure are given to the tree_admin.
      
      
      //        currTreeString           = &tree_master.getTreeString(i);
      currPartitionSize        = tree_master.getPartitionSize(i);
      startModelCurrTree       = model_master.get_model(tree_master.getDefaultModelName(i));
      
      //iss.clear();
      //iss.str(currTreeString->c_str());
      //DEBUGCODE(     cerr << "iss.str():                   " << iss.str() << endl);
      
      currTree          = new BasicTree;
      currStartsequence = new faststring;
      
      startModelCurrTree->get_random_sequence(*currStartsequence, currPartitionSize);
      
      if (global_logging)
      {
        logfile << "Start sequence: " << *currStartsequence << endl;
      }
      DEBUGCODE(   cerr << "Start sequence: " << *currStartsequence << endl;   );
      // Pass the pointer of the tree to the tree master -- could be moved to tree_master
      tree_master.set_itsBasicTree(i, currTree);
      tree_vec.push_back(currTree);
      startSeq_vec.push_back( currStartsequence );
      
      // Parse trees and link models -- could be moved to tree_master
      //    Further advantages: errors while parsing the trees are detected
      //    before the first simulations might have been done.
      try {
        currTree->set_read_model_status(true);
        currTree->read_tree(iss, tree_master.getScalingFactor(i));
        DEBUGCODE( { cerr << "Tree as in memory (no link): "; currTree->output(cerr, 1); cerr << endl; } );
        currTree->link_models_to_tree(&model_master, startModelCurrTree);
        DDEBUGCODE( { cerr << "Tree as in memory:           "; currTree->output(cerr, 1); cerr << endl; } );
        DDEBUGCODE(   cerr << "default model name of tree:  " << tree_master.getDefaultModelName(i) << endl );
      }
      catch (model_admin::indexerror x) {
        cerr << x.getReason() << endl;
        exit(0);
      }
      
      currTree->get_treeroot()->set_sequence( *currStartsequence );
      model_master.init_siterates( currPartitionSize, false);
      
      DEBUGCODE( cerr << "Tree/partition number i before evolve: " << i << endl );
      currTree->evolve_tree();
    }
    
    if (!tree_master.equal_takon_sets() )
    {
      std::cerr << "Not all trees have the same terminal taxa. Simulation has been aborded.\n"  << std::endl;
      exit(0);
    }
    
    // Collect results:
    vector<BasicTree::map_of_OTUs> OTU_data(numPartitions);
    for (i=0; i<numPartitions; ++i)
      tree_vec[i]->get_map_of_OTUs(OTU_data[i]);
    
    // Check for equal set of OTUs
    //xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    
    // Print results:
    if (rep == 0)
    {
      print_results_header(cout, repetitions, (unsigned)OTU_data[0].size(), tree_master, global_outputformat);
      if (global_logging)
        print_results_header(logfile, repetitions, (unsigned)OTU_data[0].size(), tree_master, outputformat_logfile);
    }
    
    dataTypeString = tree_master.get_dataTypeString();
    int        dataType       = tree_vec[i]->get_root_model()->get_datatype();
    
    print_results_preanalysis(cout, rep, (unsigned)OTU_data[0].size(), alignmentLength, dataTypeString, global_outputformat);
    print_results_simdata(cout, OTU_data, rep, (unsigned)OTU_data[0].size(), alignmentLength, global_outputformat, dataType);
    print_results_postanalysis(cout, rep, (unsigned)OTU_data[0].size(), alignmentLength, dataTypeString, global_outputformat);
    
    // Free memory not needed any more:
    unsigned i, N = (unsigned)tree_vec.size();
    for (i=0; i<N; ++i) {
      delete tree_vec[i];
      delete startSeq_vec[i];
    }
    OTU_data.clear();
    tree_vec.clear();
    startSeq_vec.clear();
  } // end: All repetitions
  
  if (global_logging)
    logfile.close();
  
  return 0;
}