pvarmain.cpp

/***************************************************************************************************
*  The PolyMoSim project is distributed under the following license:
*  
*  Copyright (c) 2006-2022, Christoph Mayer, Forschungsmuseum Alexander Koenig, 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 "pvartree.h"

#include "tree_admin.h"

//#include "sequence_gen.h"

#include "CRandom.h"

#include <sstream>

#include <ctime>



using namespace std;



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



void welcome(FILE *of, const char *s="")

{

  myPrint(of, s, welcome_str);

}





void welcome(ofstream &of, const char *s="")

{

  of << s << welcome_str;

}





void print_results_header(ostream& os, unsigned repetitions, unsigned ntaxa,

			  tree_admin &tchef, output_format format)

{

  if (format == 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 == fasta)  {

  }

  else if (format == phylip) {

  }

  else if (format == log_format)

  {

    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, unsigned 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 == nexus) {

    os << "Begin data;" << endl

       << "   Dimensions ntax=" << ntaxa << " nchar=" << nchar << ";" << endl

       << "   Format missing=? gap=- datatype=DNA;" << endl

       << "   Matrix" << endl;

  }

  else if (format == fasta) {

  }

  else if (format == phylip) {

    os << "\t" << ntaxa << "\t" << nchar << endl;

  }

}





void print_results_postanalysis(ostream& os, unsigned rep, unsigned ntaxa, unsigned nchar, unsigned format)

{

  UNUSED (rep);

  UNUSED (ntaxa);

  UNUSED (nchar);



  char c;

  if (format == nexus) {

    os << ";"          << endl;

    os << "end;"  << endl << endl;

  }

  else if (format == fasta) {

  }

  else if (format == phylip) {

  }

  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, unsigned format)

{

  UNUSED (os);

  UNUSED (rep);

  UNUSED (ntaxa);

  UNUSED (nchar);



  unsigned i;

  BasicTree::map_of_OTUs::iterator  it, it_end;

  it     = OTU_data[0].begin();

  it_end = OTU_data[0].end();

  

  if (format == nexus) {

    while (it != it_end) {

      cout << *(it->first) << " ";

      for (i=0; i < OTU_data.size(); ++i)

	cout << *OTU_data[i][it->first] << " ";

      cout << endl;

      ++it;

    }

  }

  else if (format == fasta) {

    while (it != it_end) {

      cout << ">" << *(it->first) << endl;

      for (i=0; i<OTU_data.size(); ++i)

	cout << *OTU_data[i][it->first] << " ";

      cout << endl;

      ++it;

    }

  } else if (format == phylip) {

    while (it != it_end) {

      cout << *(it->first) << " ";

      for (i=0; i<OTU_data.size(); ++i)

	cout << *OTU_data[i][it->first] << " ";

      cout << endl;

      ++it;

    }

  }

}



//******************************************

//pvartree

//******************************************



int main(){



  int                     i;



  ofstream                os_tree;

  ofstream                os_newm;



  ofstream                os;

  ifstream                is;

  istringstream           iss;

  model_admin             modelchef;

  BasicTree               baum;

  //  BasicTree               *stringtree;

  tree_admin              tree_master;

  model_admin             model_master;



  faststring                  rootsequence;

  //  double                  A, G, C;



  double                  (*normal_random)(double);

  void                    (*seedrandom)(unsigned long);



  normal_random   =       &(next_random_gauss_polar_box_muelle_MT19937);

  seedrandom      =       &(srandom_MT19937);



  seedrandom(911);



  

  //  sequence_gen            startsequence;



  if (!FileExists("Modell.dat"))

  {

    cerr << "File Modell.dat does not exist." << endl;

    exit(1);

  }

  if (!FileExists("Baum.dat"))

  {

    cerr << "File Baum.dat does not exist." << endl;

    exit(1);

  }

  if (!FileExists("parameter.dat"))

  {

    cerr << "File parameter.dat does not exist." << endl;

    exit(1);

  }



  try {

    model_master.create("Modell.dat");

  }

  catch(nuc_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);

  }

  try {

    tree_master.create("Baum.dat");

  }

  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);

  }





  os_tree.open("pvartrees.dat");

  os_newm.open("new_Models.dat");



  faststring dummy;

  faststring tmp_s;



  for(i = 0; i < tree_master.getNumTrees(); i++) {

    pvartree tree_to_vary;

    iss.clear();

    tmp_s = tree_master.getTreeString(i);



    iss.str(tmp_s.c_str());



    tree_to_vary.set_random_generator(normal_random);

    try{

      tree_to_vary.read_parameters_from_file("parameter.dat");

    }

    catch(pvartree::error x) {

      cerr << "Error: File " <<  "parameter.dat" << " not found." << endl;

      exit(1);

    }

    catch(pvartree::readerror x) {

      cerr << "Error while reading the parameterfile: " << "parameter.dat" << endl;

      cerr << "Line:   " << x.getLine() << endl;

      cerr << "Reason: " << x.getUnknownKeyword() << endl;

      exit(1);

    }



    dummy = "Model_default_";

    dummy.append(faststring(i+1));

       

    tree_to_vary.vary_tree(i+1, iss, "Modell.dat", tree_master.getScalingFactor(i), tree_master.getDefaultModelName(i));

    os_tree.width(10);

    os_tree << tree_master.getScalingFactor(i);

    os_tree.width(10);

    os_tree << tree_master.getPartitionSize(i);

    os_tree.width(25);

    os_tree << dummy;

    os_tree << "   ";

    tree_to_vary.print_varied_tree(os_tree);

    os_tree << endl;



    tree_to_vary.print_models(os_newm);

  }



  os_tree.close();

  os_newm.close();





#if defined(mingw32_HOST_OS) || defined(__MINGW32__) || defined(WIN32)  

  system("pause");

#endif



  return 0;

}