// This file is part of BULL, a program for phylogenetic simulations
// most of the code was written by Mark T. Holder.

//	This program is for internal use by the lab of Dr. Tandy Warnow only.
//	Do not redistribute the code.  It is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  
//
//	Some of the code is from publically available source by Paul Lewis, Ziheng Yang, 
//	John Huelsenbeck, David Swofford , and others  (as noted in the code).
//	In fact the main structure of the program was created by modifying Paul Lewis' 
//	basiccmdline.cpp from his NCL
//
//	This code was used in Mark's dissertation, some changes were made in order to 
//	get it to compile on gcc.  It is possible that this porting introduced bugs (very
//	little debugging has been done on UNIX platforms).	I would suggest checking 
//	the simulator by generating data on trees with short branches, etc.


#include <sstream>
#include <algorithm>

#include "bull.hpp"
#include "string_extensions.hpp"
#include "set_reader.hpp"

using namespace std;
using namespace bull;

extern long gseed;

int GetHighestFileNum(std::string tag,int num=0);
void ReadNexusVector(std::vector<std::string> *v, NexusToken &token);

int GetHighestFileNum(std::string tag,int num/*=0*/) {
	ifstream testf;
	std::string tempstr;
	tempstr=tag;
	tempstr+=num;
	testf.open(tempstr.c_str());
	while (testf.good()) {
		testf.close();
		num++;
		tempstr=tag;
		tempstr+=num;
		testf.open(tempstr.c_str());
	}
	testf.close();
	return num - 1;
}



void Bull::HandleGetTrees(NexusToken& token ) {
	GetTreesOpts gto = ParseGetTrees(token);
	try {
		ExecGetTrees(gto);
	}
	catch (XBull & x) {
		throw XBull(x.msg, token);
	}
}

void Bull::HandleGetTrees(NexusToken& token ) {
	CodLikeStartOpts clso;
	ParseCodLikeStartValCommand(&clso, token);
	try {
		ExecCodLikeStartValCommand(clso);
	}
	catch (XBull & x) {
		throw XBull(x.msg, token);
	}
}

void Bull::HandleSimulate(NexusToken& token ) {
	SimulateOpts sto;
	ParseSimulate(&sto, token);
	try {
		ExecSimulate(sto);
	}
	catch (XBull & x) {
		throw XBull(x.msg, token);
	}
}

void Bull::RequireTokenAndAdvance(NexusToken & token, const char *expected, const char *cmd) {
	if ( !token.Equals(expected)) {
		errormsg = "Expecting ";
		errormsg.append(cmd);
		if (cmd != NULL) {
			errormsg += " in the ";
			errormsg.append(cmd);
			errormsg += "command,";
		}
		errormsg += " but found ";
		errormsg += token.GetToken();
		errormsg += " instead";
		throw XBull(errormsg, token);
	}

}
/// Throws an XBull excpetion if token is not ";"
void Bull::RequireSemicolon(NexusToken & token, const char *cmd) const {
	if ( !token.Equals(";"))
		ThrowNoSemicolon(token, cmd);
}

void Bull::ThrowNoSemicolon(NexusToken & token, const char *cmd) const {
	errormsg = "Expecting ';'";
	if (cmd != NULL) {
		errormsg += " to terminate the ";
		errormsg.append(cmd);
		errormsg += "command,";
	}
	errormsg += " but found ";
	errormsg += token.GetToken();
	errormsg += " instead";
	throw XBull(errormsg, token);
}
/**
 * @method HandleEndblock [void:protected]
 * @param token [NexusToken&] the token used to read from in
 * @throws XBull
 *
 * Called when the END or ENDBLOCK command needs to be parsed
 * from within the Bull block.	Basically just checks to make
 * sure the next token in  the data file is a semicolon.
 */
void Bull::HandleEndblock( NexusToken& token ) {
	token.GetNextToken();
	RequireSemicolon(token, "END");
}

/**
 * @method HandleExecute [void:protected]
 * @param token [NexusToken&] the token used to read from in
 * @throws XBull
 *
 * Handles everything after the EXecute keyword and the terminating
 * semicolon.  Flushes all blocks before executing file specified,
 * and no warning is given of this
 */
void Bull::HandleExecute( NexusToken& token ) {
	bool autopurge=false;
	token.GetNextToken();
	if (token.Equals("Purge")) {
		autopurge=true;
		token.GetNextToken();
	}
	if ( inf_open && ! autopurge) {
		errormsg = "Cannot issue execute command from within a Bull block";
		throw XBull(errormsg, token);
	}

	std::string fn = token.GetToken();

	token.GetNextToken();
	if (!autopurge && token.Equals("Purge")) {
		autopurge=true;
		token.GetNextToken();
	}
	RequireSemicolon(token, "EXECUTE");

	// Before going through with this, make sure we're not going to overwrite 
	// any stored blocks
	bool stuff_stored = !taxa->IsEmpty();
	stuff_stored = (stuff_stored || !trees->IsEmpty());
	if (stuff_stored) {
		if (!autopurge || !UserSaysOk("Ok to delete?", "Data has already been read and stored")) {
			message = "\nExecute command aborted.";
			PrintMessage();
			return;
		}
		ResetBlocks();
	}

	cout << "\nOpening " << fn << "..." << endl;
	
	ifstream inf(fn.c_str(), ios::binary);
	inf_open = true;
	NexusToken ftoken(inf);
	
	try {
		Execute( ftoken );
	}
	catch( XBull x ) {
		long line = x.line;
		long col = x.col; 
		istream::pos_type pos;
		if (x.validPos)
			pos = x.pos;
		if (errormsg.length() == 0)
			NexusError( x.msg, pos, line, col);
		else
			NexusError(errormsg, pos, line, col );
	  Reset();
	}
	
	if (inf_open)
		inf.close();
	inf_open = false;
	
	if ( purged && !taxa->IsEmpty() ) {
		purged=false;
		cout << "  TAXA block found" << endl;
		taxa->Report(logf);
		FinishTaxaBlock();
		if (!trees->IsEmpty())
			FinishTreesBlock();
	}

	if (!trees->IsEmpty()) {
		if (taxa->IsEmpty())
			cout << "TREES block found, but ignored because there is not an active taxa block"<<endl;
		else {
			cout << "  TREES block found" << endl;
			trees->Report(logf);
		}
	}
	cout << "Done with file "<<fn<<endl;
}

void Bull::SkipBlock(NexusToken &ftoken, const std::string & currBlock) const {
	SkippingBlock(currBlock);
	for (;;) {
		ftoken.GetNextToken();
		if ( ftoken.Equals("END") || ftoken.Equals("ENDBLOCK") ) {
			ftoken.GetNextToken();
			RequireSemicolon(ftoken, "END");
			break;
		}
		if ( ftoken.AtEOF() ) {
			errormsg = "Encountered end of file before END or ENDBLOCK in block ";
			errormsg += currBlock;
			throw XBull( errormsg, ftoken);
		}
	}
}

void Bull::HandleExecuteBullBlocksInFile( NexusToken & tok ) {
	//only allows bull block commands to be executed (not other blocks)
	// Get the file name to execute

	tok.GetNextToken();
	std::string fn = tok.GetToken();
	// get the semicolon terminating the EXECUTE command
	tok.GetNextToken();
	
	RequireSemicolon(tok, "EXECUTE");

	ifstream inf( fn.c_str(), ios::binary );
	NexusToken ftoken(inf);
	ftoken.GetNextToken();
	if ( !ftoken.Equals("#NEXUS") ) {
		errormsg = "Expecting #NEXUS to be the first token in the file, but found ";
		errormsg += ftoken.GetToken();
		errormsg += " instead";
		throw XBull(errormsg, ftoken);
	}
	
	ExecuteStarting();
	
	ftoken.SetLabileFlagBit( NexusToken::saveCommandComments );
	ftoken.GetNextToken();
	while (!ftoken.AtEOF()) {
		if ( ftoken.Equals("BEGIN") ) {
			ftoken.GetNextToken();
			if (ftoken.Equals( "BULL") ) {
				EnteringBlock( ftoken.GetToken() );
				ftoken.GetNextToken();
				RequireSemicolon(ftoken, "BEGIN BULL");
				for (;;) {
					ftoken.GetNextToken();
					if ( ftoken.Abbreviation("ENdblock") ) {
						HandleEndblock( ftoken );
						break;
					}
					else if ( ftoken.Abbreviation("Quit") ) {
						message = "\nIgnoring Quit command from non-primary file\n";
						PrintMessage();
						break;
					}
					TryAdvancedCommands(ftoken);
				}
			}
			else 
				SkipBlock(ftoken, ftoken.GetToken());
		}
		ftoken.SetLabileFlagBit( NexusToken::saveCommandComments );
		ftoken.GetNextToken();
	}
	
	ExecuteStopping();
}


GetTreesOpts Bull::ParseGetTrees(NexusToken& token ) {
	assert(gto != NULL);
	if (!taxa) {
		errormsg="You can't get trees without an active set of taxa";
		throw XBull( errormsg, token);
	}
	bool readingMostRecent=false;
	long fromTree=-1;
	long toTree=-1
	unsigned mode=7;
	
	token.GetNextToken();
	while (token.GetToken()!=";") {
		if (token.Abbreviation("STOREBrlens")) {
			token.GetNextToken();
			gto->storeBrLensFromFile=true;
			if (token.GetToken() == "=") {
				token.GetNextToken();
				if (token.Abbreviation("No"))
					gto->storeBrLensFromFile=false;
				else if (!token.Abbreviation("Yes")) {
					errormsg="Expecting YES or NO after CurrentBranchLengths = option to LSCORE command";
					throw XBull(errormsg, token);
				}
			}
			token.GetNextToken();
		}
		if (token.Abbreviation("FRom")) {
			fromTree=GetLongInteger(token);
			if (fromTree < 0) {
				errormsg="Expecting positive integer after FROM option of GETTREES";
				throw XBull(errormsg, token);
			}	
		}
		if (token.Abbreviation("To")) {
			toTree=GetLongInteger(token);
			if (toTree < 0) {
				errormsg="Expecting positive integer after TO option of GETTREES";
				throw XBull(errormsg, token);
			}	
		}
		if (token.Abbreviation("REplace"))
			mode=3;
		if (token.Abbreviation("Mode")) {
			mode=GetLongInteger(token);
			if (mode!=3 && mode!=7) {
				errormsg="Right now bull only get trees with mode 3 or mode 7";
				throw XBull(errormsg, token);
			}	
		}
		if (token.Abbreviation("FIle"))
			filename=GetFileName(token);
		if (token.Abbreviation("PREfix")) {
			token.GetNextToken();
			if (token.GetToken() == "=")	
				token.GetNextToken();
			filePref=token.GetToken();
		}
		if (token.Abbreviation("MOStrecent"))
			readingMostRecent=true;
	}
	
	if (readingMostRecent) {
		int highNum = GetHighestFileNum(filePref);
		if (highNum < 0) {
			errormsg = "Couldn't open file "; errormsg += filePref; errormsg += "0";
			throw XBull(errormsg, token);
		}  
		filename=filePref;
		filename+=highNum;
	}
	
	gto->filename = filename;
	gto->fromTree = (fromTree < 2 ? 0 : fromTree - 1);
	gto->toTree = (toTree < 1 ? -1 : toTree - 1);
	gto->mode = mod;
	return gto
}

void Bull::ExecGetTrees(const GetTreesOpts & gto) {
	// Read the file that should contain the trees
	ifstream treefile(gto.filename.c_str());
	if (!treefile.good()) {
		string errormsg = "Couldn't open the file ";
		errormsg += gto.filename;
		throw XBull(errormsg);
	}
	
	MyNexus tempNexus(this);
	TreesBlock tempTree(*taxa);
	tempNexus.Add(&tempTree);
	NexusToken ftoken(treefile);

	try {
		tempNexus.Execute( ftoken );
	}
	catch( XBull x ) {
		const std::string & m = (errormsg.empty() ? x.msg : errormsg); 
		tempNexus.NexusError(m, x.pos, x.line, x.col );
	}
	if (treefile.good())
		treefile.close();

	if (tempTree.IsEmpty()) {
		errormsg="No Trees Block found in "; errormsg+=filename;
		throw XBull(errormsg, token);
	}
	
	const unsigned ntreesInFile = tempTree.GetNumTrees();
	if (ntreesInFile < 1) {
		errormsg="No Trees found in ";
		errormsg+=filename;
		throw XBull(errormsg);
	}
	
	const unsigned toTree = (unsigned) (gto->toTree < 0 ? ntreesInFile - 1 : gto->toTree); 
	const unsigned fromTree = gto->fromTree;
	if (fromTree > toTree) {
		errormsg="Cant get Trees from ";
		errormsg+=fromTree; 
		errormsg+=" to "; 
		errormsg+=toTree;
		throw XBull( errormsg, token);
	}
		

	// Store the specified trees
	// these trees are deleted in the for loop if there is an error, otherwise
	//	they are put under the control of the kernel
	vector<Tree *>  treesToBeAdded; 
	
	for (unsigned i = fromTree; i < toTree; i++) {
		const std::string newick = tempTree.GetTranslatedTreeDescription(i); 
		Tree *temptree = new Tree(newick, storeBrLensFromFile);
		if (!temptree->IsGood()) {
			for (int j=0;j < i-fromTree;j++)
				delete treesToBeAdded[j];
			errormsg="Problem Reading Tree Description of  "; 
			errormsg+=tempTree.GetTreeName(i);
			delete temptree;
			throw XBull( errormsg, token);
		}
		std::string s=tempTree.GetTreeName(i);
		ToUpper(s);
		temptree->SetName(s);
		treesToBeAdded.push_back(temptree);
		}

	ioObject.message.clear();
	ioObject.message += (int)(toTree-fromTree); 
	ioObject.message += " trees read from ";
	ioObject.message += filename;
	ioObject.printMessage();

	kernel->appendTrees(gto.mode, treesToBeAdded);
}		


/**
 * @method HandleLog [void:protected]
 * @param token [NexusToken&] the token used to read from in
 * @throws XBull
 *
 * Called when the LOG command needs to be parsed
 * from within the Bull block.
 */
void Bull::HandleLog( NexusToken& token )
{
	bool starting = false;
	bool stopping = false;
	bool appending = false;
	bool replacing = false;
	bool name_provided = false;
	std::string logfname;
	
	// Retrieve all tokens for this command, stopping only in the event
	// of a semicolon or an unrecognized keyword
	//
	for (;;)
	{
		token.GetNextToken();

		if ( token.Equals(";") ) {
			break;
		}
		else if ( token.Abbreviation("STOp") ) {
			stopping = true;
		}
		else if ( token.Abbreviation("STArt") ) {
			starting = true;
		}
		else if ( token.Abbreviation("Replace") ) {
			replacing = true;
		}
		else if ( token.Abbreviation("Append") ) {
			appending = true;
		}
		else if ( token.Abbreviation("File") ) {
			logfname = GetFileName(token);
			name_provided = true;
		}
		else {
			errormsg = "Unexpected keyword (";
		  errormsg += token.GetToken();
		  errormsg += ") encountered reading LOG command";
			throw XBull( errormsg, token);
		}
   }
   
   // Check for incompatible combinations of keywords
   //
   if ( stopping && ( starting || appending || replacing || name_provided ) ) {
		errormsg = "Cannot specify STOP with any of the following START, APPEND, REPLACE, FILE";
		throw XBull( errormsg, token);
   }
   
   if ( appending && replacing ) {
		errormsg = "Cannot specify APPEND and REPLACE at the same time";
		throw XBull( errormsg, token);
   }

   if ( logf_open && ( starting || name_provided || appending || replacing ) ) {
		errormsg = "Cannot start log file since log file is already open";
		throw XBull( errormsg, token);
   }

   // Is user closing an open log file?
   //
   if ( stopping ) {
		logf.close();
		logf_open = false;
		logEachStep=false;
		
		//message = "\nLog file closed";
		//PrintMessage();
		cout<<"\nLog file closed"<<endl;
		return;
   }
   
	// If this far, must be attempting to open a log file
	//
	if ( !name_provided ) {
		errormsg = "Must provide a file name when opening a log file\n";
		errormsg += "e.g., log file=doofus.txt start replace;";
		throw XBull( errormsg, token);
	}	
	
	if ( appending ) {
		logf_open = true;
		logf.open( logfname.c_str(), ios::out | ios::app );

		message = "\nAppending to log file ";
		message += logfname;
		PrintMessage();
	}
	else if ( replacing ) {
		logf_open = true;
		logf.open( logfname.c_str() );

		message = "\nReplacing log file ";
		message += logfname;
		PrintMessage();
	}
	else  {
		bool exists = FileExists( logfname.c_str() );
		bool userok = true;
		if ( exists && !UserSaysOk( "Ok to replace?", "Log file specified already exists" ) )
			userok = false;
		if ( userok ) {
			logf_open = true;
			logf.open( logfname.c_str() );
			logf<<setprecision(8);
		}

		if ( exists && userok ) {
			//message = "\nReplacing log file ";
			//message += logfname;
			cout<<"\nReplacing log file "<<logfname<<endl;
		}
		else if ( userok ) {
			//message = "\nLog file ";
			//message += logfname;
			//message += " opened";
			cout<<"\nLog file "<<logfname<<" opened"<<endl;
		}
		else {
			message = "\nLog command aborted";
		}
		PrintMessage();
	}
}

/**
 * @method PreprocessNextCommand [void:public]
 *
 * Accepts a string in the form of a Bull block containing one command
 * and processes it just like a real Bull block in a NEXUS data file.
 */
void Bull::HandleNextCommand()
{
	istringstream cmdin( next_command );

	NexusToken token(cmdin);
	try {
		Read( token );
	}
	catch( XBull x ) 
	{
		if (errormsg.length() == 0)
				NexusError( x.msg, x.pos, x.line, x.col );
		else	NexusError( errormsg, x.pos, x.line, x.col );
	  Reset();
	}
}

void ReadNexusVector(std::vector<std::string> *v, NexusToken &token) {
	if (token.Equals("(")) {
		token.GetNextToken();
		while (!token.Equals(")")) {
			v->push_back(token.GetToken());
			token.GetNextToken();
		}
	}
	else 
		v->push_back(token.GetToken());
}


EncodingType Bull::InterpretOutputType(NexusToken &token) {
	if (token.Abbreviation("PROtein") || token.Abbreviation("AMino"))
		return EncodingType(AminoAcid);
	if (token.Abbreviation("Dna"))
		return EncodingType(DNANoGap);
	errormsg="The output type ";
	errormsg+=token.GetToken(); 
	errormsg+=" is invalid";
	throw XBull( errormsg, token);
}


void Bull::HandleSimulate(NexusToken& token ) {	
	if (treelist.empty()) {
		errormsg="You can't simulate without a tree in memory";
		throw XBull( errormsg, token);
	}

	//Parse command into SimulateOpts struct
	SimulateOpts sto = ParseSimulate(token);
	ExecSimulateCommand(sto);
}

SimulateOpts Bull::ParseSimulate(NexusToken& token ) const {	
	token.GetNextToken();
	int nOutputs = 0;
	int nSimChars = 0;
	bool automatic = false;
	SimulateOpts simOpts;
	
	for (; token.GetToken() != ";"; token.GetNextToken()) { 
		bool needAnotherToken = true;
		if (token.Abbreviation("NReps"))
			simOpts.nReps = GetLongInteger(token);
		else if (token.Abbreviation("COLLapsefile"))
			simOpts.collapsedTreeFilename = GetFileName(token);
		else if (token.Abbreviation("NChars"))
			nSimChars = GetLongInteger(token);
		else if (token.Abbreviation("Concat"))
			simOpts.concatenations = GetLongInteger(token);
		else if (token.Abbreviation("NOUTput")) {
			nOutputs=GetLongInteger(token);
			if (nOutputs < 1) {
				errormsg="nOutputs must be > 0";
				throw XBull( errormsg, token);
			}
		}
		else if (token.Abbreviation("Paupblockfile")) {
			token.GetNextToken();
			RequireTokenAndAdvance(token, "=", "Simulate");
			ReadNexusVector(&(simOpts.paupBlockFilename), token);
		}
		else if (token.Abbreviation("OVERWrite"))
			simOpts.overwrite = true;
		else if (token.Abbreviation("File")) {
			token.GetNextToken();
			RequireTokenAndAdvance(token, "=", "Simulate");
			ReadNexusVector(&(simOpts.outputFilenames), token);
		}
		else if (token.Abbreviation("TAg"))
			tagname = GetFileName(token);
		else if (token.Abbreviation("AUTOmatic"))
			automatic = true;
		else if (token.Abbreviation("OUTputtypes")) {
			token.GetNextToken();
			RequireTokenAndAdvance(token, "=", "Simulate");
			if (token.Equals("(")) {
				token.GetNextToken();
				while (!token.Equals(")")) {
					simOpts.outTypes.push_back(InterpretOutputType(token));
					token.GetNextToken();
				}
			}
			else
				simOpts.outTypes.push_back(InterpretOutputType(token));
		}
		else {
			Tree * temptree;
			if(token.IsInteger()){
				long tn = token.GetLongEquivalent();
				if (tn > treeList.size() || tn < 1) {
					errormsg = "Tree ";
					errormsg += token.GetToken();
					errormsg += " unknown.  simulation aborted.";
					throw XBull( errormsg, token);
				}
				temptre = treeList[tn-1];
			}
			else {
				try {
					temptre = FindTreeFromName(token.GetToken());
				}
				catch (NoSuchTree) {
					errormsg="Tree ";
					errormsg+=token.GetToken();
					errormsg+=" unknown.";
					throw XBull( errormsg, token);
				}
			}
			simOpts.treeName.push_back(token.GetToken());//keeps track of what the user called the tree
			simOpts.treeAlias.push_back(temptre);
		}	
	}


	if (simOpts.outTypes.size() == 0){
		//dna is the default output type
		simOpts.outTypes.push_back(EncodingType(DNANoGap));
	}
	if ((nOutputs > 0) && simOpts.outTypes.size() != nOutputs) {
		errormsg = "The number of OutputTypes in a Simulate command must match the number of specified in the NOutput option";
		throw XBull( errormsg, token);
	}
	if (automatic) {
		if (simOpts.tagname.empty()) {
			errormsg = "The Tag option of the Simulate command must be used if the Automatic naming mode is requested.";
			throw XBull( errormsg, token);
		}
		if (!simOpts.outputFilenames.empty()) {
			errormsg = "The FILE option of the Simulate command cannot be used if the Automatic naming mode is requested (use the Tag instead of file).";
			throw XBull( errormsg, token);
		}
		simOpts.outputFilenames.clear();
		std::vector<EncodingType>::const_iterater ot = simOpts.outTypes.begin();
		for (; ot != simOpts.outTypes.end(); ot++) {
			string ofname = EncodingTypeToString(*ot);
			ofname.append(simOpts.tagname);
			simOpts.outputFilenames.append(ofname);
		}
	}
	else if (simOpts.outTypes.size() != simOpts.outputFilenames.size()) {
		errormsg = "The number of OutputTypes in a Simulate command must match the number of Files specified";
		throw XBull( errormsg, token);
	}
	if ((!simOpts.paupBlockFilename.empty()) && simOpts.paupBlockFilename.size() != simOpts.outTypes.size()) {
		errormsg = "The number of PaupBlockFile in a Simulate command must match the number of OutputTypes specified";
		throw XBull( errormsg, token);
	}
	

	//TEMPORARY lots of the interaction with SSettings assumes that you are using SSRFCodonSubModel
	if (modVec.empty()) {
		errormsg="You can't simulate a codon model without first declaring the models using codLikeStartVal";
		throw XBull( errormsg, token);
	}
	
	if ( simOpts.concatenations == 0) {
		if (nSimChars == 0) {
			errormsg="Either the number of characters or the number of concatenations must be specified";
			throw XBull( errormsg, token);
		}
		const unsigned naa = this->getNumCodonsInModel();
		const unsigned ndna = 3*naa;
		if (nSimChars % ndna) {
			errormsg="The number of characters must be a multiple of the number of amino acids";
			throw XBull( errormsg, token);
		}
		simOpts.concatenations = ndna/nSimChars;
	}

	return simOpts;
}
void Bull::HandleCodLikeStartVal(NexusToken& token )
{	
	int currbrlen=-1,nst=-1,naa=-1;
	double *sharedModParam,kapp=-1.0;
	workingOnBranches=false;
	sharedModParam=new double [10];
	for (int i=0;i < 10;i++)
		sharedModParam[i]=-1.0;
	double **aafreq,*mults;
	double treesc=1.0;
	int code;
//#ifdef ALLOWMULTIHITS
	double mh=0.0;
//#endif	
	aafreq=NULL;
	mults=NULL;
	HandleCodLikeStartVal(token,currbrlen,nst,naa,sharedModParam,kapp,aafreq,mults,treesc,code,mh);
	delete []sharedModParam;
	delete []mults;
	free_RectMats(aafreq);
}

void Bull::HandleCodLikeStartVal(NexusToken& token,int& currbrlen,int& nst,int& naa,double *&sharedModParam,
									double& kapp,double **&aafreq,double *&mults,double& treesc,int &code,double &mh)
{	delete[] previousLikelihoods;
	previousLikelihoods=NULL;
	nParamDirectionsToSkip=0;
	nParamGroupsToSkip=0;
	nBranchesToSkip=0;
	
	paramImprovThisRound=0.0;
	ParseCodLikeStartValCommand(token,currbrlen,nst,naa,sharedModParam,kapp,aafreq,mults,treesc,code,mh);
	if (previousLikelihoods)
		CenterOverflowMultipliers(mults,previousLikelihoods,naa,1610);

	//TEMPORARY should be modified to allow parameter/options to be persistent
	if (nst == 1) 
		{sharedModParam[3]=sharedModParam[4]=sharedModParam[5]=sharedModParam[6]=sharedModParam[7]=sharedModParam[8]=1.0;}
	else if (nst == 2)	
		{if(kapp<-1)	
			{if(sharedModParam[3]<0.0) {
				delete []mults; 
				delete [] previousLikelihoods; 
				previousLikelihoods=NULL;
				free_RectMats(aafreq);
				throw XBull("Kappa must be specified  option to CodLikeStartVal command");
			}
			else if (fabs(sharedModParam[3]-sharedModParam[5]) > 0.00001 || fabs(sharedModParam[3]-sharedModParam[6]) > 0.00001|| fabs(sharedModParam[3]-sharedModParam[8]) > 0.00001 || fabs(sharedModParam[4]-sharedModParam[7]) > 0.00001) {
				delete []mults;
				delete [] previousLikelihoods;
				previousLikelihoods=NULL;
				free_RectMats(aafreq);
				throw XBull("The specified rmatrix doenst conform to K2P in CodLikeStartVal command");
			}
			sharedModParam[4]/=sharedModParam[3];
			sharedModParam[7]/=sharedModParam[3];
			sharedModParam[3]=sharedModParam[4]=sharedModParam[5]=sharedModParam[8]=1.0;
			}
		else
			{sharedModParam[3]=sharedModParam[5]=sharedModParam[6]=sharedModParam[8]=1.0;
			sharedModParam[4]=sharedModParam[7]=kapp;
			}
		}
	else if (nst == 6) {
		if(sharedModParam[3]<0.0) {
			delete []mults;
			delete [] previousLikelihoods;
			previousLikelihoods=NULL;
			free_RectMats(aafreq);
			throw XBull("RMatrix must be specified	 option to CodLikeStartVal command");
		}
	}
	else if (kapp < 0.0) {
		if(sharedModParam[3]<0.0) {
			sharedModParam[3]=sharedModParam[4]=sharedModParam[5]=sharedModParam[6]=sharedModParam[7]=sharedModParam[8]=1.0;
		}
	}
	else {
		sharedModParam[3]=sharedModParam[5]=sharedModParam[6]=sharedModParam[8]=1.0;
		sharedModParam[4]=sharedModParam[7]=kapp;
	}
	if (sharedModParam[1]<0.0) {
		sharedModParam[0]=sharedModParam[1]=sharedModParam[2]=0.25;
	}
	if (naa < 1) {
		delete []mults; 
		delete [] previousLikelihoods;
		previousLikelihoods=NULL;
		free_RectMats(aafreq);
		throw XBull("You must specify the number of amino acid freqs in CodLikeStartVal command");
	}
#ifdef	ALLOWMULTIHITS
	sharedModParam[9]=mh;
#endif

Model **arrOfMods;
	arrOfMods=new Model*[naa];
	modArr=new SSRFCodonSubMod*[naa];
	sizeOfModArr=naa;
	CreateSetOfSSRFCodonSubMods(naa,sharedModParam,aafreq,modArr,mults,treesc,code);
	for (int k=0;k < naa;k++)
		arrOfMods[k]=(Model *)modArr[k];
	delete []arrOfMods;
	
}


void CenterOverflowMultipliers(double *mults,double *previousLikelihoods,int naa,int ntax)
{	for (int i=0;i < naa;i++)
		{double y=exp((*previousLikelihoods++)/((double) ntax));
		*mults++=y;
		cout<<y<<endl;
		}
}

void Bull::ParseCodLikeStartValCommand(NexusToken& token,int& currbrlen,int& nst,int& naa,
	double *&sharedModParam,double& kapp,double **&aafreq,double *&mults,double& treesc,int &code,double &mh)
{	code=GenCode(Mito);//default is Mito
	token.GetNextToken();
	while (token.GetToken()!=";")
		{bool needsAnother=true;
		if (token.Abbreviation("CURREntbranchlengths") || token.Abbreviation("CURRBranchlengths"))
			{token.GetNextToken();
			if (token.GetToken() == "=")
				{token.GetNextToken();
				if (token.Abbreviation("No"))
					currbrlen=0;
				else if (token.Abbreviation("Yes"))
							currbrlen=1;
						else	{while (token.GetToken()!=";")	token.GetNextToken();
								delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
								throw XBull("Expecting YES or NO after CurrentBranchLengths = option to CodLikeStartVal command");
								}
				token.GetNextToken();
				}
			else	{currbrlen=1;needsAnother=false;}
			}
		else if (token.Abbreviation("LASTParamimprovement"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			lastParamImprov=atof(token.GetToken().c_str());
			if (lastParamImprov < 0.0)
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); 
				throw XBull("lastParamImprov  must be  >= 0.0 CodLikeStartVal command");
				}
			}
		else if (token.Abbreviation("PARAMIMPROVEThisround"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			paramImprovThisRound=atof(token.GetToken().c_str());
			if (paramImprovThisRound < 0.0)
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); 
				throw XBull("paramImprovThisRound	must be  >= 0.0 CodLikeStartVal command");
				}
			}
		else if (token.Abbreviation("LASTBranchimprovement"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			lastBranchImprov=atof(token.GetToken().c_str());
			if (lastBranchImprov < 0.0)
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); 
				throw XBull("lastBranchImprov	must be  >= 0.0 CodLikeStartVal command");
				}
			}
		else if (token.Abbreviation("GRoupskip"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			nParamGroupsToSkip=(int) token.GetLongEquivalent();
			if (nParamGroupsToSkip < 0)
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
					throw XBull("Expecting positive number before Groupskip option to CodLikeStartVal command");
					}
			}
		else if (token.Abbreviation("PARAMskip"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			nParamDirectionsToSkip=(int) token.GetLongEquivalent();
			if (nParamDirectionsToSkip < 0)
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
					throw XBull("Expecting positive number before Paramskip option to CodLikeStartVal command");
					}
			}
		else if (token.Abbreviation("BRanchskip"))
			{
			throw MTHException("branchskip not supported");
			}
		else if (token.Abbreviation("GEneticCode"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (token.Abbreviation("Mitochondrial"))
				code=GenCode(Mito);
			else
				if (token.Abbreviation("Nuclear"))
					code=GenCode(Nuclear);
				else 
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
					throw XBull("Expecting either Mito or Nuclear code");
					}
				
			}
		else if (token.Equals("nst"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (token.Equals("1") || token.Equals("2") || token.Equals("6"))
				nst=token.GetLongEquivalent();
			else	{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
					throw XBull("Expecting 1 , 2 or 6 after nst option to CodLikeStartVal command");
					}
			}
		else if (token.Abbreviation("BASEfreq"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (!token.Equals("(")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); 
			throw XBull("Expecting ( after basefreq option to CodLikeStartVal command");}
			for (int i=0;i < 3;i++)
				{token.GetNextToken();
				sharedModParam[i]=atof(token.GetToken().c_str());
				if (sharedModParam[i] <= 0.0) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq);
				throw XBull("basefreqs must be > 0.0 CodLikeStartVal command");}
				}
			if (sharedModParam[1]+sharedModParam[2]+sharedModParam[0] >= 1.0)
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Sum of A, C, and G must be <1.0 CodLikeStartVal command");}
			token.GetNextToken();
			if (!token.Equals(")")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) after basefreq option to CodLikeStartVal command");}
			}
		else if (token.Abbreviation("RMATrix"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (!token.Equals("(")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ( after rmatrix option to CodLikeStartVal command");}
			for (int i=0;i < 6;i++)
				{token.GetNextToken();
				sharedModParam[3+i]=atof(token.GetToken().c_str());
				if (sharedModParam[3+i] <= 0.0)	
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); 
					throw XBull("rmatrix  must be >0.0 CodLikeStartVal command");}
				}
			token.GetNextToken();
			if (!token.Equals(")")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) after rmatrix	 option to CodLikeStartVal command");}
			}
		else if (token.Abbreviation("KAPpa"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			kapp=atof(token.GetToken().c_str());
			if (kapp < 0.0) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("kappa	 must be >0.0 CodLikeStartVal command");}
			}
#ifdef	ALLOWMULTIHITS
		else if (token.Abbreviation("DOublehit"))//not a great name, but allows multiplier abbreviation to work
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			mh=atof(token.GetToken().c_str());
			if (mh < 0.0)	{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("MultipleHitProb  must be >= 0.0 CodLikeStartVal command");}
			}
#endif
		else if (token.Abbreviation("AAFreq"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (!token.Equals("(")) 
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ( after aafreq option to CodLikeStartVal command");}
			token.GetNextToken();
			if (!token.IsInteger()) 
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting number of codons after aafreq=( option to CodLikeStartVal command");}
			naa=token.GetLongEquivalent();
			if (naa < 1)	
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("number of codons must be >0 in  CodLikeStartVal command");}
			aafreq=new_RectMats(naa,20);
			for (int ii=0;ii < naa;ii++)
				{token.GetNextToken();
				if (!token.Equals("(")) 
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ( for next site in aafreq option to CodLikeStartVal command");}
				for (int j=0;j < 20;j++)
					{token.GetNextToken();
					aafreq[ii][j]=atof(token.GetToken().c_str());
					if (aafreq[ii][j]<0.0)	
						{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("each amino acid freq	must be >0.0 CodLikeStartVal command");}
					}
				token.GetNextToken();
				if (!token.Equals(")")) 
					{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) for next site in aafreq option to CodLikeStartVal command");}
				}
			token.GetNextToken();
			if (!token.Equals(")")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) after aafreq	option to CodLikeStartVal command");}
			}
		else if (token.Abbreviation("MUlt"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (!token.Equals("(")) 
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ( after mult option to CodLikeStartVal command");}
			if (naa < 1)	
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("number of codons must be specified before the mult option is used in CodLikeStartVal command");}
			mults=new double[naa];
			for (int ii=0;ii < naa;ii++)
				{token.GetNextToken();
				mults[ii]=atof(token.GetToken().c_str());
				if (mults[ii]<0.0)	
						{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("multiplier  must be >0.0 CodLikeStartVal command");}
				}
			token.GetNextToken();
			if (!token.Equals(")")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) after mult  option to CodLikeStartVal command");}
			}
		else if (token.Abbreviation("PRevlike"))
			{if(naa < 1)	
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("number of codons must be specified before the mult option is used in CodLikeStartVal command");}
			if (previousLikelihoods) delete[] previousLikelihoods;
			previousLikelihoods=new double [naa];
			token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			if (!token.Equals("(")) 
				{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ( after prevlike option to CodLikeStartVal command");}
			if (workingOnBranches && previousLikelihoods)
				{for (int ii=0;ii < naa;ii++)
					token.GetNextToken();
				}
			else
				{previousLikelihoods=new double[naa];
				for (int ii=0;ii < naa;ii++)
					{token.GetNextToken();
					previousLikelihoods[ii]=atof(token.GetToken().c_str());		
					}
				}
			token.GetNextToken();
			if (!token.Equals(")")) {delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("Expecting ) after prevLike  option to CodLikeStartVal command");}
			}
		else if (token.Abbreviation("TReescale"))
			{token.GetNextToken();
			if (token.Equals("="))	token.GetNextToken();
			treesc=atof(token.GetToken().c_str());
			if (treesc < 0.0)	{delete []mults; delete [] previousLikelihoods; previousLikelihoods=NULL;free_RectMats(aafreq); throw XBull("The Tree scaling	must be >0.0 CodLikeStartVal command");}
			}
		if (needsAnother)	token.GetNextToken();
		}
}