#include "../headers/markers.h"

using namespace std;

/**-----------------------------------------------------------------------------
	## FUNCTION:
	bool stl_comp_mrkpos(const mrkpos & m1, const mrkpos & m2)
	-----------------------------------------------------------------------------
	## RETURN:
	bool : true if m1.pos < m2.pos, false otherwise.
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ const mrkpos & m1 : a const reference to the 1st mrkpos struct
	@ const mrkpos & m2 : a const reference to the 2nd mrkpos struct
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	function used by the STL sort() and lower_bound() functions to compare 2 mrkpos values.
	-----------------------------------------------------------------------------
*/

bool stl_comp_mrkpos(const mrkpos & m1, const mrkpos & m2) {
	return m1.pos < m2.pos;
}


void fill_qtl(s_params* params, s_mk & marker, char** const qp_fields, double qtlpos){
	
	strcpy(marker.name, qp_fields[0]);	// qtl name
	marker.chr = atoi(qp_fields[1]);	// Index of the chromosome
	marker.num_qtl = params->nb_qtl;	// Index of the QTL
	marker.pos = qtlpos;				// Position of the marker/QTL
	
	// We split each favorable allele
	char **line2 = str_split(qp_fields[5], "/");
	marker.nb_all_fav = count(line2);

	// Memory allocation for the array of favorable alleles for each marker/QTL
	marker.all_fav = NULL;
	malloc2d(marker.all_fav, marker.nb_all_fav, params->geno_loci_max, char);

	// Loop for each favorable allele
	for(int j=0; j<marker.nb_all_fav; j++){
		// The favorable alleles names are copied
		if( (int)strlen(line2[j]) >= params->geno_loci_max ){
			fprintf(stderr, "Error in qtlpos file, at qtl %s:\nthe favorable allele symbol must not exceed %d character(s)!\n", qp_fields[0], params->geno_loci_max - 1);
			fflush(stderr);
			exit(1);
		}
		strcpy(marker.all_fav[j], line2[j]);
	}
	free(line2);
}


void fill_mk(s_params* params, s_mk & marker, char** const qp_fields, const mrkpos & mp){
	
	strcpy(marker.name, mp.mrk);		// Marker name
	marker.chr = atoi(qp_fields[1]);	// Index of the chromosome
	marker.num_qtl = params->nb_qtl;	// Index of the QTL
	marker.pos = mp.pos;				// Position of the marker/QTL

	marker.nb_all_fav = 0;
	marker.all_fav = NULL;
}


s_mk* assign_mrks_to_qtlw(s_params* params, s_mk* markers, const vector<mrkpos>& v_mrkpos, char** const qp_fields, char** const qw_fields){

	char* err = NULL;
	double qtlpos = strtod(qp_fields[2], &err);
	double posmin = qtlpos - atof(qw_fields[1]);
	double posmax = qtlpos + atof(qw_fields[1]);
	int nmrk = 0; // effective number of mrk assigned to the qtl

	if(*err != '\0') {
		fprintf(stderr, "Error: in the qtlpos file, : qtl %s position is incorrect.\n", qp_fields[0]);
		fflush(stderr);
		exit(1);
	}
	params->nb_qtl++;

	int ilast=-1, i;
	bool no_mrk_in_w = true;

	for(i=0; i<(int)v_mrkpos.size(); i++){
		if(v_mrkpos[i].pos < qtlpos) ilast=i;
		if (posmin <= v_mrkpos[i].pos && v_mrkpos[i].pos <= posmax){
			no_mrk_in_w = false;
			
			if(i-ilast==1){ // Insert the QTL
				realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
				fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

				realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
				params->tab_index_qtl[params->nb_mks_qtl] = 1;
				
				params->nb_mks_qtl++;
			}
			// insert the MARKER
			realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
			fill_mk(params, markers[params->nb_mks_qtl], qp_fields, v_mrkpos[i]);

			realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
			params->tab_index_qtl[params->nb_mks_qtl] = 0;

			params->nb_mks_qtl++;
			params->nb_mks++;
			nmrk++;
		}
	}
	if (no_mrk_in_w){ // there is no mrk within the window.
		fprintf(stderr, "Error : no mrk within the window could be selected for the qtl %s!\nA larger window could be necessary.", qp_fields[0]);
		fflush(stderr);
		exit(1);
	}
	else if(i-ilast==1){ // Insert the QTL, case where every mrk pos < qtlpos
		realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
		fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

		realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
		params->tab_index_qtl[params->nb_mks_qtl] = 1;

		params->nb_mks_qtl++;
	}
	return markers; // return the reallocated address
}

s_mk* assign_mrks_to_qtln(s_params* params, s_mk* markers, const vector<mrkpos>& v_mrkpos, char** const qp_fields, char** const qn_fields){

	double qtlpos = atof(qp_fields[2]);
	int mrknb = atoi(qn_fields[1]); // wanted nb of flanking mrks
	int nmrk = 0; // effective number of mrk assigned to the qtl
	size_t vsize = v_mrkpos.size(); // nb of mrks in the map belonging to the current chr

	params->nb_qtl++;

	mrkpos s_qtlpos; s_qtlpos.pos = qtlpos; 
	
	vector<mrkpos>::const_iterator beg = v_mrkpos.begin();
	// return an iterator to the first element (mrk pos) >= qtlpos, or v_mrkpos.end() if all mrkpos < qtlpos.
	vector<mrkpos>::const_iterator ge = lower_bound(beg, v_mrkpos.end(), s_qtlpos, stl_comp_mrkpos);

	int ige = ge - beg;

	if(mrknb>(int)vsize) mrknb=vsize; // warning, nb of requested mrks is greater than the nb of mrks available on this chr.

	set<int> closest_mrk_indexes; // stores the indexes of the n nearest mrk to the qtl (sorted in ascending order)
	set<int>::iterator it;

	double deltal, deltar;
	int il=ige, ir=ige, iclosest;
	il--;

	for(int i=0; i<mrknb; i++){

		if(il<0) { // left limit reached => the remaining closest mrk are on the right
			closest_mrk_indexes.insert(ir++); continue;
		}
		if(ir==(int)vsize) { // right limit reached => the remaining closest mrk are on the left
			closest_mrk_indexes.insert(il--); continue;
		}
		deltal = qtlpos - v_mrkpos[il].pos; //delta left
		deltar = v_mrkpos[ir].pos - qtlpos; //delta right

		if(deltal<deltar) iclosest = il--;	// the closest is the left one 
		else iclosest = ir++;				// the closest is the right one
		
		closest_mrk_indexes.insert(iclosest); // append the index of the current closest.
	}
	for(it=closest_mrk_indexes.begin(); it!=closest_mrk_indexes.end(); ++it){

		if(*it == ige){ // insert the QTL
			realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
			fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

			realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
			params->tab_index_qtl[params->nb_mks_qtl] = 1;
			
			params->nb_mks_qtl++;
		}
		// insert the MARKER
		realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
		fill_mk(params, markers[params->nb_mks_qtl], qp_fields, v_mrkpos[*it]);

		realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
		params->tab_index_qtl[params->nb_mks_qtl] = 0;

		params->nb_mks_qtl++;
		params->nb_mks++;
		nmrk++;
	}
	// Insert the QTL, case where every mrk pos < qtlpos
	if(ige == (int)vsize){
		realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
		fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

		realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
		params->tab_index_qtl[params->nb_mks_qtl] = 1;

		params->nb_mks_qtl++;
	}
	return markers; // return the reallocated address
}

s_mk* assign_mrks_to_qtll(s_params* params, s_mk* markers, const vector<mrkpos>& v_mrkpos, char** const qp_fields, char** const ql_fields){

	double qtlpos = atof(qp_fields[2]);
	int mrknb = count(ql_fields);

	set<string> set_mrk_liste(ql_fields+1, ql_fields+mrknb); // mrks are lexicographically ordered.
	set<string>::iterator it;
	set_mrk_liste.erase(""); // remove potential empty fields

	params->nb_qtl++;
	
	int ilast=-1, i;
	for(i=0; i<(int)v_mrkpos.size(); i++){ // iterate on markers (sorted in ascendant order)
		if(v_mrkpos[i].pos < qtlpos) ilast=i; // memorize the last index when mrk pos < qtlpos
		if ((it=set_mrk_liste.find(v_mrkpos[i].mrk)) != set_mrk_liste.end()){
			
			if(i-ilast==1) { // insert the QTL here (previous mrk pos < qtlpos, and present mrk pos >= qtl pos)
				realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
				fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

				realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
				params->tab_index_qtl[params->nb_mks_qtl] = 1;
				
				params->nb_mks_qtl++;
			}
			// insert the MARKER
			realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
			fill_mk(params, markers[params->nb_mks_qtl], qp_fields, v_mrkpos[i]);

			realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
			params->tab_index_qtl[params->nb_mks_qtl] = 0;

			params->nb_mks_qtl++;
			params->nb_mks++;
			set_mrk_liste.erase(it);
		}
	}
	// at least one marker is not in the map !
	if(!set_mrk_liste.empty()){
		fprintf(stderr, "Error at qtl %s, following mrks are not in the map file:\n",qp_fields[0]);
		for(it=set_mrk_liste.begin(); it!=set_mrk_liste.end(); ++it) fprintf(stderr, " %s", (*it).c_str());
		fflush(stderr);
		exit(1);
	}
	// insert the QTL, case where every mrk pos < qtlpos
	if(i-ilast==1){
		realloc1d(markers, params->nb_mks_qtl + 1, s_mk);
		fill_qtl(params, markers[params->nb_mks_qtl], qp_fields, qtlpos);

		realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);
		params->tab_index_qtl[params->nb_mks_qtl] = 1;

		params->nb_mks_qtl++;
	}
	return markers; // return the reallocated address
}

/**-----------------------------------------------------------------------------
	## FUNCTION:
	s_mk *init_tab_mks(s_params *params)
	-----------------------------------------------------------------------------
	## RETURN:
	@ s_mk *markers : array of markers (QTL included)
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_params *params : structure containing all data parameters
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Initializing an array of markers (parsing data).
	-----------------------------------------------------------------------------
*/

s_mk *init_tab_mks1(s_params *params){

	s_mk *markers = NULL;

	int cpt_lines = 0, nb_columns_map = 5;
	char *temp = (char *)malloc(NB_LINE_MEAN*sizeof(char));
	int last_chr = 0, last_qtl = 0;
	double last_pos = 0.0;

	// Open the marker file with the reading mode
	FILE *file = open_file(params->file_mks, "r") ;

	// We read the input file (markers file) line by line
	while(fgets(temp, NB_LINE_MEAN, file) && (strcmp(temp, "\n") != 0) && (strcmp(temp, "\r\n") != 0) ){

		search(temp, '\n'); // New line at the end of the string is removed
		search(temp, '\r');

		//We split each column of each line in an array
		char **line = str_split(temp, "\t");

		// The number of columns for each line is checked
		if( count(line) != nb_columns_map && count(line) != nb_columns_map - 1 ){
			fprintf(stderr, "Error: number of information l.%d of the map file\n", cpt_lines + 1);
			fflush(stderr);
			exit(1);
		}

		// The first line is checked [geno/chr/qtl/pos/all+]
		if( cpt_lines == 0 ){
			if( (strncmp(line[0], "Locus", 5) != 0) || (strncmp(line[1], "Chr", 3) != 0) || (strncmp(line[2], "QTL", 3) != 0) || (strncmp(line[3], "Pos", 3) != 0) || (strncmp(line[4], "All+", 4) != 0) ){
				fprintf(stderr, "Error: the header of the map file is not conform!\n");
				fflush(stderr);
				exit(1);
			}
		}
		else{
			// Memory allocation for each marker
			realloc1d(markers, params->nb_mks_qtl + 1, s_mk);

			// We take information for each marker
			strcpy(markers[params->nb_mks_qtl].name, line[0]); // Marker name
			markers[params->nb_mks_qtl].chr = atoi(line[1]); // Index of the chromosome
			markers[params->nb_mks_qtl].num_qtl = atoi(line[2]); // Index of the QTL
			markers[params->nb_mks_qtl].pos = atof(line[3]); // Position of the marker/QTL

			/*
			if( markers[params->nb_mks_qtl].chr < last_chr ){
				fprintf(stderr, "Error: line %d of the map file:\nChromosomes indexes have to be in ascendant order\n", cpt_lines + 1);
				fflush(stderr); exit(1);
			}
			*/
			if( markers[params->nb_mks_qtl].num_qtl < last_qtl ){
				fprintf(stderr, "Error: line %d of the map file:\nQTLs indexes have to be in ascendant order\n", cpt_lines + 1);
				fflush(stderr); exit(1);
			}
			if( markers[params->nb_mks_qtl].num_qtl == last_qtl && markers[params->nb_mks_qtl].pos < last_pos ){
				fprintf(stderr, "Error: line %d of the map file:\nThe QTL position and its Flanking Markers positions, have to be in ascendant order\n", cpt_lines + 1);
				fflush(stderr); exit(1);
			}
			last_chr = markers[params->nb_mks_qtl].chr;
			last_qtl = markers[params->nb_mks_qtl].num_qtl;
			last_pos = markers[params->nb_mks_qtl].pos;

			// Memory allocation for each marker
			realloc1d(params->tab_index_qtl, params->nb_mks_qtl + 1, int);

			// QTL position: keep its index
			if( strncmp(markers[params->nb_mks_qtl].name, "qtl", 3) == 0 ){

				// The number of columns for each QTL line is checked
				if( count(line) != nb_columns_map ){
					fprintf(stderr, "Error: number of information for the QTL %s l.%d of the map file\n", line[0], cpt_lines + 1);
					fflush(stderr);
					exit(1);
				}

				// We split each favorable allele
				char **line2 = str_split(line[4], "/");
				markers[params->nb_mks_qtl].nb_all_fav = count(line2);

				// Memory allocation for the array of favorable alleles for each marker/QTL
				markers[params->nb_mks_qtl].all_fav = NULL;
				malloc2d(markers[params->nb_mks_qtl].all_fav, markers[params->nb_mks_qtl].nb_all_fav, params->geno_loci_max, char);

				// Loop for each favorable allele
				for(int j=0; j<markers[params->nb_mks_qtl].nb_all_fav; j++){
					// The favorable alleles names are copied
					if( (int)strlen(line2[j]) >= params->geno_loci_max ){
						fprintf(stderr, "Error l.%d of the map file: the favorable allele name must be only one character (i.e a, b, c...) !\n", cpt_lines + 1);
						fflush(stderr);
						exit(1);
					}
					strcpy(markers[params->nb_mks_qtl].all_fav[j], line2[j]);
				}
				free(line2);

				params->nb_qtl++;
				params->vect_qtl_names.push_back(markers[params->nb_mks_qtl].name);
				params->tab_index_qtl[params->nb_mks_qtl] = 1;
			}
			// Marker position
			else{
				markers[params->nb_mks_qtl].nb_all_fav = 0;
				markers[params->nb_mks_qtl].all_fav = NULL;
				params->nb_mks++;
				params->tab_index_qtl[params->nb_mks_qtl] = 0;
			}
			params->nb_mks_qtl++;
		}
		free(line);
		cpt_lines++;
	}

	// Check if the QTL selected is not out of range
	if( (params->qtl_selected != 0) && (params->qtl_selected > params->nb_qtl) ){
		fprintf(stderr, "Error: the QTL selected for the analysis is out of range regarding the map file !\n");
		fflush(stderr);
		exit(1);
	}

	if( fclose(file) == EOF ){
		fprintf(stderr, "Error closing file %s !\n", params->file_mks);
		fflush(stderr);
		exit(1);
	}
	else {
		//printf("The file %s was read correctly...\n",params->file_mks); // RECAP NB MKS / QTL IN THE INTERFACE
	}
	free(temp);

	return markers;
}

/**-----------------------------------------------------------------------------
	## FUNCTION:
	s_mk *init_tab_mks2(s_params *params, const char* file_qp)
	-----------------------------------------------------------------------------
	## RETURN:
	@ s_mk *markers : array of markers (QTL included)
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_params *params    : structure containing all data parameters
	@ const char* file_qp : name of the main QTL input file (.qtlpos)
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Initializing an array of markers (parsing data).
	-----------------------------------------------------------------------------
*/

s_mk *init_tab_mks2(s_params *params, const char* file_qp){

	int cpt_qtl = 0, cpt_mrk = 0, nb_columns_qtlpos = 6, nb_columns_mrk = 3;
	s_mk *markers = NULL;
		
	char file_qa[NB_PATH_MAX]; // annex QTL input file
	strcpy(file_qa, file_qp);
	
	char header_qa[NB_LINE_MEAN];
	char method = 0;
	char* ext = strrchr(file_qa, '.');
	*ext = '\0';

	// test the existence of the annex QTL file (priority: qtll > qtlw > qtln)
	if (access(strcat(file_qa, ".qtll"), F_OK) == 0){ // file_qtll exist
		method = 'l';
		strcpy(header_qa, QTLL_HEADER);
	}	
	else if (!(*ext='\0') && access(strcat(file_qa, ".qtlw"), F_OK) == 0){ // file_qtlw exists
		method = 'w';
		strcpy(header_qa, QTLW_HEADER);
	}	
	else if (!(*ext='\0') && access(strcat(file_qa, ".qtln"), F_OK) == 0){ // file_qtln exists
		method = 'n';
		strcpy(header_qa, QTLN_HEADER);
	}
	else { // none of these files exists !
		fprintf(stderr, "Error : annex QTL input file (either .qtll or .qtlw or .qtln) is missing\n");
		fflush(stderr);
		exit(1);
		//// TODO : automatic marker selection (TAGSNP algorithm)
	}
	map< string, vector<mrkpos> > chr_mrkpos;
	map< string, vector<mrkpos> >::iterator it;
	char line[NB_LINE_MEAN];
	char qa_line[NB_LINE_MEAN];

	// open map file in readonly mode (default)
	ifstream fmap(params->file_mks);
	if(!fmap) {fprintf(stderr,"Error opening %s.\n",params->file_mks); exit(1);}

	// map file header
	fmap.getline(line,NB_LINE_MEAN);
	chomp(line);
	if( strcmp(line,"mrk\tchr\tpos") != 0 ){
		fprintf(stderr, "Error: the header of the map file is not conform!\n");
		fflush(stderr);
		exit(1);
	}
	// map file data are stored in the hash map chr_mrkpos : chr_i => [ struct{mrk1_name, mrk1_pos}, struct{mrk2_name, mrk2_pos}, ... ]
	while(!fmap.getline(line,NB_LINE_MEAN).eof()){

		chomp(line); if(*line == '\0') continue; // ignore empty line
		char** map_fields = str_split(line,"\t");

		if( count(map_fields) != nb_columns_mrk ){
			fprintf(stderr, "Error: in the map file, line %d : incorrect number of columns\n", cpt_mrk+2);
			fflush(stderr);
			exit(1);
		}
		char* err = NULL;
		mrkpos mp; // struct containing the mrk name and position.
		strcpy(mp.mrk, map_fields[0]);
		mp.pos = strtod(map_fields[2], &err);
		if(*err != '\0') {
			fprintf(stderr, "Error: in the map file, line %d: incorrect mrk position.\n", cpt_mrk+2);
			fflush(stderr);
			exit(1);
		}
		chr_mrkpos[map_fields[1]].push_back(mp);
		free(map_fields);
		//cpt_mrk++;
	}
	fmap.close();
	
	for(it = chr_mrkpos.begin(); it != chr_mrkpos.end(); ++it){
		vector<mrkpos>& v_mrkpos = it->second;
		// markers are sorted (ascendent) according to their positions (see stl_comp_mrkpos definition).
		sort(v_mrkpos.begin(), v_mrkpos.end(), stl_comp_mrkpos);
	}
	ifstream fqp(file_qp); // one main qtl positions file with qtlpos extention.
	if(!fqp) {fprintf(stderr,"Error opening %s.\n",file_qp); exit(1);}
	
	ifstream fqa(file_qa); // one anexe qtl file among 3 possible formats, with .qtlw or .qtln or qtll extention.
	if(!fqa) {fprintf(stderr,"Error opening %s.\n",file_qa); exit(1);}

	// qtlpos file header
	fqp.getline(line,NB_LINE_MEAN);
	chomp(line);
	if( strcmp(line,"QTL\tchr\tpos\tCI min\tCI max\tall+") != 0 ){
		fprintf(stderr, "Error: the header of the qtlpos file is not conform!\n");
		fflush(stderr);
		exit(1);
	}
	// annex qtl file header
	fqa.getline(qa_line,NB_LINE_MEAN);
	chomp(qa_line);

	if( strcmp(qa_line,header_qa) != 0 ){
		fprintf(stderr, "Error: the header of the qtl%c file is not conform!\n", method);
		fflush(stderr);
		exit(1);
	}

	// parsing both qtlpos and annex qtl files simultaneously !
	while(!fqp.getline(line,NB_LINE_MEAN).eof() && !fqa.getline(qa_line,NB_LINE_MEAN).eof()){

		chomp(line); if(*line == '\0') continue; // ignore empty line

		char** qp_fields = str_split(line,"\t");
		if( count(qp_fields) != nb_columns_qtlpos ){
			fprintf(stderr, "Error: in the qtlpos file, line %d : incorrect number of columns.\n", cpt_qtl+2);
			fflush(stderr);
			exit(1);
		}
		char** qa_fields = str_split(qa_line,"\t");

		if(strcmp(qp_fields[0], qa_fields[0]) != 0) {
			fprintf(stderr, "Error : in qtlpos and qtl%c files, line %d: qtls names does not match.\n", method, cpt_qtl+2);
			exit(1);
		}
		// v_mrkpos contains all the mrks from the chr this qtl belongs to.
		const vector<mrkpos>& v_mrkpos = chr_mrkpos[qp_fields[1]];
		if(v_mrkpos.empty()){
			fprintf(stderr,"Error : the chr %s is not present in the map file.\n", qp_fields[1]);
			exit(1);
		}
		// 3 different methods : w (flanking window in cM), n (nb of flanking closest mrk), l (list of mk).
		if(method == 'w')
			markers = assign_mrks_to_qtlw(params, markers, v_mrkpos, qp_fields, qa_fields);
		else if(method == 'n')
			markers = assign_mrks_to_qtln(params, markers, v_mrkpos, qp_fields, qa_fields);
		else if(method == 'l')
			markers = assign_mrks_to_qtll(params, markers, v_mrkpos, qp_fields, qa_fields);
		
		params->vect_qtl_names.push_back(qp_fields[0]);
		cpt_qtl++;
		
		free(qp_fields);
		free(qa_fields);
	}
	//// test if fqa is also eof ?

	fqp.close();
	fqa.close();

	// Check if the QTL selected is not out of range
	if( (params->qtl_selected != 0) && (params->qtl_selected > params->nb_qtl) ){
		fprintf(stderr, "Error: the QTL selected for the analysis is out of range regarding the map file !\n");
		fflush(stderr);
		exit(1);
	}
	return markers;
}


s_mk *init_tab_mks(s_params *params){

	char file_qtlpos[NB_PATH_MAX]; // QTLs input file (with qtlpos extension)
	strcpy(file_qtlpos, params->file_mks);
	
	char* ext = strrchr(file_qtlpos, '.'); *ext = '\0';
	strcat(file_qtlpos, ".qtlpos");

	// test the existence of file_qtlpos
	int ret = access(file_qtlpos, F_OK);
	
	if (ret == -1 && errno == ENOENT) // file_qtlpos does not exist
		return init_tab_mks1(params); // => assume the .map conforms the old format

	params->qtlpos = true; //// attribut provisoire ? (voir dump_tab_mks())	
	// file_qtlpos exists => new format, 3 files (map + qtlpos + one among qtlw, qtln and qtll)
	return init_tab_mks2(params, file_qtlpos);
}

/**-----------------------------------------------------------------------------
	## FUNCTION:
	free_markers(s_mk *markers, int nb_mks)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_mk *markers : array of markers (QTL included)
	@ int nb_mks_qtl : total number of markers (QTL included)
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Freeing memory for the names of favorable alleles for each marker.
	-----------------------------------------------------------------------------
*/

void free_markers(s_mk *markers, int nb_mks_qtl){

	for(int i=0; i<nb_mks_qtl; i++){
		free_alloc2d(markers[i].all_fav);
	}
	free(markers);
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void print_tab_mks(s_mk *markers, int nb_mks)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_mk *markers : array of markers (QTL included)
	@ int nb_mks_qtl : total number of markers (QTL included)
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Print table markers.
	-----------------------------------------------------------------------------
*/

void print_tab_mks(s_mk *markers, int nb_mks_qtl){
	int i, j;

	printf("AFFICHAGE DU TABLEAU DE MARQUEURS:\n");

	for(i=0; i<nb_mks_qtl; i++){
		printf("%s | K: %d | n: %d | pos: %1.1f | nb_a_f: %d | all_fav: ", markers[i].name, markers[i].chr, markers[i].num_qtl, markers[i].pos, markers[i].nb_all_fav);
		for(j=0; j<markers[i].nb_all_fav; j++){
			printf("%s ", markers[i].all_fav[j]);
		}
		printf("|\n");
	}
	printf("\n");
}

/**-----------------------------------------------------------------------------
	## FUNCTION:
	void print_tab_mks(s_mk *markers, int nb_mks)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_mk *markers : array of markers (QTL included)
	@ int nb_mks_qtl : total number of markers (QTL included)
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Print table markers.
	-----------------------------------------------------------------------------
*/

void dump_tab_mks(s_mk *markers, s_params *params){
	int i, j;

	char file_qtlmks[NB_PATH_MAX];		
	char* sep = strrchr(params->file_mks, '/');

	// file is in the current directory.
	if (sep == NULL) sep=params->file_mks;
	else sep = sep+1;

	sprintf(file_qtlmks, "%s/%s\n", params->file_root, sep);

	//strcpy(file_qtlmks, params->file_mks);

	char* ext = strrchr(file_qtlmks, '.');
	strcpy(ext, ".mrkqtl.map");
		
	FILE* f = open_file(file_qtlmks, "w");

	fprintf(f,"%s\n",MRK_QTL_HEADER);

	for(i=0; i < params->nb_mks_qtl; i++){
		fprintf(f, "%s\t%d\t%d\t%lf", markers[i].name, markers[i].chr, markers[i].num_qtl, markers[i].pos);
		int nb_all_fav = markers[i].nb_all_fav;
		if(nb_all_fav > 0) {
			fprintf(f,"\t"); 
			for(j=0; j< nb_all_fav; j++){
				fprintf(f,"%s", markers[i].all_fav[j]);
				if (j+1 < nb_all_fav) fprintf(f,"/");
			}
		}
		fprintf(f,"\n");
	}
	fclose(f);
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void find_index(s_mk *markers, s_params *params, int q)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_mk *markers : table of markers
	@ s_params *params : structure containing all data parameters
	@ int q : current QTL
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Assignment of the lower and upper marker flanking the current QTL.
	-----------------------------------------------------------------------------
*/

void find_index(s_mk *markers, s_params *params, int q){

	/* We find the index for the lower and upper markers */
	for(int m=0; m<(params->nb_mks_qtl); m++){
		/* We check if the marker QTL index is the same as the current QTL */
		if( markers[m].num_qtl == q ){
			/* We find the lower marker corresponding to that QTL */
			if( params->nb_locus == 0 ){
				params->num_mk_inf = m;
			}
			params->nb_locus++;
		}
	}
	/* Index of the upper marker at the current QTL position */
	params->num_mk_sup = params->nb_locus + params->num_mk_inf - 1;
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void compute_rec_frac(s_mk *markers, s_params *params)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_mk *markers : table of markers
	@ s_params *params : structure containing all data parameters
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Compute the recombination fraction r (and R) for each possible interval
	between num_mk_inf and num_mk_sup.
	Store the results in vect_rec_frac (attribute of params).
	-----------------------------------------------------------------------------
*/

void compute_rec_frac(s_mk* markers, s_params* params){
	
	double delta_pos;
	
	s_rec_frac rf;
	params->vect_rec_frac.clear();
	params->vect_rec_frac.resize( params->nb_locus-1, vector<s_rec_frac>(params->nb_locus) );
	
	// j > i
	for(int i=0; i<params->nb_locus; i++){
		for(int j=i+1; j<params->nb_locus; j++){
			
			delta_pos = markers[params->num_mk_inf + j].pos - markers[params->num_mk_inf + i].pos;

			rf.r = 0.5 * (1.0 - exp(- 2.0 * delta_pos / 100));	// r
			rf.R = (2.0 * rf.r) / (1.0 + (2.0 * rf.r));			// R

			params->vect_rec_frac[i][j] = rf;
		}
	}
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	int is_fav(char *allele, s_mk *markers, int index_mk)
	-----------------------------------------------------------------------------
	## RETURN:
	[1] if the tested allele is favorable, [0] otherwise.
	-----------------------------------------------------------------------------
	## PARAMETRES:
	@ char *allele : allele searched
	@ s_mk *markers : array of markers
	@ int index_mk : index of the current marker
	-----------------------------------------------------------------------------
	## SPECIFICATION:
	Evaluates the presence or absence of a favorable allele at a QTL.
	-----------------------------------------------------------------------------
*/

int is_fav(const char *allele, s_mk *markers, int index_mk){

	// We read the table of string containing the list of favorable alleles
	for(int i=0; i<markers[index_mk].nb_all_fav; i++){
		if( strcmp(markers[index_mk].all_fav[i], allele) == 0 ){
			return 1;
		}
	}
	return 0; // The allele was not found in the list
}


/**-----------------------------------------------------------------------------
	## FUNCTION:
	void check_fav_alleles_input_files(s_indiv *tab_indiv, s_mk *markers, s_params *params)
	-----------------------------------------------------------------------------
	## RETURN:
	-----------------------------------------------------------------------------
	## PARAMETERS:
	@ s_indiv *tab_indiv : array of individuals
	@ s_mk *markers : array of all markers
	@ s_params *params : structure containing all data parameters
	-----------------------------------------------------------------------------
	## SPECIFICATION: check the correspondence of the favorable alleles present
	in the 2 input files after table initialization.
	-----------------------------------------------------------------------------
*/

void check_fav_alleles_input_files(s_mk *markers, s_params *params){

	int qtl = 1;

	// Loop on all the markers
	for(int m = 0; m < params->nb_mks_qtl; m++){
		// We are at a QTL position
		if( params->tab_index_qtl[m] == 1 ){
			// We read the table of string containing the list of favorable alleles
			for(int i=0; i<markers[m].nb_all_fav; i++){
				if( params->str2coded.find(markers[m].all_fav[i]) == params->str2coded.end() ){
					fprintf(stderr, "Error: favorable alleles at QTL %d mentioned in the .map file don't match the parental alleles of the genotypes/pedigre file (.dat)!\n", qtl);
					exit(1);
				}
			}
			qtl++;
		}
	}
	// The different favorable alleles correspond in both files
}
