/*  Copyright by Matthias Hoechsmann (C) 2002-2004  =====================================                                     You may use, copy and distribute this file freely as long as you  - do not change the file,  - leave this copyright notice in the file,  - do not make any profit with the distribution of this file  - give credit where credit is due  You are not allowed to copy or distribute this file otherwise  The commercial usage and distribution of this file is prohibited  Please report bugs and suggestions to <mhoechsm@TechFak.Uni-Bielefeld.DE>*/#include <deque>#include <functional>#include <fstream>#include <iomanip>#include <iostream>#include <list>#include <sstream>#include <string>#include <map>//#include <sys/timeb.h>
#include <sys/times.h>
#include <unistd.h>
#ifndef WIN32#include "config.h"#endif#include "Arguments.h"#include "alignment.h"#include "debug.h"
//#include "global_alignment.h"
#include "treeedit.h"#include "misc.h"
#include "progressive_align.h"
#include "rna_alignment.h"#include "rnaforest.h"
#include "rnaforestsz.h"#include "rnafuncs.h"#include "rna_profile_alignment.h"#include "rna_algebra.h"#include "rnaforester_options.h"#include "alignment.t.cpp"
//#include "global_alignment.t.cpp"
#include "treeedit.t.cpp"//#include "ppforest.t.cpp"using namespace std;/* ****************************************** *//*          Definitions and typedefs          *//* ****************************************** */struct ToLower : public unary_function<char,char> {	char operator()(char a)	{		return tolower(a);	};};template <class L>void makeDotFileAli(const PPForest<L> &ppf, const RNAforesterOptions &options){	if(options.has(RNAforesterOptions::OutputAlignmentDotFormat))	{		string filename;		options.get(RNAforesterOptions::OutputAlignmentDotFormat,filename,string("ali.dot"));		ofstream s(filename.c_str());		ppf.printDot(s);	}}template <class L>void makeDotFileInp(const PPForest<L> &ppf, const RNAforesterOptions &options, Uint count){	if(options.has(RNAforesterOptions::MakeDotForInputTrees))	{		ostringstream ss;		ofstream os;		ss << "input" << count << ".dot";		os.open(ss.str().c_str());		ppf.printDot(os);		os.close();	}}static const string RNAFORESTER_VERSION = "1.5";static const string PROMPT = "Input string (upper or lower case); & to end for multiple alignments, @ to quit\n";static const string SCALE = "....,....1....,....2....,....3....,....4....,....5....,....6....,....7....,....8\n";void alignMultiple(deque<RNAProfileAlignment*> &alignList, Score &score,const RNAforesterOptions &options);void alignPairwise(deque<RNAForest*> &inputListPW,Score &score,const RNAforesterOptions &options);void cutAfterChar(string &s,char c);void editPairwise(list<RNAForestSZ*> &inputListSZ,Score &score,RNAforesterOptions &options);void alignPairwiseSimple(deque<RNAForest*> &inputListPW,Score &score,RNAforesterOptions &options);
static void showversion(const char *prog){	cout << prog << ", version " << RNAFORESTER_VERSION << endl;	cout << "Copyright Matthias Hoechsmann 2001-2004," << endl << "mhoechsm@techfak.uni-bielefeld.de" << endl;}int main(int argc, const char **argv){		string buffer;	string baseStr,viennaStr,nameStr;	Ulong basePairCount,maxDepth;	deque<RNAForest*> inputListPW;	deque<RNAProfileAlignment*> alignList;	bool showScale=true,multipleAlign=false;	istream *inputStream=NULL;	ifstream *inputFile=NULL;	Uint structure_count=1;	int suboptPercent=100;	double minPairProb=0.25;	list<RNAForestSZ*> inputListSZ;	try	{		RNAforesterOptions options(argc,argv);		// check options		if(options.has(RNAforesterOptions::Help))		{			options.help();			exit(EXIT_SUCCESS);		}		if(options.has(RNAforesterOptions::SecretHelp))		{			options.secretHelp();			exit(EXIT_SUCCESS);		}		if(options.has(RNAforesterOptions::Version))		{			showversion(argv[0]);			exit(EXIT_SUCCESS);      		}		// read score values		Score score(options);		if(!options.has(RNAforesterOptions::ShowOnlyScore))			score.print();		// check option suboptimals		if(options.has(RNAforesterOptions::LocalSubopts))		{			options.get(RNAforesterOptions::LocalSubopts,suboptPercent,100);			if(suboptPercent<0 || suboptPercent>100)			{				cerr << "error: value for parameter --subopts must be in range from 0 to 100" << endl;				exit(EXIT_FAILURE);			}		}#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library				options.get(RNAforesterOptions::PredictMinPairProb,minPairProb,0.25);		if(options.has(RNAforesterOptions::PredictProfile))		  cout << "Minumim required basepair probability (-pmin): " << minPairProb << endl;#endif				// show if suboptimals		if(!options.has(RNAforesterOptions::ShowOnlyScore))			if(options.has(RNAforesterOptions::LocalSimilarity) && options.has(RNAforesterOptions::LocalSubopts))			  cout << "calculate suboptimals within " << suboptPercent << "% of global optimum" << endl << endl;		// profile search		if(options.has(RNAforesterOptions::ProfileSearch))		  {		    string filename;		    		    options.get(RNAforesterOptions::ProfileSearch,filename,string(""));		    if(filename=="")		      {			cerr << "no profile filename" << endl;			exit(EXIT_FAILURE);		      }		    		    RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(filename);		    alignList.push_back(rnaProfileAli);		  }		if(options.has(RNAforesterOptions::ReadFromFile))		{			string filename;			options.get(RNAforesterOptions::ReadFromFile,filename,string(""));			inputFile=new ifstream(filename.c_str());			if(inputFile->fail())			{				cerr << "cannot open file: \"" << filename << "\"" << endl;				exit(EXIT_FAILURE);			}			inputStream=inputFile;		}		else
		{			inputStream=&cin;
		}

		if(showScale)		  {		    if(!options.has(RNAforesterOptions::NoScale) && !options.has(RNAforesterOptions::ReadFromFile))		      cout << endl << PROMPT << SCALE;		    showScale=false;		  }		for(;;)		{			getline(*inputStream,buffer);			if(inputStream->eof())			  {			    			    if(options.has(RNAforesterOptions::Multiple) && !options.has(RNAforesterOptions::ProfileSearch))				buffer="&";			    else			      exit(EXIT_SUCCESS);			  }			if(buffer.empty())				continue;			// quit if character is @			if(buffer[0]=='@')				break;			// delete '\r' at line end from non unix files			if(buffer[buffer.size()-1]=='\r')				buffer.erase(buffer.size()-1);			// check for name of structure			if(buffer[0]=='>')			{				nameStr=&buffer[1];				continue;			}			// cut after blank			cutAfterChar(buffer,' ');			// check for aligning multiple structures			// if input is read from file the eof has the same meaning as &			if( buffer[0]=='&')				multipleAlign=true;			else			{				// check for base string			  			  if(RNAFuncs::isRNAString(buffer))				{				  baseStr=buffer;				  // convert to small letters				  transform(baseStr.begin(),baseStr.end(),baseStr.begin(),ToLower());				  // t -> u				  replace(baseStr.begin(),baseStr.end(),'t','u');				  // delete '.'  and '-' from alignment files				  remove(baseStr.begin(),baseStr.end(),'.');				  remove(baseStr.begin(),baseStr.end(),'-');#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library				  				  if(options.has(RNAforesterOptions::PredictProfile))				    {				      ostringstream ss;		 				      string constraint;				      // if there is no name given (> ...) use a counter				      if(nameStr=="")					ss << "> " << structure_count;				      else					ss << nameStr;					      				      cout << "Predicting structure profile for sequence: " << ss.str() << endl;				      RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(baseStr,ss.str(),constraint,minPairProb);				      alignList.push_back(rnaProfileAli);				      makeDotFileInp(*rnaProfileAli,options,structure_count);				      structure_count++;				    }#endif				  				  //				  continue;				}			  else			    {				// check for vienna string				  if(RNAFuncs::isViennaString(buffer,basePairCount,maxDepth))			    {#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library			      			      // skip structure lines if structures are predicted			      if(options.has(RNAforesterOptions::PredictProfile))				{				  cout << "ignoring structure: " << buffer << endl;				  continue;				}#endif			      			     viennaStr=buffer;	       			    }			  else   				{ 					cerr << "The input sequence is neither an RNA/DNA string nor in vienna format." << endl;
					cerr << "line: " << buffer << endl;					showScale=true;					exit(EXIT_FAILURE);				}				// add structure to input list				if(options.has(RNAforesterOptions::Multiple))				{					ostringstream ss;		 					// if there is no name given (> ...) use a counter					if(nameStr=="")					  ss << "> " << structure_count;					else					  ss << nameStr;						RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(baseStr,viennaStr,ss.str());					makeDotFileInp(*rnaProfileAli,options,structure_count);										alignList.push_back(rnaProfileAli); 					}				else				{
				  if(options.has(RNAforesterOptions::TreeEdit))				    {				      RNAForestSZ *rnaForestSZ=new RNAForestSZ(baseStr,viennaStr,nameStr);				      inputListSZ.push_back(rnaForestSZ);				    }				  else				    {				      RNAForest *rnaForest=new RNAForest(baseStr,viennaStr,nameStr);							      nameStr="";				      makeDotFileInp(*rnaForest,options,structure_count);				      inputListPW.push_back(rnaForest);				    }				}				structure_count++;      				showScale=true;			}			}			// ***** multiple alignment			if((options.has(RNAforesterOptions::Multiple) && multipleAlign) || (options.has(RNAforesterOptions::ProfileSearch) && alignList.size()==2))			{				alignMultiple(alignList,score,options);				multipleAlign=false;				structure_count=1;								if(options.has(RNAforesterOptions::ProfileSearch))				  {				    string filename;		    				    options.get(RNAforesterOptions::ProfileSearch,filename,string(""));		    		    				    RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(filename);				    alignList.push_back(rnaProfileAli);				  }				else				  break;				//				break;			}			// ***** pairwise alignment			if(inputListPW.size()==2)
			{			  if(options.has(RNAforesterOptions::GlobalAlignment))			    alignPairwiseSimple(inputListPW,score,options);			  else			    alignPairwise(inputListPW,score,options);			  break;			}						if(inputListSZ.size()==2)			{			  editPairwise(inputListSZ,score,options);			  break;			}		}		// free dynamic allocated memory		deque<RNAForest*>::const_iterator it;		for(it = inputListPW.begin(); it!=inputListPW.end(); it++)			delete *it;		DELETE(inputFile);		//	  getchar();	// only for testing		return (0);	}	catch(RNAforesterOptions::IncompatibleException e)	{		e.showError();		return(EXIT_FAILURE);	} 	catch(RNAforesterOptions::RequiresException e)	{		e.showError();		return(EXIT_FAILURE);	} }void cutAfterChar(string &s,char c){	string::size_type pos=s.find(c);	if(pos!=string::npos)		s.erase(pos);}void alignMultiple(deque<RNAProfileAlignment*> &alignList, Score &score,const RNAforesterOptions &options){	DoubleScoreProfileAlgebraType *alg;	deque<pair<double,RNAProfileAlignment*> > resultList;//	double optScore;	Uint clusterNr=1;	double minPairProb;	options.get(RNAforesterOptions::ConsensusMinPairProb,minPairProb,0.5);		// distance or similarity	if(options.has(RNAforesterOptions::CalculateDistance))		alg=new DoubleDistProfileAlgebra(score);	else		alg=new DoubleSimiProfileAlgebra(score);	cout << endl;		progressiveAlign(alignList,resultList,alg,options);	cout << endl;	cout << "*** Results ***" << endl << endl; 	cout << "Minimum basepair probability for consensus structure (-cmin): " << minPairProb << endl << endl;	deque<pair<double,RNAProfileAlignment*> >::const_iterator it;	for(it=resultList.begin();it!=resultList.end();it++)	  {	    cout << "RNA Structure Cluster Nr: " << clusterNr << endl;	    cout << "Score: " << it->first << endl;	    cout << "Members: " << it->second->getNumStructures() << endl << endl;	    if(options.has(RNAforesterOptions::FastaOutput))