/*
  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>
#include <string.h>

#ifndef WIN32
#include "config.h"
#endif

#include "Arguments.h"
#include "alignment.h"
#include "debug.h"

//#include "global_alignment.h"
#include "treeedit.h"

#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
#include <libxml++/libxml++.h>
#include <libxml/xmlschemas.h>
#include "xsd.h"
#endif
#endif

#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,RNAFuncs::AddXmlInfos &xmlInfos);
void alignPairwise(deque<RNAForest*> &inputListPW,Score &score,const RNAforesterOptions &options,RNAFuncs::AddXmlInfos &xmlInfos);
void cutAfterChar(string &s,char c);

void editPairwise(list<RNAForestSZ*> &inputListSZ,Score &score,RNAforesterOptions &options);
void alignPairwiseSimple(deque<RNAForest*> &inputListPW,Score &score,RNAforesterOptions &options);

#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
extern "C" {
  bool validateXSD(string filename);
}
#endif
#endif

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;
	istringstream inputString;
	string tmpInput = "";
	ifstream *inputFile=NULL;
	Uint structure_count=1;
	int suboptPercent=100;
	double minPairProb=0.25;
#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
	xmlpp::Document* xmlDoc;
#endif
#endif
	string xmlOrig;
	map<int,string> comments;
	map<int,string> idMapping;
	map<int,string> descriptions;
	map<int,string> names;
	map<int,string> synonyms;
	int seqID = 0;
	stringstream ss;
	RNAFuncs::AddXmlInfos xmlInfos;
	
	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 << "Minimum 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);
			}
#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2 							
			// getline extracts all characters of the first line until '\n' to check whether inputFile is an xml file
			getline(*inputFile,buffer);
			
			// input file is an xml file 
			if(buffer.find("<?xml",0)==0){
			  buffer="";
			  
			  // validation
			  //if(!validateXSD(filename)){
			  //  exit(EXIT_FAILURE);
			  //}
			  
			  // create Dom parser
			  xmlpp::DomParser domParser;
			  domParser.parse_file(filename);
			 			  
			  xmlDoc = domParser.get_document();
			  xmlpp::Element* elemRoot = xmlDoc->get_root_node();
			  xmlpp::Node::NodeList nodeListRnaStructure = elemRoot->get_children();
			  
			  // for each rnastructure element
			  xmlpp::Node::NodeList::iterator it1;
			  xmlpp::Node::NodeList::iterator it2;
			  xmlpp::Node::NodeList::iterator it3;
			  for(it1=nodeListRnaStructure.begin();it1!=nodeListRnaStructure.end();it1++){
			    if((*it1)->get_name() == "rnastructure"){
			      seqID++;
			      tmpInput.append(">");
			      ss << seqID;
			      tmpInput.append(ss.str());
			      ss.str("");
			      tmpInput.append("\n");
			    }
			   			       
			    xmlpp::Node::NodeList nodeListSeq = (*it1)->get_children();
			    for(it2=nodeListSeq.begin();it2!=nodeListSeq.end();it2++){
			      if((*it2)->get_name()=="sequence"){
				xmlpp::Element* elemSeq = (xmlpp::Element*)(*it2);
				
				// map internal id to seqID :: JK debug
				cout << "seqID: " << seqID << endl;
				cout << "idMapping[seqID]: " << elemSeq->get_attribute("seqID")->get_value() << endl;
				//
				idMapping[seqID] = elemSeq->get_attribute("seqID")->get_value();

				xmlpp::Node::NodeList nodeListName = elemSeq->get_children();
				for(it3=nodeListName.begin();it3!=nodeListName.end();it3++){
				  if((*it3)->get_name()=="name" && ((xmlpp::Element*)(*it3))->get_child_text()){
				    // map the name to seqID
				    names[seqID] = ((xmlpp::Element*)(*it3))->get_child_text()->get_content();
				  }
				  if((*it3)->get_name()=="synonyms" && ((xmlpp::Element*)(*it3))->get_child_text()){
				    // map synonyms to seqID
				    synonyms[seqID] = ((xmlpp::Element*)(*it3))->get_child_text()->get_content();
				  }
				  if((*it3)->get_name()=="description" && ((xmlpp::Element*)(*it3))->get_child_text()){
				    // map the description to seqID
				    descriptions[seqID] = ((xmlpp::Element*)(*it3))->get_child_text()->get_content();
				  }
				  if((*it3)->get_name()=="freeSequence" || (*it3)->get_name()=="nucleicAcidSequence"){
				    xmlpp::Element* elemFreeSequence = (xmlpp::Element*)(*it3);
				    tmpInput.append(elemFreeSequence->get_child_text()->get_content());
				    tmpInput.append("\n");
				  }
				} // end for it3
			      }
			      
			      // get comment if available
			      if((*it2)->get_name()=="comment" && ((xmlpp::Element*)(*it2))->get_child_text()){
				comments[seqID] = ((xmlpp::Element*)(*it2))->get_child_text()->get_content();
			      }
			      			    			      
			      // get structure
			      if((*it2)->get_name()=="structure" && ((xmlpp::Element*)(*it2))->get_child_text()){
				xmlpp::Element* elemStr = (xmlpp::Element*)(*it2);
				tmpInput.append(elemStr->get_child_text()->get_content());
				tmpInput.append("\n");
			      }
			    } // end for it2
			    
			  } // end for it1
			  xmlOrig = xmlDoc->write_to_string();

			  // struct with additional xml infos
			  xmlInfos.idmapping = idMapping;
			  xmlInfos.comments = comments;
			  xmlInfos.descriptions = descriptions;
			  xmlInfos.names = names;
			  xmlInfos.synonyms = synonyms;
			  xmlInfos.xmlInput = true;
				  
			  inputString.str(tmpInput.c_str());
			  inputStream=&inputString;
			  
			} else {

			  xmlInfos.xmlInput = false;
			  // no xml file write the first line from buffer back to stream
			  for(int i=buffer.size();i>=0;i--){
			    inputFile->putback(buffer[i]);
			  }
#endif
#endif
			  inputStream=inputFile;

#ifdef HAVE_LIBXMLPLUSPLUS	
#ifdef HAVE_LIBXML2	  
			}
#endif
#endif
		} 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]=='>')
			{
			  /*if(buffer.find(" ",0) != string::npos){
			    int t = buffer.find(" ",0);
			    nameStr=&buffer[t];
			  } else {
			    nameStr=&buffer[1];
			    }*/
			  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,xmlInfos);
				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,xmlInfos);
			  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);

		return (0);
	}
	catch(RNAforesterOptions::IncompatibleException e)
	{
		e.showError();
		return(EXIT_FAILURE);
	} 
	catch(RNAforesterOptions::RequiresException e)
	{
		e.showError();
		return(EXIT_FAILURE);
	}
#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
	catch(xmlpp::validity_error ve) 
	{
	  cout << ve.what() << endl;
	//return(EXIT_FAILURE);
	} 
	catch(xmlpp::parse_error pe)
	{
	  cout << pe.what() << endl;
	  //return(EXIT_FAILURE);
	} 
	catch(xmlpp::exception e)
	{
	  cout << e.what() << endl;
	  //return(EXIT_FAILURE);
	} 
#endif
#endif
	  	
}



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,RNAFuncs::AddXmlInfos &xmlInfos)
{
	DoubleScoreProfileAlgebraType *alg;
	deque<pair<double,RNAProfileAlignment*> > resultList;
//	double optScore;
	Uint clusterNr=1;
	double minPairProb;
	string outputFile;

	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))
	      {
		it->second->printFastaAli(false);
		cout << endl;
	      }
	    else
	      {  

		it->second->printSeqAli();

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library		
		// print alignment
		it->second->printStrAli();
		cout << endl;
#endif		
	      }

	    // save profile
	    if(options.has(RNAforesterOptions::SaveProfile))
	    {
	      string filename;
	      filename=options.generateFilename(RNAforesterOptions::SaveProfile,".pro", "rna.pro",clusterNr);
      it->second->save(filename);
	    }
	    

		// print consensus structure
		cout << "Consensus sequence/structure:" << endl;
		it->second->printConsensus(minPairProb);
		cout << endl;

#ifdef HAVE_LIBG2  // This features require the g2 library		
		// generate squiggle plots
			if(options.has(RNAforesterOptions::MakeSquigglePlot))
			{
				RNAProfileAlignment::SquigglePlotOptions sqOptions;
				string filename;

				// plot showing full base information
				filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,".ps", "rnaprofile.ps",clusterNr);
				sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
				sqOptions.minPairProb=minPairProb;
				sqOptions.mostLikelySequence=false;
				it->second->squigglePlot(filename,sqOptions);

				// plot showing consensus sequence
				filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,"_cons.ps", "rnaprofile_cs.ps",clusterNr);
				sqOptions.mostLikelySequence=true;
				it->second->squigglePlot(filename,sqOptions);
			}
#endif		      

		clusterNr++;
	  }
#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
	if(options.has(RNAforesterOptions::XmlOutputFile)){
	  options.get(RNAforesterOptions::XmlOutputFile,outputFile,string(""));
	} else {
	  outputFile = "result.xml";
	}

	if(options.has(RNAforesterOptions::GenerateXML))
        { 
		RNAFuncs::printMAliXML(resultList,options,minPairProb,xmlInfos,outputFile);	
	}
#endif
#endif

	/*
	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
		cout << "global optimal score: ";
	}
	cout << optScore << endl;

	list<RNAProfileAlignment*>::iterator it=inputMapProfile.begin();
	RNAProfileAlignment *ppf=it;
	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
		// generate dot file
		makeDotFileAli(*ppf,options);

		ppf->print();
		cout << endl;
		ppf->printConsensus();
	}
	*/

	// save profile alignment to binary file
	/*	      if(options.has(RNAforesterOptions::SaveMultipleAliFile))
	{
	string filename;

	filename=generateFilename(options,RNAforesterOptions::SaveMultipleAliFile,".mta", "unknown.dot");
	ofstream s(filename.c_str());
	f1->save(s);
	}
	*/

	/*
	// generate squiggle plots
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{
		RNAProfileAlignment::SquigglePlotOptions sqOptions;
		string filename;

		// plot showing full base information
		filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,".ps", "rnaprofile.ps");
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		sqOptions.mostLikelySequence=false;
		ppf->squigglePlot(filename,sqOptions);

		// plot showing consensus sequence
		filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,"_cons.ps", "rnaprofile_cons.ps");
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		sqOptions.mostLikelySequence=true;
		ppf->squigglePlot(filename,sqOptions);
	}
	*/

	
	delete alg;
}


void alignPairwise(deque<RNAForest*> &inputListPW,Score &score,const RNAforesterOptions &options,RNAFuncs::AddXmlInfos &xmlInfos)
{
	IntScoreRNA_AlgebraType *alg;
	Uint xbasepos,ybasepos,xlen,ylen;
	list<pair<Uint,Uint> > xsubopts;
	list<pair<Uint,Uint> > ysubopts;
	string seq1,seq2,str1,str2;
	char s[8];
	int suboptPercent;
	Uint count=1;
	RNAFuncs::SquigglePlotOptions sqOptions;
	tms tmsStart, tmsEnd;
	string outputFile;

	// read options
	options.get(RNAforesterOptions::LocalSubopts,suboptPercent,100);

#ifdef HAVE_LIBG2  // This features require the g2 library	
	// generate squiggle plot
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{	      
		// get sq options
		sqOptions.hideBaseNumbers=options.has(RNAforesterOptions::SquiggleHideBaseNumbers);
		options.get(RNAforesterOptions::SquiggleBaseNumberInterval,sqOptions.baseNumInterval,(Uint)20);
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		options.get(RNAforesterOptions::SquiggleScaleFactor,sqOptions.scale,1.0);
		sqOptions.generateFIG=options.has(RNAforesterOptions::SquiggleGenerateFIG);		
#ifdef HAVE_LIBGD
		sqOptions.generatePNG=options.has(RNAforesterOptions::SquiggleGeneratePNG);
		sqOptions.generateJPG=options.has(RNAforesterOptions::SquiggleGenerateJPG);
#endif
		
	}
#endif
	
	// distance or similarity
	if(options.has(RNAforesterOptions::CalculateDistance))
		alg=new IntDistRNA_Algebra(score);
	else
	{
		if(options.has(RNAforesterOptions::RIBOSUMScore))
			alg=new RIBOSUM8560(score);
		else
			alg=new IntSimiRNA_Algebra(score);
	}

	RNAForest *f1=inputListPW.front();
	RNAForest *f2=inputListPW.back();

	if(options.has(RNAforesterOptions::SpaceTimeInfo))
	  {
	    IntScore_AlgebraType *algGlobClassic;

	    algGlobClassic=new IntDist_Algebra(score);

	    if(!options.has(RNAforesterOptions::ShowOnlyScore))
	      {
		cout << "F1_NUMNODES" << ";";
		cout << "F2_NUMNODES" << ";";
		cout << "F1_DEGREE" << ";";
		cout << "F2_DEGREE" << ";";
		cout << "F1_LEAVES" << ";"; 
		cout << "F2_LEAVES" << ";";
		cout << "F1_DEPTH" << ";"; 
		cout << "F2_DEPTH" << ";";
		cout << "F1_NUMCSFS" << ";";
		cout << "F2_NUMCSFS" << ";";
		cout << "TABLE_SIZE_4D" << ";";
		cout << "TABLE_SIZE_2D" << ";";
		cout << "GLOBALI_TIME" << ";";
		cout << "GLOBALI_TIME_SPEEDUP" << ";";
		cout << "LOCALALI_TIME" << ";";
		cout << "LOCALALI_TIME_SPEEDUP" << ";";
		
		cout << endl;
	      }

	    cout << f1->size() << "\t";
	    cout << f2->size() << "\t"; 
	    cout << f1->maxDegree() << "\t";
	    cout << f2->maxDegree() << "\t";
	    cout << f1->numLeaves() << "\t";
            cout << f2->numLeaves() << "\t";
	    cout << f1->maxDepth() << "\t";
	    cout << f2->maxDepth() << "\t";
	    cout << f1->getNumCSFs() << "\t";
	    cout << f2->getNumCSFs() << "\t";
	    cout << f1->size()*f2->size()*f1->maxDegree()*f1->maxDegree() << "\t";
	    cout << f1->getNumCSFs()*f2->getNumCSFs() << "\t";

 	    // global alignment
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,false,true);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getGlobalOptimum() << "#";
	    }

 	    // global alignment speedup
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,false,false);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getGlobalOptimum() << "#";
	    }

 	    // local alignment
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,true,true);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getLocalOptimum() << "#";
	    }

 	    // local alignment speedup
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,true,false);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getLocalOptimum() << "#";
	    }
	    

	    cout << endl;

	    exit(EXIT_SUCCESS);
	  }
	
	Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg);
	RNA_Alignment ppfali;
	ppfali.setStructureNames(f1->getName(),f2->getName());
	double optScore;


	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
	  if(options.has(RNAforesterOptions::SmallInLarge))
	    {
	      cout << "small-in-large optimal score: ";
	    }
	  else
	    {
	      if(options.has(RNAforesterOptions::LocalSimilarity))
		cout << "local optimal score: ";
	      else
		cout << "global optimal score: ";
	    }
	}

	if(options.has(RNAforesterOptions::SmallInLarge))
	  {
	    cout << ali.getSILOptimum() << endl;
	    optScore = ali.getSILOptimum();
	  }
	else
	  {
	    if(options.has(RNAforesterOptions::LocalSimilarity))
	      {
		cout << ali.getLocalOptimum() << endl;
		optScore = ali.getLocalOptimum();
	      }  
	    else
	      {
		cout << ali.getGlobalOptimum() << endl;
		optScore = ali.getGlobalOptimum();
		if(options.has(RNAforesterOptions::RelativeScore))
		  cout << ali.getGlobalOptimumRelative() << endl;
	      }
	  }


	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
	  if(options.has(RNAforesterOptions::SmallInLarge))
	    {
	      ali.getOptSILAlignment(ppfali,ybasepos);
	      cout << "starting at position: " << ybasepos << endl << endl; 
	    }
	  else
	    {	      
	      if(options.has(RNAforesterOptions::LocalSimilarity))
		{
		  ali.resetOptLocalAlignment(suboptPercent);
		  ali.getOptLocalAlignment(ppfali,xbasepos,ybasepos);		      		    
		  cout << "starting at positions: " << xbasepos << "," << ybasepos << endl << endl; 
		  xmlInfos.xbasepos = xbasepos;
		  xmlInfos.ybasepos = ybasepos;
		}
	      else		  
		ali.getOptGlobalAlignment(ppfali);		  
	    }

		// generate dot file
		makeDotFileAli(ppfali,options);

		// print alignment
		ppfali.getSequenceAlignments(seq1,seq2);				
		ppfali.getStructureAlignment(str1,true);
		ppfali.getStructureAlignment(str2,false);

		if(options.has(RNAforesterOptions::FastaOutput)){
		  cout << ppfali.getStructureNameX() << endl;
		  cout << seq1 << endl;
		  cout << str1 << endl;
		  cout << ppfali.getStructureNameY() << endl;
		  cout << seq2 << endl;
		  cout << str2 << endl;
		  cout << endl;
		}
		else {
		  RNAFuncs::printAli(ppfali.getStructureNameX(),ppfali.getStructureNameY(),seq1,seq2,str1,str2);
		}

		xlen=seq1.size();
		ylen=seq2.size();
		if(options.has(RNAforesterOptions::LocalSimilarity))
		{
			xsubopts.push_back(pair<Uint,Uint>(xbasepos,xlen));
			ysubopts.push_back(pair<Uint,Uint>(ybasepos,ylen));
		}

#ifdef HAVE_LIBG2  // This features require the g2 library		
		if(options.has(RNAforesterOptions::MakeSquigglePlot))
		  {
		    sprintf(s,"%d",count);
		    ppfali.squigglePlot(s,sqOptions);      
		  }
#endif	       

		// suboptimal alignments
		if(options.has(RNAforesterOptions::LocalSubopts))
		{
			while(ali.nextLocalSuboptimum())
			{
				count++;
				cout << "local optimal score: ";
				cout << ali.getLocalOptimum() << endl;
				ali.getOptLocalAlignment(ppfali,xbasepos,ybasepos);
				cout << "starting at positions: " << xbasepos << "," << ybasepos << endl << endl; 

				// print alignment
				ppfali.getSequenceAlignments(seq1,seq2);
				ppfali.getStructureAlignment(str1,true);
				ppfali.getStructureAlignment(str2,false);
				RNAFuncs::printAli(ppfali.getStructureNameX(),ppfali.getStructureNameY(),seq1,seq2,str1,str2);  
				xlen=seq1.size();
				ylen=seq2.size();
				xsubopts.push_back(pair<Uint,Uint>(xbasepos,xlen));
				ysubopts.push_back(pair<Uint,Uint>(ybasepos,ylen));	

#ifdef HAVE_LIBG2  // This features require the g2 library		
		                if(options.has(RNAforesterOptions::MakeSquigglePlot))		  
				{
				  sprintf(s,"%d",count);				
				  ppfali.squigglePlot(s,sqOptions);      	    	  
				}
#endif			       
			}
		}

	}

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library	
#ifdef HAVE_LIBXMLPLUSPLUS
#ifdef HAVE_LIBXML2
	// generate xml
	if(options.has(RNAforesterOptions::XmlOutputFile)){
	  options.get(RNAforesterOptions::XmlOutputFile,outputFile,string(""));
	} else {
	  outputFile = "result.xml";
	}
	if(options.has(RNAforesterOptions::GenerateXML))
	{
	  RNAFuncs::printPAliXML(ppfali.getStructureNameX(),ppfali.getStructureNameY(),seq1,seq2,str1,str2,optScore,options,xmlInfos,outputFile);
	}
#endif	
#endif
#endif

#ifdef HAVE_LIBG2  // This features require the g2 library
	// generate squiggle plot
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{
		f1->plot2d("x_str", xsubopts, sqOptions);
		f2->plot2d("y_str", ysubopts, sqOptions);  
	}
#endif	

	// clear input list
	deque<RNAForest*>::const_iterator it;
	for(it = inputListPW.begin(); it!=inputListPW.end(); it++)
		delete *it;
	inputListPW.clear();
	delete alg;
}


void editPairwise(list<RNAForestSZ*> &inputListSZ,Score &score,RNAforesterOptions &options)
{
  //  timeb t1,t2;
  IntScoreSZAlgebraType *alg;

//  if(options.has(RNAforesterOptions::CalculateDistance))
    alg=new IntDistSZAlgebra(score);
//  else
//    alg=new IntSimiSZAlgebra(score);
  
  RNAForestSZ *f1=inputListSZ.front();
  RNAForestSZ *f2=inputListSZ.back();

  //  ftime(&t1);
  Mapping<int,RNA_Alphabet> mapping(f1,f2,*alg);
  //  ftime(&t2);

//	f1->printParameters("F1");
//	f2->printParameters("F2");

  cout << "Global optimum: " << mapping.getGlobalOptimum() << endl;
  //  cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << endl;
}



void alignPairwiseSimple(deque<RNAForest*> &inputListPW,Score &score,RNAforesterOptions &options)
{
  //  timeb t1,t2;
  IntScore_AlgebraType *alg;

  //  if(options.has(RNAforesterOptions::CalculateDistance))
    alg=new IntDist_Algebra(score);
//  else
//    alg=new ScoreAlgebraSimple(score);

  RNAForest *f1=inputListPW.front();
  RNAForest *f2=inputListPW.back();

  //  ftime(&t1);
  Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg,false);
  //  ftime(&t2);

//	f1->printParameters("F1");
  //f2->printParameters("F2");

  cout << "Global optimum: " << ali.getGlobalOptimum() << endl;
  //  cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << endl;
}

#ifdef HAVE_LIBXML2
#ifdef HAVE_LIBXMLPLUSPLUS
extern "C" {
  bool validateXSD(string filename){
    xmlSchemaParserCtxtPtr ctxt;
    xmlSchemaValidCtxtPtr valCtxt;
    xmlSchemaPtr schema;
    int val;

    ctxt = xmlSchemaNewMemParserCtxt(XSD_STRING,sizeof(XSD_STRING));
        
    schema = xmlSchemaParse(ctxt);
    
    valCtxt = xmlSchemaNewValidCtxt(schema);

    val = xmlSchemaValidateFile(valCtxt,filename.c_str(),0);

    if(val==0){
      return true;
    } else {
      return false;
    }
  } 
}
#endif
#endif