caferror.py

#!/usr/bin/python
#############################################################################
# Genome Assembly and Annotation Error Estimation Algorithm Utilizing CAFEv3.0
# Written by: Gregg Thomas, Summer 2012
# Hahn Lab, Indiana University
# Contact: grthomas@indiana.edu
#
# This version of the script first minimizes scores for varying error models across all species 
# and then (if specified) incrementally adds or subtracts 10% of the global error to each species 
# to minimize scores and estimate error for each species individually.
#
# May 2016: Added an error message that displays if a lambda command with more than one search
# parameter is input. We don't recommend doing this. We do recommend estimating error with a
# 1 lambdda model and then subsequently using that error on a multi-lambda run of CAFE.
# Streamlined the output code with the printWrite function.
# Removed getAllSpec function and replaced it with regex.
#############################################################################

import sys, argparse, os, random, datetime, time, re

############################################
#Function Definitions
############################################
def errorOut(errnum, errmsg):
# Formatting for error messages.
	fullmsg = "|**Error " + str(errnum) + ": " + errmsg + " |";
	border = " " + "-" * (len(fullmsg)-2);
	print "\n" + border + "\n" + fullmsg + "\n" + border + "\n";

############################################
def optParse(errorflag):
# This function handles the command line options.

	parser = argparse.ArgumentParser();

	parser.add_argument("-i", dest="input_file", help="A CAFE shell script with the full CAFE path in the shebang line, the load, tree, and lambda commands. These lines will be read and incorporated into the caferror shell script.");
	parser.add_argument("-e", dest="user_err_start", help="The starting point for the grid search. Should be between 0 and 1. Default: 0.4", type=float, default=0.4);
	parser.add_argument("-d", dest="user_tmp_dir", help="A directory in which all caferror files will be stored. If none is specified, it will default to caferror_X, with X being some integer one higher than the last directory.", default="");
	parser.add_argument("-f", dest="first_run", help = "Boolean option to perform a pre-error model run (1) or not (0). Default: 0", type=int, default=1);
	parser.add_argument("-c", dest="curve_option", help="Boolean option. caferror can either perform the grid search (0) or search a pre-specified space (1). Default: 0", type=int, default=0);
	parser.add_argument("-t", dest="error_tries", help="A list of error values to search over. Note: -c MUST be set to 1 to use these values. Enter as a comma delimited string, ie -t 0.1,0.2,0.3", default="")
	parser.add_argument("-l", dest="user_log_file", help="Specify the name for caferror's log file here. Default: caferrorLog.txt", default="caferrorLog.txt");
	parser.add_argument("-o", dest="output_file", help="Output file which stores only the error model and score for each run. Default: caferror_default_output.txt", default="caferror_default_output.txt");
	parser.add_argument("-s", dest="ind_min", help="Boolean option to specify whether to perform only the global error search (0) or continue with individual species minimizations (1). Default: 0", type=int, default=0);
	parser.add_argument("-v", dest="verbose", help="Boolean option to have detailed information for each CAFE run printed to the screen (1) or not (0). Default: 1", type=int, default=1);
	parser.add_argument("-m", dest="run_mode", help=argparse.SUPPRESS, type=int, default=0);

	args = parser.parse_args();

	if errorflag == 0:

		if args.input_file == None:
			errorOut(1, "-i must be defined");
			optParse(1);

		if not all(op in [0,1] for op  in [args.first_run, args.curve_option, args.ind_min, args.verbose]):
			errorOut(2, "-f, -c, -s, and -v must all take values of either 0 or 1");
			optParse(1);

		if args.user_err_start > 1 or args.user_err_start < 0:
			errorOut(3, "-e must take values of between 0 and 1");
			optParse(1);

		if args.curve_option == 0 and args.error_tries != "":
			errorOut(4, "With -t specified -c must also be set to 1");
			optParse(1);

		if args.user_tmp_dir != "" and args.user_tmp_dir[len(args.user_tmp_dir) - 1] != "/":
			args.user_tmp_dir = args.user_tmp_dir + "/";

		if args.curve_option == 1 and args.error_tries == "":
			args.error_tries = [0.0, 0.001, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95];
		elif args.curve_option == 1 and args.error_tries != "":
			args.error_tries = args.error_tries.split(",");
		
		### -m: Run mode for individual species minimzations (for debugging purposes only... not supported).
		###		0 - Default: No shuffle, background constant.
		###		1 - No shuffle, background updated.
		###		2 - Shuffle, background constant.
		###		3 - Shuffle, background updated.

		return args.input_file, args.output_file, args.user_err_start, args.user_tmp_dir, args.curve_option, args.error_tries, args.user_log_file, args.ind_min, args.run_mode, args.verbose, args.first_run;

	elif errorflag == 1:
		parser.print_help();
		sys.exit();

############################################
def printWrite(o_name, v, o_line1, o_line2="", pad=0):
# Function to print a string AND write it to the file.
	if o_line2 == "":
		outline = o_line1;
	else:
		outline = o_line1 + " "*(pad-len(o_line1)) + o_line2;
	if v == 1 or v == -2:
		print outline;
	f = open(o_name, "a");
	f.write(outline + "\n");
	f.close();

############################################
def cafeRun(families, newick, lamb, path, tempDir, cafLog, specDict, spectomin, logrm, fflag, nr, vb):
#This is the function which creates an error file based on the current error (by calling maxfamsize), writes the appropriate caferror.sh script
#and executes the current caferror.sh script.

	printWrite(cafLog, 1, "--------------------------");

	outline = "CAFE run number " + str(nr);
	if spectomin[0] == "final":
		outline = outline + "(FINAL RUN)";
	printWrite(cafLog, 1, outline);

	if spectomin[0] == "all":
		printWrite(cafLog, 1, str(specDict[errSpec[0]]) + " Error Model Run for all");
		logFile = tempDir + "cafe_" + str(specDict[errSpec[0]]) + "_all_log.txt";
	elif spectomin[0] == "none":
		printWrite(cafLog, 1, "Pre run (no error models)");
		logFile = tempDir + "cafe_pre_log.txt";
	elif spectomin[0] == "final":
		printWrite(cafLog, 1, "With minimized error models");
		logFile = tempDir + "cafe_final_log.txt";
		repFile = tempDir + "cafe_final_report";
	else:
		printWrite(cafLog, 1, str(specDict[spectomin[0]]) + " Error Model Run for " +  spectomin[0]);
		logFile = tempDir + "cafe_" + str(specDict[spectomin[0]]) + "_" + spectomin[0] + "_log.txt";
	printWrite(cafLog, 0, "-----");
	# The above logic statements name the log file for the current CAFE run. The format for the log files is: cafe_[errormodel]_[species]_log.txt.
	# Special cases for species include 'all' in which the current errormodel is applied to all during global error estimation species and 'final' 
	# which is the last CAFE run with the minimized errormodels applied to each individual species.

	if logrm == 1:
		os.system("rm " + logFile);
	#In the event that CAFE initializes incorrectly, it must be re-run with the same parameters. When that happens, the old log file is removed here
	#and will be replaced with one of the same name. The contents of the logfile are saved to the 'badRuns.txt' file in the getScore function.

	if len(spectomin) == 1:
	#For global error prediction (spectomin = ['all']), the final CAFE run (spectomin = ['final']), and individual species minimizations (spectomin = [the
	#current species] the spectomin list will contain just one element.
		if spectomin[0] == "none":
			printWrite(cafLog, 0, "## No error models generated");	
		else:
			for key in specDict:
				ErrString = "cafe_errormodel_" + str(specDict[key]) + ".txt";
				ErrFile = tempDir + ErrString

				genErrFile = errFileCheck(key, ErrString, specDict, tempDir, cafLog)
				#This line calls the errFileCheck function to see if the error model file for this error distribution and this species is already created.

				if genErrFile == 1:
					maxfamsize(families, specDict[key], ErrFile);
					#This line calls the maxfamsize function below to create the error file.

	# else:
	# #When -g is specified in the command line, spectomin will contain more than one element. [NOTE: -g not yet developed]
	# 	for key in spectomin:

	# 		ErrString = "cafe_errormodel_" + str(specDict[key]) + ".txt";
	# 		ErrFile = tempDir + ErrString
		
	# 		genErrFile = errFileCheck(key, ErrString, specDict, tempDir, cafLog)
	# 		#This line calls the errFileCheck function to see if the error model file for this error distribution and this species is already created.

	# 		if genErrFile == 1:
	# 			maxfamsize(families, specDict[key], ErrFile);
	# 			#This line calls the maxfamsize function below to create the error file.
	# # The above blocks create the error file for the current run (if necessary) by calling the errFileCheck and maxfamsize functions.

	#####
	printWrite(cafLog, 0, "-----");
	printWrite(cafLog, 1, "Rewriting CAFE shell script...");

	##########
	# These lines write the caferror.sh script.
	shellfile = "caferror.sh";

	# shellfileList = os.listdir(path);
	# shellcount = 1;
	# while shellfile in shellfileList:
	#	shellfile = "caferror" + str(shellcount) + ".sh";
	#	shellcount = shellcount + 1;
	# cLog.write("\nshellfile: " + shellfile + "\n");
	# cLog.write("logfile: " + logFile + "\n");
	# Some lines for debugging. Uncomment as desired.

	cafeFile = open(shellfile, "w");

	cafeFile.write(path + "\n" + newick + "\n");

	if fflag == 0:
		loadLine = "load -i " + families + " -t 10 -l " + logFile + "\n"
	elif fflag == 1:
		loadLine = "load -i " + families + " -t 10 -l " + logFile + " -filter\n"

	cafeFile.write(loadLine);

	if spectomin[0] != "none":
		if len(spectomin) == 1:
			for key in specDict:
				ErrString = "cafe_errormodel_" + str(specDict[key]) + ".txt";
				ErrFile = tempDir + ErrString;
				errLine = "errormodel -model " + ErrFile + " -sp " + key + "\n";
				cafeFile.write(errLine);	
	
		else:
			for key in spectomin:
				ErrString = "cafe_errormodel_" + str(specDict[key]) + ".txt";
				ErrFile = tempDir + ErrString;
				errLine = "errormodel -model " + ErrFile + " -sp " + key + "\n";
				cafeFile.write(errLine);	

	cafeFile.write(lamb + "\n");
	if spectomin[0] == "final":
		repline = "report " + repFile;
		cafeFile.write(repline);

	cafeFile.close();
	# End caferror.sh writing.
	##########

	os.system("chmod +x " + shellfile);
	if vb == 0:
		printWrite(cafLog, 1, "Running CAFE [silently] with error models listed above...");
		os.system("./" + shellfile + " >> " + tempDir + "cafe.out 2>&1");
	else:
		printWrite(cafLog, 1, "Running CAFE with error models listed above...");
		os.system("./" + shellfile);
	nr = nr + 1;
	#These lines give the caferror.sh script permissions and execute it to run cafe given the current parameters.

	print("CAFE run complete! Retrieving Score..........");
	#This is printed here. Score retrieval is actually handled by getScore (directly below), which is always called immediately after 
	#this function in the main block.

	return nr;

############################################
def errFileCheck(errspec_check, errstring_check, specdict_check, thedir, cafLog):
#This function checks if, for a given species and error model, the error model file has already been created. It returns 1 the error model file
#has not yet been created and 0 if it already has.

	fileList = os.listdir(thedir);

	if errstring_check in fileList:
		printWrite(cafLog, 0, "## " + str(specdict_check[errspec_check]) + " Error File Already Created for " + errspec_check);
		return 0;

	else:
		printWrite(cafLog, 0, "## Generating " + str(specdict_check[errspec_check]) + " Error File for " + errspec_check);	
		return 1;

############################################
def getScore(error, spectomin, tempDir, wout, cafLog):
#This function retrieves the score calculated by CAFE from a given log file, specified by the amount of error and species used in that run.
#This function now also checks if the previous CAFE run was initialized properly. If so, the program will continue and if not it will send
#a signal back to the call telling it to re-run CAFE with the same parameters.

	init = "good";
	lcount = 0;
	linfcount = 0;

	if spectomin[0] == "all":
		scoreFile = tempDir + "cafe_" + str(error) + "_all_log.txt";
	elif spectomin[0] == "final":
		scoreFile = tempDir + "cafe_final_log.txt";
	elif spectomin[0] == "none":
		scoreFile = tempDir + "cafe_pre_log.txt";
	else:
		scoreFile = tempDir + "cafe_" + str(error) + "_" + spectomin[0] + "_log.txt";

	sFile = open(scoreFile, "r");
	sLines = sFile.readlines();
	sFile.close();

	initFile = open(tempDir + "InitFile.txt", "a");
	initFile.write(scoreFile);
	initFile.write("\n")

	for s in xrange(len(sLines)):
		if sLines[s][:7] == "Poisson":
			initFile.write(sLines[s]);
		#if sLines[s][:7] == "Poisson" and sLines[s].find("inf") != -1:
		#	init = "bad";
		if sLines[s][:7] == ".Lambda":
			lcount = lcount + 1;
		if sLines[s][:7] == ".Lambda" and sLines[s].find("inf") != -1:
			linfcount = linfcount + 1;

		if sLines[s].find("Lambda Search Result:") != -1:
			lamval = sLines[s+1][9:sLines[s+1].index("&") - 1];
			score = sLines[s+1][sLines[s+1].index("Score: ") + 7:].replace("\n","");
			#score = score.replace('\n', '');

	if lcount == linfcount:
		init = "bad";

	initFile.write(str(score));
	initFile.write("\n");
	initFile.write(init);
	initFile.write("\n\n");
	initFile.close();


	printWrite(cafLog, 1, "Score with above error models:", str(score), 35);
	printWrite(cafLog, 1, "Lambda with above error models:", str(lamval), 35);


	if init == "bad":
		printWrite(cafLog, 1, "++WARNING: CAFE failed to converge or initialize. This run will be re-done. Check badRuns.txt for more info.\n");

		bFile = open(tempDir + "badRuns.txt", "a");
		bFile.write("\n************\n\n")
		bFile.write(scoreFile);
		bFile.write("\n");
		for sline in sLines:
			bFile.write(sline);
		bFile.close();

	if init == "good":
		if wout == 1:
			outFile = open(outFilename, "a");
			line = str(error) + "\t" + score + "\n"
			outFile.write(line);
			outFile.close();

	return float(score), lamval, init;

############################################
def maxfamsize(inFilename, totError, errFilename):
# This function creates an error file in the proper format for CAFE by reading the input gene family file and
# finding the gene family with the largest number of genes.

	negAsym = 0.5;
	posAsym = 1 - negAsym;
	# I used this to run some simulations with asymmetric error distributions. Right now it is set to split the error evenly
	# between +1 and -1, but go ahead and change negAsym if needed.

	inFile = open(inFilename, "r");
	lines = inFile.readlines();
	inFile.close();

	maxfs = 0;

	for i in range(len(lines)):
		if i == 0:
			continue;		
		cline = lines[i].replace("\n","").split("\t");	
		k = 2;
	
		while k <= len(cline) - 1:
			if int(cline[k]) > maxfs:
				maxfs = int(cline[k]);			
			k = k + 1;

	erroutFile = open(errFilename, "w");
	erroutFile.write("maxcnt:");
	erroutFile.write(str(maxfs));
	erroutFile.write("\n");
	erroutFile.write("cntdiff -1 0 1\n");

	j = 0;

	while j <= maxfs:
		if j == 0:
			pline = str(j) + " 0.00 " + str((totError / 2) + (1 - totError)) + " " + str(totError / 2);	
			erroutFile.write(pline);
		else:
			pline = str(j) + " " + str(totError * negAsym) + " " + str(1 - totError) + " " + str(totError * posAsym);		
			erroutFile.write(pline);

		if j != maxfs:
			erroutFile.write("\n");
		j = j + 1;
	
	erroutFile.close();

############################################
#Main Block
############################################

startsec = time.time();
start = datetime.datetime.now().strftime("%m.%d.%Y | %I:%M:%S");

inFilename, outFilename, initError, tmpDir, wholeCurveOpt, errTries, caferrorLog, indSpecMin, Mode, vOpt, firstRun = optParse(0);
# The first step is to retrieve the proper values for options and filenames based on the user's command line specifications.
# errTries = [0.0, 0.001, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95];

initCheck = "";
spec_to_min = [];
rmlog = 0;
FilterFlag = 0;

inFile = open(inFilename, "r");
inLines = inFile.readlines();
inFile.close();

for x in range(len(inLines)):
	if inLines[x][:2] == "#!":
		CafePath = inLines[x].replace("\n","");

	if inLines[x][:4] == "load":
		FamFile = inLines[x][(inLines[x].index("-i") + 3):];
		FamFile = FamFile[:FamFile.index(" ")];

		if inLines[x].find("-filter") != -1:
			FilterFlag = 1;

	if inLines[x][:4] == "tree":
		Tree = inLines[x].replace('\n', '');

	if inLines[x][:6] == "lambda":
		LamStruct = inLines[x].replace('\n', '');

		if "-t" in LamStruct:
			t_struct = LamStruct[LamStruct.index("-t"):];
			if not all(c in '-t \t(1,);' for c in t_struct):
				errorOut(5, "Only single lambda searches are recommended for estimating error. For more info, see the CAFE Manual.");
				sys.exit();
# Obtains input variables from an existing CAFE script specified with -i.

treestring = Tree[Tree.index("("):];
treestring = re.sub('[)][\d.eE-]+:[\d.eE-]+', ')', treestring);
treestring = re.sub(':[\d.eE-]+', '', treestring);
treestring = re.sub('<[\d]+>', '', treestring);
errSpec = treestring.replace("(","").replace(")","").replace(";","").split(",");
# This block extracts the species names from the input phylogeny.

mainSpecDict = {};
for eachspec in errSpec:
	mainSpecDict[eachspec] = initError;
spec_to_min.append("all");
# Some initializations of important variables above. 
# errSpec gets a list of all species in the input tree. 
# Throughout the program, mainSpecDict will keep track of the species and their corresponding minimized errors. It follows the [key]:[value] format
# of [species]:[minimized error model].
# spec_to_min is what will be passed as the species to minimize for the global error estimation algorithm. Usually this will be passed simply 
# as ['all'].

#########
# This block creates a new default directory for each run of caferror.py
if tmpDir == "":
	prevList = os.listdir(os.getcwd());

	isDirMade = 0;
	x = 1

	while isDirMade == 0:
		dirName = "caferror_" + str(x);
		if dirName not in prevList:		
			tmpDir = os.getcwd() + "/" + dirName + 	"/";
			isDirMade = 1;
		else:
			x = x + 1;

	if tmpDir.find("//") != -1:
		tmpDir = tmpDir.replace("//", "/");

	os.system("mkdir " + tmpDir);
	# If the user didn't define a temp directory with -d, the default value is set here. This also handles a possible error with the double slashes.

else:
	os.system("mkdir " + tmpDir);

# End default dir block
##########

caferrorLog = tmpDir + caferrorLog;
if caferrorLog.find("//") != -1:
	caferrorLog = caferrorLog.replace("//", "/");

clogfile = open(caferrorLog, "w");
clogfile.write("");
clogfile.close();
# If the user didn't define a log file with -l, the default value is set here. This also handles a possible error with the double slashes.

outFilename = tmpDir + outFilename;
outFile = open(outFilename, "w");
outFile.write("ErrorModel\tScore\n");
outFile.close();
# This primes the output file with the column headers.

####################################
### Begin input info block!
pad = 35

printWrite(caferrorLog, 1, "# =========================================================================");
printWrite(caferrorLog, 1, "#\t\tAssembly/annotation error estimation");
printWrite(caferrorLog, 1, "#\t\t\t" + str(start));
printWrite(caferrorLog, 1, "# Using CAFE shell file:", inFilename, pad);
printWrite(caferrorLog, 1, "# --------------------------");
printWrite(caferrorLog, 1, "#\t\t\tINPUT INFO");
printWrite(caferrorLog, 1, "# CAFE path set as:", CafePath, pad);
printWrite(caferrorLog, 1, "# Using gene family file:", FamFile, pad);
printWrite(caferrorLog, 1, "# Using tree command:", Tree, pad);
printWrite(caferrorLog, 1, "# Using lambda command", LamStruct, pad);
printWrite(caferrorLog, 1, "# --------------------------");
printWrite(caferrorLog, 1, "#\t\t\tOPTIONS INFO");
printWrite(caferrorLog, 1, "# -c " + str(wholeCurveOpt), "Global grid search option.", pad);
printWrite(caferrorLog, 1, "# -e " + str(initError), "Starting error estimate of " + str(initError) + ".", pad);
printWrite(caferrorLog, 1, "# -f " + str(firstRun), "Initial CAFE run option.", pad);
printWrite(caferrorLog, 1, "# -v " + str(vOpt), "CAFE verbosity option.", pad);
printWrite(caferrorLog, 1, "# -s " + str(indSpecMin), "Individual species minimization option.", pad);
printWrite(caferrorLog, 1, "# --------------------------");
printWrite(caferrorLog, 1, "#\t\t\tOUTPUT INFO");
printWrite(caferrorLog, 1, "# Putting all files in:", tmpDir, pad);
printWrite(caferrorLog, 1, "# Log file for all runs:", caferrorLog, pad);
printWrite(caferrorLog, 1, "# Output file:", outFilename, pad);
printWrite(caferrorLog, 1, "# =========================================================================");
### End input info block!
####################################

printWrite(caferrorLog, 1, "# Beginning Global Error Prediction...");

numruns = 1;

if firstRun == 1:
	spec_to_min[0] = "none";
	while initCheck == "bad" or initCheck == "":
		numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
		lastScore, lastLam, initCheck = getScore(initError, spec_to_min, tmpDir, 1, caferrorLog);
		if initCheck == "bad":
			rmlog = 1;
	initCheck = "";
	rmlog = 0;
	preScore = lastScore;
	preLam = lastLam;
	spec_to_min[0] = "all";
#The pre-error model run if -f 1.

####################################
#Global error estimation begins here for -c 0 (default)

if wholeCurveOpt == 0:
#If -c is set to 0, which is the default value, caferror enters this block to predict error in an iterative fashion that uses information from the 
#previous run until a minimum score is achieved. Caferror starts at the error value specified by -e (0.4 by default).
	errList = [];
	errList.append(initError);

	errMin = initError;
	#Initializations of some variables. The 'Holder' variables are used to keep the place of parameters which are un-needed for any one cafeRun call.

	while initCheck == "bad" or initCheck == "":
		numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
		lastScore, lastLam, initCheck = getScore(initError, spec_to_min, tmpDir, 1, caferrorLog);
		if initCheck == "bad":
			rmlog = 1;
	initCheck = "";
	rmlog = 0;
	#The initial Cafe run using the error valued specified by -e (0.4 by default).

	minScore = lastScore;
	errMin = initError;
	lastError = initError;

	nextError = 0.0;
	posLimit = 1.0;
	posLimScore = lastScore;
	negLimit = 0.0;
	negLimScore = 0.0; 

	currentError = initError / 2;
	errList.append(currentError);

	for eachspec in mainSpecDict:
		mainSpecDict[eachspec] = currentError;

	tally = 0;
	keepGoing = 1;

	while keepGoing == 1:
	#This will keep guessing error models until one of two termination scenarios are met. These scenarios are checked for at the
	#end of each CAFE run below.

		tally = tally + 1;

		errList.append(currentError);
		while initCheck == "bad" or initCheck == "":
			numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
			currentScore, currentLam, initCheck = getScore(currentError, spec_to_min, tmpDir, 1, caferrorLog);
			if initCheck == "bad":
				rmlog = 1;
		initCheck = "";
		rmlog = 0;

#		outFile = open(outFilename, "a");
#		line = str(currentError) + "\t" + str(currentScore) + "\n"
#		outFile.write(line);
#		outFile.close();
#		print("++++++++++";
#		print errMin, minScore;
#		print negLimit, posLimit;
#		print("++++++++++";
#		cLog.write("errMin and minScore before logic: ");
#		cLog.write(str(errMin));
#		cLog.write(", ");
#		cLog.write(str(minScore));
#		cLog.write("\n");
#		cLog.write("negLimit and posLimit before logic: ");
#		cLog.write(str(negLimit));
#		cLog.write(", ");
#		cLog.write(str(posLimit));
#		cLog.write("\n");
		#Some print statements to track values for debugging. Uncomment as desired.

		##################
		#The main logic statements of the code are below. These take the score and error values of the most recent CAFE run
		#and evaluate them against some information obtained from previous runs, eventually arriving at a point 
		#where the termination scenarios are met with a minimum score value.
		if currentScore < minScore:
			if currentError < lastError:
				nextError = currentError / 2;

				posLimit = lastError * 2;
				posLimScore = lastScore;

			elif currentError > lastError:
				nextError = currentError * 2;

				negLimit = lastError / 2;
				megLimScore = lastScore;

			minScore = currentScore;	
			minLam = currentLam;
			errMin = currentError;

		elif currentScore > minScore:
			if currentError < errMin:
				nextError = errMin * 2

				negLimit = currentError;
				negLimScore = currentScore;

			elif currentError > errMin:
				nextError = errMin / 2;

				posLimit = currentError;
				posLimScore = currentScore;

#		print("----------";
#		print errMin, minScore;
#		print negLimit, posLimit;
#		print("----------";
#		cLog.write("errMin and minScore after logic: ");
#		cLog.write(str(errMin));
#		cLog.write(", ");
#		cLog.write(str(minScore));
#		cLog.write("\n");
#		cLog.write("negLimit and posLimit after logic: ");
#		cLog.write(str(negLimit));
#		cLog.write(", ");
#		cLog.write(str(posLimit));
#		cLog.write("\n");
		#Some print statements to track values for debugging. Uncomment as desired.
	
		if nextError > posLimit:
			nextError = ((currentError + posLimit) / 2);
		if nextError < negLimit:
			nextError = ((currentError + negLimit) / 2);
		if nextError in errList:
			nextError = ((currentError + nextError) / 2);

#		print("***************";
#		print currentError;
#		print len(str(currentError));
#		print("***************";
		#Some print statements to track values for debugging. Uncomment as desired.

		lastError = currentError;

		#####
		#These are the termination scenarios. If the error parameter reaches a certain level of precision (12 decimal places), or if
		#CAFE has been run 20 times, error estimation is terminated. This may be improved.
		if len(str(currentError)) >= 15:
			keepGoing = 0;
		elif tally > 20:
			keepGoing = 0;
		#End termination block.
		#####
		else:
			currentError = nextError;
			for eachspec in mainSpecDict:
				mainSpecDict[eachspec] = nextError;
			lastScore = currentScore;
		#End logic block.
		##################

#Global error estimation ends here for -c 0 (default)
####################################



####################################
#Global error estimation begins here for -c 1

elif wholeCurveOpt == 1:
#If -c is set to 1 by the user caferror enters this block to calculate scores for a pre-set set of values [0,1).
#This will be less precise than the method above, but will give the user an easy to visualize output when plotted.

	minScore = 1000000000.0;
	errMin = 1.0;
	#Initializations of some variables.


	for currentError in errTries:
	#This runs Cafe using all error models in the list 'errTries' defined above. Feel free to add or remove values as desired.

		for eachspec in mainSpecDict:
			mainSpecDict[eachspec] = currentError;
	
		while initCheck == "bad" or initCheck == "":
			numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
			currentScore, currentLam, initCheck = getScore(currentError, spec_to_min, tmpDir, 1, caferrorLog);
			if initCheck == "bad":
				rmlog = 1;
		initCheck = "";
		rmlog = 0;

		if float(currentScore) < minScore and minScore != float("infinity"):
			minScore = float(currentScore);
			minLam = currentLam;
			errMin = currentError;


#Global error estimation ends here for -c 1
####################################


####################################
#Individual species error estimation begins here if -s 1
#
#After the above algorithm has minimized error scores across all species (by applying the same error model to all species simultaneously), the user has
#the option to continue minimization for each individual species. This is done by iteratively adding and subtracting 10% of the global min error from each
#current species's error model until the score ceases to improve.

#errMin = 0.04765625;
#minScore = 50384.197476;
#minLam = 0.27626649636157;

finalSpecDict = {};

if Mode == 2 or Mode == 3:
	random.shuffle(errSpec);

if Mode == 1 or Mode == 3:
	for eachspec in mainSpecDict:
		mainSpecDict[eachspec] = errMin;
	overallMinScore = minScore;

for eachspec in mainSpecDict:
	finalSpecDict[eachspec] = 0.0

if indSpecMin == 1:
	printWrite(caferrorLog, 1, "***********************************************************************");
	printWrite(caferrorLog, 1, "Global error prediction complete. Beginning error minimization for individual species.");

	moreErrSpec = [];

	for curSpec in errSpec:
		printWrite(caferrorLog, 1, "--------------------------");
		printWrite(caferrorLog, 1, "**Attempting to Minimize " + curSpec + " by adding more error");

		if Mode == 0 or Mode == 2:
			specMinScore = minScore;
		if Mode == 1 or Mode == 3:
			specMinScore = overallMinScore;
		oldspecMin = minScore;

		zcheck = 0;
		minimized = 0;
		currentError = errMin + (0.1 * errMin);
		#Some initializations.

		if Mode == 1 or Mode == 3:
			specMinErr = errMin;

		if Mode == 0 or Mode == 2:
			for eachspec in mainSpecDict:
				mainSpecDict[eachspec] = errMin;

		spec_to_min[0] = curSpec;

		while minimized == 0:
		#For each species, 10% of the global min error will be iteratively added in this loop until the score ceases to decrease OR the error has reached 1.0. At which point minimized will be set
		#to 1 and the loop will exit.

			mainSpecDict[curSpec] = currentError;

			if mainSpecDict[curSpec] <= 1.0:
				while initCheck == "bad" or initCheck == "":
					numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
					currentScore, currentLam, initCheck = getScore(currentError, spec_to_min, tmpDir, 0, caferrorLog);
					if initCheck == "bad":
						rmlog = 1;
				initCheck = "";
				rmlog = 0;

			else:
				minimized = 1;
				continue;

			if float(currentScore) > float(specMinScore):
				check = "higher than";
				check2 = "will not be further minimized.";

				if finalSpecDict[curSpec] == 0.0:
					finalSpecDict[curSpec] = errMin;

				oldspecMin = specMinScore;
				minimized = 1;

			elif float(currentScore) == float(specMinScore):
				check = "the same as";
				check2 = "will not be further minimized.";
				oldspecMin = specMinScore;
				minimized = 1;

			elif float(currentScore) < float(specMinScore):
				check = "lower than";
				check2 = "will continue to be minimized...";

				if finalSpecDict[curSpec] == 0.0:
					finalSpecDict[curSpec] = errMin;

				oldspecMin = specMinScore;
				specMinScore = currentScore;

				if zcheck == 1:
					minimized = 1;

				finalSpecDict[curSpec] = currentError;

				specMinErr = currentError;

### 				if Mode == 1 or Mode == 3:
###					overallMinScore = currentScore;
###				Toggle this line to either keep the background constant (commented) or update it each time a species is minimized.

				currentError = currentError + (0.1 * errMin);

				if currentError >= 1:
					currentError = 1.0
					zcheck = 1;

				if curSpec not in moreErrSpec:
					moreErrSpec.append(curSpec);

		if Mode == 1 or Mode == 3:
			mainSpecDict[curSpec] = specMinErr;
			
	for curSpec in errSpec:
	#This block iteratively subtracts 10% of the global min error to each species' error model to attempt to find out which species contain LESS error than
	#the global min, and to get a general guess of what that error might be.
		if curSpec not in moreErrSpec:
			printWrite(caferrorLog, 1, "--------------------------");
			printWrite(caferrorLog, 1, "**Attempting to Minimize " + curSpec + " by adding less error");

			if Mode == 0 or Mode == 2:
				specMinScore = minScore;
			if Mode == 1 or Mode == 3:
				specMinScore = overallMinScore;
			oldspecMin = minScore;

			zcheck = 0;
			minimized = 0;
			currentError = errMin - (0.1 * errMin);

			if Mode == 1 or Mode == 3:
				specMinErr = errMin;

			if Mode == 0 or Mode == 2:
				for eachspec in mainSpecDict:
					mainSpecDict[eachspec] = errMin;

			spec_to_min[0] = curSpec;

			while minimized == 0:
			# For each species, 10% of the global min error will be iteratively subtracted in this loop until the score ceases to decrease OR the error has reached 0.0.
			# At which point minimized will be set to 1 and the loop will exit.
				mainSpecDict[curSpec] = currentError;
				
				if mainSpecDict[curSpec] >= 0.0:
					while initCheck == "bad" or initCheck == "":
						numruns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
						currentScore, currentLam, initCheck = getScore(currentError, spec_to_min, tmpDir, 0, caferrorLog);
						if initCheck == "bad":
							rmlog = 1;
					initCheck = "";
					rmlog = 0;
				else:
					minimized = 1;
					continue;
				
				if float(currentScore) > float(specMinScore):
					check = "higher than";
					check2 = "will not be further minimized.";

					if finalSpecDict[curSpec] == 0.0:
						finalSpecDict[curSpec] = errMin;

					oldspecMin = specMinScore;
					minimized = 1;

				elif float(currentScore) == float(specMinScore):
					check = "the same as";
					check2 = "will not be further minimized.";

					if finalSpecDict[curSpec] == 0.0:
						finalSpecDict[curSpec] = errMin;

					oldspecMin = specMinScore;
					minimized = 1;

				elif float(currentScore) < float(specMinScore):
					check = "lower than";
					check2 = "will continue to be minimized...";
					oldspecMin = specMinScore;
					specMinScore = currentScore;

					if zcheck == 1:
						minimized = 1;

					finalSpecDict[curSpec] = currentError;

## #					if Mode == 1 or Mode == 3:
					specMinErr = currentError;
###					Toggle this line to either keep the background constant (commented) or update it each time a species is minimized.

		
					currentError = currentError - (0.1 * errMin);
	
					if currentError <= 0:
						currentError = 0.0;
						zcheck = 1;
					
			if Mode == 1 or Mode == 3:
				mainSpecDict[curSpec] = specMinErr;

# Individual species error estimation ends here if -s 1
####################################
	#####
	# The final CAFE run which generates a report with the minimized error.
	printWrite(caferrorLog, 0, "\n***********************************************************************");

	spec_to_min[0] = "final";

	while initCheck == "bad" or initCheck == "":
		numrnuns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, finalSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
		finalScore, finalLam, initCheck = getScore("final", spec_to_min, tmpDir, 0, caferrorLog);
		if initCheck == "bad":
			rmlog = 1;
	initCheck = "";
	rmlog = 0;
	#####
####################################
if indSpecMin == 1:
# Output block if -s 1
	printWrite(caferrorLog, 1, "# =======================================================================");
	printWrite(caferrorLog, 1, "# Final error estimates by species:");
	
	for key in finalSpecDict:
		printWrite(caferrorLog, 1, "# " + key, str(finalSpecDict[key]), pad);

	printWrite(caferrorLog, 1, "# Score with individual errors:", str(finalScore), pad);
	printWrite(caferrorLog, 1, "# Lambda with individual errors:", str(finalLam), pad);

else:
# Output block if -s 0
	spec_to_min[0] = "final";

	for eachspec in mainSpecDict:
		mainSpecDict[eachspec] = errMin;

	while initCheck == "bad" or initCheck == "":
		numrnuns = cafeRun(FamFile, Tree, LamStruct, CafePath, tmpDir, caferrorLog, mainSpecDict, spec_to_min, rmlog, FilterFlag, numruns, vOpt);
		finalScore, finalLam, initCheck = getScore("final", spec_to_min, tmpDir, 0, caferrorLog);
		if initCheck == "bad":
			rmlog = 1;
	initCheck = "";
	rmlog = 0;

####################################
# Main output block
printWrite(caferrorLog, 1, "# =======================================================================");
if firstRun == 1:
	printWrite(caferrorLog, 1, "# ************************************");
	printWrite(caferrorLog, 1, "# Score with no errormodel:", str(preScore), pad);
	printWrite(caferrorLog, 1, "# Lambda with no errormodel:", str(preLam), pad);
printWrite(caferrorLog, 1, "# ************************************");
printWrite(caferrorLog, 1, "# Global Error Estimation:", str(errMin), pad);
printWrite(caferrorLog, 1, "# Score with global errormodel:", str(minScore), pad);
printWrite(caferrorLog, 1, "# Lambda with global errormodel:", str(minLam), pad);
printWrite(caferrorLog, 1, "# ************************************");
printWrite(caferrorLog, 1, "# =======================================================================");

endsec = time.time();
runtime = (float(endsec) - float(startsec)) / float(60.0);
end = datetime.datetime.now().strftime("%m.%d.%Y | %I:%M:%S");

printWrite(caferrorLog, 1, "# Caferror finished at:", str(end), pad);
printWrite(caferrorLog, 1, "# Runtime:", str(runtime) + " minutes", pad);
#############################################################################
#End Main Block
#############################################################################