#include <string>
#include <sstream>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#include <math.h>
#include <stdio.h>

using namespace std;


// Interface to gsl_cdf_hypergeometric_Q with some add-on behaviour:
// gsl_cdf_hypergeometric_Q returns P(x>k) this interface returns P(x>=k)
double hypgeom_Q(unsigned int N, unsigned int K, unsigned int n, unsigned int k){
	if(k==0)		//special case P(x>=k|k=0)= 1
		return 1;
	else			// normal case, should give P(x>=k)
		return gsl_cdf_hypergeometric_Q(k-1,K,N-K,n);
}

// Generates random permutation of length N, containing K true-values and N-K false-values
bool* get_random_genome(unsigned int N, unsigned int K, const gsl_rng* rng){
	bool* perm= new bool[N];
	unsigned int i;
	for(i=0; i<K; i++)
		perm[i]= true;
	for(i=K; i<N; i++)
		perm[i]= false;
	int r;
	bool temp;
	for(unsigned int i=0; i<N-1; i++){

		r= (int)floor(gsl_rng_uniform(rng) * (N-i));
		temp= perm[r];
		perm[r]= perm[i];
		perm[i]= temp;
	}
	return perm;
}

/* Scans binary array for clusters of consecutive true-values and returns an array structure
 specifying the segments found.
 parameters:
	1. reference to binary array to be scanned
	2. maximum gap length
 returns:
	N x 2 matrix C specified as an array of array pointers, where each subarray @{$C[i]} is
	a pair of numbers defining the starting, $C[i]->[0], and ending, $C[i]->[1], indices of
	segment i.
*/
vector<int*> get_segments(bool* b, int N, int maxgap){
	vector<int*> segms;
	int* segm;
	int i1=0;
	int i2=0;
	int i3=0;
	while(i1<N){ // loop1
		//print "$i1\n";
		if(!b[i1]){
			i1++;
		}
		else{
			i2=i1;
			while(i2<N){ // loop2
				//print "   $i2\n";
				if(b[i2]){
					if(i2==(N-1)){
						segm= new int[2];
						segm[0]= i1; segm[1]= i2;
						segms.push_back(segm);
						i1=N; // iterated till the end of array -> breaking loop2 and loop1
						break;
					}
					else{
						i2++;
					}
				}
				else{
					i3=i2;
					while(i3<N){//loop3
						// print "      $i3\n";
						if((i3-i2)>maxgap){
							segm= new int[2];
							segm[0]= i1; segm[1]= i2-1;
							segms.push_back(segm);
							i1=i3;
							i2=N;  // back to loop1 -> breaking loop3 and loop2
							break;
						}
						else if(b[i3]){
							i2=i3;
							break; // back to loop2 -> beaking loop3
						}
						else if(i3==N-1){
							segm= new int[2];
							segm[0]=i1; segm[1]= i2-1;
							segms.push_back(segm);
							i2=N; // iterated till the end of array -> breaking loop3, loop2 and loop1
							i1=N;
							break;
						}
						else{
							i3++;
						}
					} //close loop3
				}
			} //close loop2
		}
	} //close loop1
	return segms;
}


vector<unsigned int*> get_nk_values(bool* b, vector<int*> segms){
	vector<unsigned int*> nks;
	unsigned int* nk;
	unsigned int n,k;
	for(unsigned int i=0; i<segms.size(); i++){
		k=0;
		for(int j=segms[i][0]; j<=segms[i][1]; j++)
			if(b[j]) k++;
		n = segms[i][1]-segms[i][0]+1;

		nk= new unsigned int[2];
		nk[0]= n;
		nk[1]= k;
		nks.push_back(nk);
	}
	return nks;
}


vector<double> bootstrap(int N,
						 int K,
						 vector<double> Ps,
						 int I,
						 int maxgap,
						 unsigned long seed){

	gsl_rng * rng = gsl_rng_alloc (gsl_rng_taus);
	gsl_rng_set (rng, seed);

	bool* genome;
	vector<int*> segms;
	vector<int*> segms_temp;
	vector<unsigned int*> nk; // (n,k) value pairs obtained from segments
	vector<double> p_sim;
	vector<double> p_sim_min(I,0);

	for(int i=0; i<I; i++){
		genome= get_random_genome(N,K,rng);
		segms.clear();
		segms_temp.clear();
		nk.clear();
		p_sim.clear();

		for(int j=0; j<=maxgap; j++){
			segms_temp= get_segments(genome,N,j);
			segms.insert(segms.end(),segms_temp.begin(),segms_temp.end());
		}
		segms= get_segments(genome,N,maxgap);
		nk= get_nk_values(genome,segms);

		for(unsigned int j=0; j< nk.size(); j++){
			p_sim.push_back( hypgeom_Q(N,K,nk[j][0],nk[j][1]) );
		}
		if(p_sim.size() == 0)
			p_sim_min[i] = 0;
		else
			p_sim_min[i]= *min_element( p_sim.begin(), p_sim.end() );

	}

	sort(p_sim_min.begin(), p_sim_min.end());

	vector<double> Ps_adj(Ps.size(), 0 );

	for(unsigned int i=0; i<Ps.size(); i++){
		//int smaller_then= int( lower_bound(p_sim.begin(), p_sim.end(), p_val.at(i), >) );
		int j= 0;
		double p= Ps[i];
		while(j<I && p_sim_min[j]<=p){j++;}
		Ps_adj[i]= double(j)/I;
	}


	return Ps_adj;
}



/*************************************************************/
/*
/* TESTS
/*
/*************************************************************/

/* Tests generation of random boolean vectors (in this context genomes)
/*/
void test_get_random_genome(int N, int K, int I){
	gsl_rng * rng = gsl_rng_alloc (gsl_rng_taus);
	gsl_rng_set (rng, time(NULL));

	cout << "# Testing generation of randome genomes, with N=" << N << " & K=" << K << "\n";
	cout << "Iteration\tGenome\n";
	for(int i=0; i<I; i++){
		cout << (i+1) << "\t";
		bool* gen= get_random_genome(N,K,rng);
		for(int j=0; j<N; j++){
			cout << " " << gen[j];
		}
		cout << "\n";
	}
}

/* Tests detection of clusters from boolean vectors (get_clusters) &
/* calculation of cluster length and k-values.
/*/
void test_get_clusters(int N, int K, int k_th, int I, int maxgap){
	bool* gen;
	vector<int*> segms;
	vector<unsigned int*> nks;
	gsl_rng * rng = gsl_rng_alloc (gsl_rng_taus);
	gsl_rng_set (rng, time(NULL));

	cout << "# Testing clustering, with maxgap=" << maxgap << "\n";
	cout << "Iteration\tVector\tSegments\t(n,k)\n";
	cout << "\t";
	for(int i=0; i<N; i++){
	 	printf("%4u",i);
	}

	int count=0;
	bool counted=false;
	cout << "\n";
	for(int i=0; i<I; i++){
		gen= get_random_genome(N,K,rng);
		segms= get_segments(gen,N,maxgap);
		nks= get_nk_values(gen,segms);

		//printing
		cout << (i+1) << "\t";
		for(int j=0; j<N; j++)
			printf("%4u", gen[j]);
		cout << "\n";

		for(unsigned int j=0; j<segms.size(); j++)
			cout << " (" << segms[j][0] << "," << segms[j][1] << ")";

		cout << "\n";
		for(unsigned int j=0; j<nks.size(); j++)
			cout << " (" << nks[j][0] << "," << nks[j][1] << ")";
		counted=false;
		for(unsigned int j=0; j<nks.size(); j++)
			if(nks[j][1]>=k_th){
				if(!counted){
					count++;
					counted=true;
				}
			}

		cout << "\n";
	}
	cout << "# random genomes with at least one window with k>=" << k_th << "\n" << count << "/" << I << "\n";
}



int main (int argc, char** argv) {
	// TESTING
	//test_get_random_genome(20,10,10); 	// test_gen_random_genome(N,K,Iterations)
	//test_get_random_genome(20,0,5);
	//test_get_clusters(2751,76,5,100,5); 		//test_get_clusters(N,K,k_th,Iterations,maxgap)
	//test_get_clusters(20,10,10,3);


	//checking input
	if(argc < 3){
		cerr << "USAGE:\n"
			   << argv[0] << " filein fileout\n"
			   << "Where filein contains two blocks of data delimited by an empy line.\n"
			   << "First block should contain lines composed of a header a tab and a \n"
			   << "header value. Required headers are N,K and maxgap. Lines with required\n"
			   << "headers can be in any order, lines with other headers are ignored.\n"
			   << "The second block should start with a line of tab-delimited headers,\n"
			   << "including, but not limiting to, \'n\' and \'k\'. The corresponding\n"
			   << "unsigned integer values should be delimited by tabs and layed out\n"
			   << "starting from the second line of this block.\n";
		exit(1);
	}


	unsigned int N,K,maxgap;
	char* fin= argv[1];
	ifstream in (fin , ifstream::in);
	char buffer[256];
	string header;


	//Reading data: Block 1
	while(in.getline(buffer,256)){
		string line(buffer);
		if(line == "")
			break;	//end of Block 1

		istringstream iss (line,istringstream::in);
		iss >> header;
		if(header == "N"){
			iss >> N;
		}else if(header == "K"){
			iss >> K;
		}else if(header == "maxgap"){
			iss >> maxgap;
		}
	}


	//Reading data: Block 2: headers
	in.getline(buffer,256);
	string line(buffer);
	istringstream iss (line,istringstream::in);
	vector<int> map;
	unsigned int col_count=0;
	while(iss >> header){
		if(header == "n"){ map.push_back(0);}
		else if(header == "k"){ map.push_back(1);}
		else{ map.push_back(-1); }
		col_count++;
	}
	//Reading data: Block 2: data
	vector<unsigned int> ns;
	vector<unsigned int> ks;
	vector<unsigned int> data[]= {ns,ks};
	unsigned int counter=0;
	int col;
	unsigned int value;
	string rubbish;
	do{
		col= map[counter++%col_count];
		if(col>=0){
			in >> value;
			if(in)
				data[col].push_back(value);
		}
		else{
			in >> rubbish;
		}
	}while(in);
	in.close();
	ns= data[0];
	ks= data[1];

	// Calculating P-values and P_adj-values
	unsigned int I= 1000;
	unsigned long seed= time(NULL);
	vector<double> Ps(ns.size(),0);
	for(unsigned int i=0; i<ns.size(); i++){
		Ps[i]= hypgeom_Q(N,K,ns[i],ks[i]);
	}
	vector<double> Ps_adj;
	Ps_adj= bootstrap(N,K,Ps,I,maxgap,seed);


	// Writing output: Block 1
	char* fout= argv[2];
	ofstream out (fout);
	out << "N\t" << N << "\n";
	out << "K\t" << K << "\n";
	out << "maxgap\t" << maxgap << "\n\n";

	// Writing output: Block 2
	out << "n\tk\tP\tP_adj\n";
	for(unsigned int i=0; i<Ps.size(); i++){
		out << ns[i] << "\t" << ks[i] << "\t" << Ps[i] << "\t" << Ps_adj[i] << "\n";
	}
	//
	return 0;

}