/* io.cpp
Member function definitions for IO class
*/

#include <iostream>
#include <fstream>
using std::ifstream;
using std::ofstream;
using std::cout;
using std::endl;
using std::ios;

#include <stdexcept>
using std::runtime_error;
using std::out_of_range;

#include <string>
using std::string;

#include <vector>
using std::vector;

#include "io.h"
#include "coaltree.h"
#include "series.h"

IO::IO() {

	/* these will eventually be taken care of with a Parameters class */
	inputFile = "in.trees";
	outputPrefix = "out";

	ifstream inStream;
	inStream.open( inputFile.c_str(),ios::out);

	string line;
	int pos;
	if (inStream.is_open()) {
		while (! inStream.eof() ) {
			getline (inStream,line);
			
			if (line.size() > 0) {
				if (line[0] != '#') {
			
					// Catching log probabilities of trees
					string annoString;
					
					// migrate annotation
					annoString = "log likelihood = ";
					pos = line.find(annoString);
		
					if (pos >= 0) {
						string thisString;
						thisString = line.substr(pos+annoString.size());
						thisString.erase(thisString.find(']'));
						double ll = atof(thisString.c_str());
						problist.push_back(ll);
					}
					
					// beast annotation
					annoString = "[&lnP=";
					pos = line.find(annoString);
		
					if (pos >= 0) {
						string thisString;
						thisString = line.substr(pos+annoString.size());
						thisString.erase(thisString.find(']'));
						double ll = atof(thisString.c_str());
						problist.push_back(ll);				
					}
					
					// find first occurance of '(' in line
					pos = line.find('(');
					if (pos >= 0) {
					
						string paren = line.substr(pos);
						
						CoalescentTree ct(paren);
						ct.pushTimesBack(2007,2007);
						treelist.push_back(ct);
						cout << "tree " << treelist.size() << " read" << endl;
							
					}

				}
			}			
		}
		inStream.close();
	}
	else {
		throw runtime_error("input file not found");
	}
		
}

/* go through problist and treelist and print highest posterior probability tree */
void IO::printHPTree() {

	string outputFile = outputPrefix + ".rules";
	cout << "Printing tree with highest posterior probability to " << outputFile << endl;

	// output tree with highest likelihood
	// if likelihoods don't exist, output final tree
	
	int index;
	if (problist.size() == treelist.size()) {
		
		double ll = problist[0];
		index = 0;
		for (int i = 1; i < problist.size(); i++) {
			if (problist[i] > ll) {
				ll = problist[i];
				index = i;
			}
		}
		
	}
	else {
		index = treelist.size() - 1;
	}
	
	treelist[index].printRuleList(outputFile);

}

/* go through treelist and summarize coalescent statistics */
void IO::printStatistics() {

	string outputFile = outputPrefix + ".stats";
	cout << "Printing coalescent statistics to " << outputFile << endl;

	/* initializing output stream */
	ofstream outStream;
	outStream.open( outputFile.c_str(),ios::out);
	
	outStream << "statistic\t95% lower\tmean\t95% upper" << endl; 
	
	int n = treelist[0].getMaxLabel();

	// COALESCENCE /////////////////////
	for (int label = 1; label <= n; label++) {
	
		Series s;
		for (int i = 0; i < treelist.size(); i++) {
			double n = treelist[i].getCoalRate(label);
			s.insert(n);
		}
		outStream << "coal_" << label << "\t";
		outStream << s.quantile(0.025) << "\t" << s.mean() << "\t" << s.quantile(0.975) << endl;
	
	}
	
	// MIGRATION ///////////////////////
	for (int from = 1; from <= n; from++) {
		for (int to = 1; to <= n; to++) {	
			if (from != to) {

				Series s;
				for (int i = 0; i < treelist.size(); i++) {
					double n = treelist[i].getMigRate(from,to);
					s.insert(n);
				}
				outStream << "mig_" << from << "_" << to << "\t";
				outStream << s.quantile(0.025) << "\t" << s.mean() << "\t" << s.quantile(0.975) << endl;

			}
		}	
	}
	
	// TRUNK PROPORTIONS //////////////
	for (int label = 1; label <= n; label++) {

	Series s;
	for (int i = 0; i < treelist.size(); i++) {
		CoalescentTree ct = treelist[i];
		ct.pruneToTrunk();
		double n = ct.getLabelPro(label);
		s.insert(n);
	}
	outStream << "pro_" << label << "\t";
	outStream << s.quantile(0.025) << "\t" << s.mean() << "\t" << s.quantile(0.975) << endl;
	
	}
	
	
	outStream.close();

}