#include "../headers/results.h"

using namespace std;


/**-----------------------------------------------------------------------------
	## FUNCTION:
	double ** init_tab_ae(s_params *params, vector<string> & experiments)
	-----------------------------------------------------------------------------
	## RETURN:
	@ double ** tab: 2D array contening the allelic effect value
	for each qtl (line) and each parental allele (col). nqtls x nfounders
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_params *params : structure containing all data parameters
	@ vector<string> & experiments : vector storing the "experiences" names (e.g. TixEj)
   -----------------------------------------------------------------------------
   ## SPECIFICATION: read the file "allelic_effects.txt", store the values in
   the returned array, and keep the "experiences" names in a vector.
   -----------------------------------------------------------------------------
*/


double ** init_tab_ae(s_params *params, vector<string> & experiments){

	double ** tab = NULL;
	
	char file_ae[NB_PATH_MAX]; //// allelic effects file name (TODO: change only the extension to .ae)
	strcpy(file_ae, params->file_mks);
	
	char* ext = NULL;
	// get the input directory path : an empty string if current directory, "/path/to/input/" otherwise. 
	if ( (ext=strrchr(file_ae, PATH_SEP))==NULL ) { file_ae[0]='\0'; ext = file_ae; }
	*(ext+1) = '\0';

	// Allelic effects file does not exist
	if ( access(strcat(file_ae, "allelic_effects.txt"), F_OK) != 0 ) return NULL;
	// open allelic effects file in readonly mode (default)
	ifstream fae(file_ae);
	if(!fae) { fprintf(stderr,"Warning ! allelic effects file could not be read.\n"); fflush(stderr); return NULL; }

	char line[NB_LINE_MEAN]; // line buffer
	
	// header #########################
	fae.getline(line, NB_LINE_MEAN);
	chomp(line);
	char** ae_fields = str_split(line,"\t");
	// skip the file if the header is empty ?
	if(ae_fields == NULL || ae_fields[0] == NULL) { free(ae_fields); fae.close(); return NULL; }
	int fields_nb = count(ae_fields);
	experiments.clear(); // vector storing the "experience" names (e.g. TixEj)
	experiments.insert(experiments.begin(), ae_fields+1, ae_fields+fields_nb);
	//int exp_nb = experiments.size();

	if(strcmp(ae_fields[0], "QTL,allele") != 0) {
		fprintf(stderr,"Error in allelic effects file:\nThe header first column has to be QTL,allele.\n");
		fflush(stderr);	exit(1);
	}
	free(ae_fields);
	ae_fields = NULL;
	// end header #########################

	// the order in which parental allele occur in the allelic file may not follow the order seen in the genotype file.
	// In tab, it has to.
	int l_nb = params->nb_qtl * params->nb_founders;
	calloc2d(tab, l_nb, fields_nb-1, double); // default value is 0

	int cpt = 0;
	char* err = NULL;
	char** qtl_allele = NULL;
	map<string, geno_coded>::const_iterator it; // iterator on params->str2coded
/*
 * allelic_effects example file.
QTL,allele	T1xE1	T1xE2	T2xE1      line index of each parental allele in str2coded (d:0, f:1, s:2, x:3)
QTL1,d	0.84	0.84	0.84           0  	0
QTL1,f	0.56	0.56	0.56           1  	1
QTL1,s	-0.36	-0.36	-0.36          2  	2
QTL1,x	-1.04	-1.04	-1.04          3  	3
...
QTL4,s	1.26	1.26	1.26          12  	2
QTL4,f	-0.79	-0.79	-0.79         13  	1
QTL4,x	-0.13	-0.13	-0.13         14  	3
QTL4,d	-0.34	-0.34	-0.34         15  	0
*/		
	while(!fae.getline(line,NB_LINE_MEAN).eof()){

		chomp(line); if(*line == '\0') continue; // ignore empty line
		ae_fields = str_split(line,"\t");
		
		int n = count(ae_fields);
		if (n != fields_nb){
			fprintf(stderr,"Error in allelic effects file:\nline %d, the number of columns differes from that in the header\n", cpt+2);
			fflush(stderr);
			exit(1);
		}
		// split QTL name and allele (first column);
		qtl_allele = str_split(ae_fields[0],",");
		
		if( (it=params->str2coded.find(qtl_allele[1])) == params->str2coded.end() ){
			fprintf(stderr,"Error in allelic effects file:\nline %d, '%s' does not match any parental allele.\n", cpt+2, qtl_allele[1]);
			fflush(stderr);
			exit(1);
		}
		//cout<<"cpt: "<<cpt<<" ; itab: "<<cpt-cpt%params->nb_founders + it->second<<" ; index of "<<it->first<<":  "<<(int)(it->second)<<endl; // it->second is the index of this parental allele in str2coded	
		int itab = cpt - cpt%params->nb_founders + it->second; // it->second is the index of this parental allele in str2coded	
		for(int j=1; j<n; j++){
			tab[itab][j-1] = strtod(ae_fields[j], &err);
			if(*err != '\0') {
				if(strcmp(ae_fields[j],"-") == 0) tab[itab][j-1] = 0.0;
				else {
					fprintf(stderr, "Error: in allelic effects file:\nline %d : the value (%s) is incorrect.\n", cpt+2, ae_fields[j]);
					fflush(stderr);
					exit(1);
				}
			}
		}
		cpt++;
		free(ae_fields); ae_fields = NULL;
		free(qtl_allele); qtl_allele = NULL;
	}
	fae.close();

	return tab;
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void crea_results_final_tabs(s_indiv *tab_indiv, s_params *params)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_indiv *tab_indiv : array of individuals
	@ s_params *params : structure containing all data parameters
   -----------------------------------------------------------------------------
   ## SPECIFICATION: All the result files/tables are written (Molecular Score,
   homozygous favorable/unfavorable heterozygous, founders alleles at all/each
   QTL for each individual.
   -----------------------------------------------------------------------------
*/

void crea_results_final_tabs(s_indiv *tab_indiv, s_params *params){

	vector<string> experiments;
	//double ** tab_ae = init_tab_ae(params, experiments); // Allelic effects are loaded from the allelic_effects.txt file
	double ** tab_ae = NULL; // allelic effects stuff is now done in the GUI.

	// Predicted Molecular Score(s) and phenotype(s) strings concatenation 
	string str_pms("");
	if(tab_ae != NULL){
		for(size_t iexp=0; iexp<experiments.size(); iexp++) str_pms += "\t" + experiments[iexp];
	}
	string str_pheno("");
	for(size_t ipheno=0; ipheno<params->vect_pheno_name.size(); ipheno++) str_pheno += "\t" + params->vect_pheno_name[ipheno];

	int index_cycles = 2;
	double uc = 0.0; // Utility criterion

	char name_file_scores[NB_PATH_MAX], name_file_parents[NB_PATH_MAX], name_file_homo_hetero[NB_PATH_MAX], name_file_check_haplo[NB_PATH_MAX];

	sprintf(name_file_scores, "%s/tab_scores.txt", params->file_root);
	sprintf(name_file_parents, "%s/tab_parents.txt", params->file_root);
	sprintf(name_file_homo_hetero, "%s/tab_homo_hetero.txt", params->file_root);
	sprintf(name_file_check_haplo, "%s/tab_check_diplo.txt", params->file_root);

	FILE *file_scores = open_file(name_file_scores, "w");
	FILE *file_parents = open_file(name_file_parents, "w");
	FILE *file_homo_hetero = open_file(name_file_homo_hetero, "w");
	FILE *file_check_haplo = open_file(name_file_check_haplo, "w");

	// Headers
	for(int i=0; i<params->nb_qtl; i++){
		if(i == 0){
			//~ fprintf(file_scores, "Id\tP1\tP2\tCycle\tGroup\tMS\tWeight\tUC\tNo.(+/+)\tNo.(-/-)\tNo.(+/-)\tNo.(?)");
			//~ fprintf(file_scores, "Id\tP1\tP2\tCycle\tGroup%s\tMS\tWeight\tUC\tNo.(+/+)\tNo.(-/-)\tNo.(+/-)\tNo.(?)", str_pms.c_str()); // PMS (predicted molecular scores)
			fprintf(file_scores, "Id\tP1\tP2\tCycle\tGroup%s\tMS\tMS_Weight\tMS_UC\tNo.(+/+)\tNo.(-/-)\tNo.(+/-)\tNo.(?)", str_pms.c_str()); // PMS (predicted molecular scores)
			fprintf(file_homo_hetero, "Id\tCycle\tGroup\tMS\tAll(+/+)\tAll(-/-)\tAll(+/-)");
			fprintf(file_parents, "Id\tCycle\tGroup\tMS");
			fprintf(file_check_haplo, "Id\tMS\tAll");

			// First individuals are founders
			for(int j=0; j<params->nb_founders; j++){
				fprintf(file_parents, "\tAll(%s)", params->coded2str[j]); // ex indiv->parental_allele
			}
		}
		//~ fprintf(file_scores, "\tQTL%d", i + 1);
		//~ fprintf(file_check_haplo, "\tQTL%d", i + 1);
		//~ fprintf(file_homo_hetero, "\tQTL%d(+/+)\tQTL%d(-/-)\tQTL%d(+/-)", i + 1, i + 1, i + 1);
		//~ for(int j=0; j<params->nb_founders; j++){
			//~ fprintf(file_parents, "\tQTL%d(%s)", i + 1, params->coded2str[j]); // ex indiv->parental_allele
		//~ }
		const char* qtl_cstr = params->vect_qtl_names[i].c_str();
		fprintf(file_scores, "\t%s", qtl_cstr);
		fprintf(file_check_haplo, "\t%s", qtl_cstr);
		fprintf(file_homo_hetero, "\t%s(+/+)\t%s(-/-)\t%s(+/-)", qtl_cstr, qtl_cstr, qtl_cstr);
		for(int j=0; j<params->nb_founders; j++){
			fprintf(file_parents, "\t%s(%s)", qtl_cstr, params->coded2str[j]); // ex indiv->parental_allele
		}
	}
	fprintf(file_scores, "%s\n", str_pheno.c_str());
	fprintf(file_homo_hetero, "\n");
	fprintf(file_parents, "\n");
	fprintf(file_check_haplo, "\n");

	int tab_nb_homo_hetero[4] = {0,0,0,0};
	int tab_nb_homo_hetero_uc[3] = {0,0,0}; // The genotype with the higher probability (++/--/+-)
	double val_homo_hetero[3] = {0.0, 0.0, 0.0};
	vector<double> pms(experiments.size(), 0); // PMS (predicted molecular scores, computed from allelic effects table)

	for(int i=0; i<params->nb_indiv; i++){
		fprintf(file_scores, "%s\t%s\t%s\t%s\t%s", tab_indiv[i].name, tab_indiv[i].par1, tab_indiv[i].par2, tab_indiv[i].cycle, tab_indiv[i].group);
		fprintf(file_homo_hetero, "%s\t%s\t%s", tab_indiv[i].name, tab_indiv[i].cycle, tab_indiv[i].group);
		fprintf(file_parents, "%s\t%s\t%s", tab_indiv[i].name, tab_indiv[i].cycle, tab_indiv[i].group);
		fprintf(file_check_haplo, "%s", tab_indiv[i].name);
		
		memset(tab_nb_homo_hetero, 0, sizeof(tab_nb_homo_hetero)); // reset to zero
		memset(tab_nb_homo_hetero_uc, 0, sizeof(tab_nb_homo_hetero_uc));
		pms.assign(pms.size(), 0.0);

		// Computation of No.(+/+)/No.(-/-)/No.(+/-)/No(?) regarding the table of homo/hetero (no more based on MS)
		for(int j=1; j< (params->nb_qtl + 1); j++){
			// Fixed homozygous favorable
			if( params->tab_homo_hetero_final[i][j][0] >= params->cut_off_locus_fixed_fav ){
				tab_nb_homo_hetero[0] += 1;
			}
			// Fixed homozygous unfavorable
			else if ( params->tab_homo_hetero_final[i][j][1] >= params->cut_off_locus_fixed_unfav ){
				tab_nb_homo_hetero[1] += 1;
			}
			// Fixed heterozygous
			else if( params->tab_homo_hetero_final[i][j][2] >= params->cut_off_locus_fixed_hetero ){
				tab_nb_homo_hetero[2] += 1;
			}
			// Uncertain genotypes
			else{
				tab_nb_homo_hetero[3] += 1;
			}
			val_homo_hetero[0] = params->tab_homo_hetero_final[i][j][0];
			val_homo_hetero[1] = params->tab_homo_hetero_final[i][j][1];
			val_homo_hetero[2] = params->tab_homo_hetero_final[i][j][2];

			// Computation of the number of homo fav/unfav, hetero for the UC method (discretization via the genotype with the higher probability)
			tab_nb_homo_hetero_uc[maxIndex(val_homo_hetero, 3)] += 1;

			//printf("%f %f %f -> %d\n", val_homo_hetero[0], val_homo_hetero[1], val_homo_hetero[2], maxIndex(val_homo_hetero, 3));

			// if tab_ae exists, add the predicted molecular score (pms) of this qtl to each experiment pms (pms[iexp]). One dose.
			if(tab_ae != NULL){ 
				for(int k=0; k<params->nb_founders; k++){
					for(int iexp=0; iexp<(int)pms.size(); iexp++){
						//pms[iexp] += params->tab_founders_final[i][j][k] * tab_ae[j-1][k]; // i = indiv ; j = qtl ; k = parental allele
						int itab = (j-1)*params->nb_founders + k; // exemple avec blanc : j=5 et k=3 => qtl5 et allele x => itab = 4*4+3 = 19
						pms[iexp] += params->tab_founders_final[i][j][k] * tab_ae[itab][iexp]; // i = indiv ; j = qtl ; k = parental allele
					}
				}
			}
		}
		str_pms = "";
		stringstream ss;
		ss.setf(ios::fixed, ios::floatfield); // default float precision = 6
		if(tab_ae != NULL) {
			for(size_t iexp=0; iexp<pms.size(); iexp++) ss<<"\t"<<(pms[iexp]*=2); // 2 doses => x2
			str_pms = ss.str(); // 
		}
		
		for(int j=0; j< (params->nb_qtl + 1); j++){
			// Molecular score (All) and Weight (if classification with allelic affects)
			if( j == 0 ){
				// Utility criterion
				uc = 2 * params->tab_scores_final[i][j] + index_cycles * sqrt(tab_nb_homo_hetero_uc[2] * 0.25);
				// PMS, MS and weight
				fprintf(file_scores, "%s\t%f\t%f", str_pms.c_str(), params->tab_scores_final[i][j] / params->nb_qtl, params->tab_scores_final[i][j] / params->nb_qtl);
				fprintf(file_scores, "\t%f\t%d\t%d\t%d\t%d", uc / 2.0, tab_nb_homo_hetero[0], tab_nb_homo_hetero[1], tab_nb_homo_hetero[2], tab_nb_homo_hetero[3] );

				fprintf(file_check_haplo, "\t%f\t%f", params->tab_scores_final[i][j] / params->nb_qtl, params->tab_check_diplo[i][j] / params->nb_qtl);

				fprintf(file_homo_hetero, "\t%f\t%f\t%f\t%f", params->tab_scores_final[i][j] / params->nb_qtl, params->tab_homo_hetero_final[i][j][0] / params->nb_qtl, params->tab_homo_hetero_final[i][j][1] / params->nb_qtl, params->tab_homo_hetero_final[i][j][2]/ params->nb_qtl);
				fprintf(file_parents, "\t%f", params->tab_scores_final[i][j] / params->nb_qtl);

				for(int k=0; k<params->nb_founders; k++){
					fprintf(file_parents, "\t%f", params->tab_founders_final[i][j][k] / params->nb_qtl);
				}
			}
			else{
				fprintf(file_scores, "\t%f", params->tab_scores_final[i][j]);
				fprintf(file_check_haplo, "\t%f", params->tab_check_diplo[i][j]);
				fprintf(file_homo_hetero, "\t%f\t%f\t%f", params->tab_homo_hetero_final[i][j][0], params->tab_homo_hetero_final[i][j][1], params->tab_homo_hetero_final[i][j][2]);
				for(int k=0; k<params->nb_founders; k++){
					fprintf(file_parents, "\t%f", params->tab_founders_final[i][j][k]);
				}
			}
		}
		// Phenotypes
		for(int ipheno=0; ipheno<params->nb_pheno; ipheno++) fprintf(file_scores, "\t%f", tab_indiv[i].pheno[ipheno]);
		fprintf(file_scores, "\n");
		fprintf(file_homo_hetero, "\n");
		fprintf(file_parents, "\n");
		fprintf(file_check_haplo, "\n");
	}

	free_alloc2d(tab_ae);
	fclose(file_scores);
	fclose(file_parents);
	fclose(file_homo_hetero);
	fclose(file_check_haplo);
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void crea_summary_events(s_indiv *tab_indiv, s_params *params)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_indiv *tab_indiv : array of individuals
	@ s_params *params : structure containing all data parameters
   -----------------------------------------------------------------------------
   ## SPECIFICATION: questionable results entered in the events.log file are sorted
   by individuals regarding the QTL/mks that can have genotyping errors (if .cmap
   file is used).
   -----------------------------------------------------------------------------
*/

void crea_summary_events(s_indiv *tab_indiv, s_params *params){
	// Creation of the new events.log file (sorted by Id)
	char path_file_sum_log[NB_PATH_MAX];
	sprintf(path_file_sum_log, "%s/events_summary.log", params->file_root);

	FILE *file_sum_log = open_file(path_file_sum_log, "w");
	fprintf(file_sum_log, "SUMMARY [Questionable results at markers/QTL sorted by Id]\n--\n");

	for(int i=0; i<params->nb_indiv; i++){
		if( params->map_ind_qtl_events.find(string(tab_indiv[i].name)) != params->map_ind_qtl_events.end() ){
			fprintf(file_sum_log, "Id: %s\tP1: %s\t P2: %s\tCycle:%s\tQTL:", tab_indiv[i].name, tab_indiv[i].par1, tab_indiv[i].par2, tab_indiv[i].cycle);
			for(int j=0; j<(int)params->map_ind_qtl_events[tab_indiv[i].name].size(); j++){
				fprintf(file_sum_log, " %d", params->map_ind_qtl_events[tab_indiv[i].name][j]);
			}
			fprintf(file_sum_log, "\n");
		}
	}
	fclose(file_sum_log);
}
