#include <iostream>
#include <fstream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <omp.h>

#include "org.h"
#include "amino_acids.h"
#include "util.h"

using namespace std;

struct avg
{
  double avg;
  int index;
};

struct Mproteome
{
  vector<string> proteins;
  //vector<double> bit_scores;
  double bit_scores[20];
  vector<string> files;
  int ordinal;
  string org_name;
  double total;
};

struct Fproteome
{
  vector<string> files;
  int ordinal;
  string org_name;
  double total;
};

struct v
{
  int id1;
  int id2;
  vector<int> counts;
  string most_unlikely;
  int uscore;
};

void do_counts(Mproteome p1, Mproteome p2, int tag_length, vector<int>& array, int& max, int& min,string& most_unlikely, int& uscore,int SOS_CUTOFF);
avg calc_running_avgs(int array[], int max, int min);
bool check_composition(string s,int tag_length,amino_acids& aas,int SOSCUTOFF);
int get_num_bins(string path,string bfn);
int calc_sos_cutoff(string path,string bfn);

int main(int argc, char* argv[])
{
  if(argc<2)
    {
      cout<<"Wrong use of function. ./dbc <full_path>"<<endl;
      exit(1);
    }
  
  string full_path=argv[1];
  util u;
  string bfn;
  string path;
  u.extract_paths(full_path,path,bfn);
  string FM = "M";
  int num_bins = get_num_bins(path,bfn);
  cout<<"num_bins = "<<num_bins<<endl;
  string PS = "P";
  int SOS_CUTOFF = calc_sos_cutoff(path,bfn);

  //string FM = argv[3]; //always set to true
  //  int num_bins=atoi(argv[4]); //set it from within
  //  string PS = argv[5]; //always set to P
  //int SOS_CUTOFF=atoi(argv[6]); //set it from within
  
  amino_acids aas;
  
  //1. read in info_file 
  ifstream instream;
  string info_file = path+"/"+bfn+"/"+bfn+"_info.txt";
  instream.open(info_file.c_str());
  if(instream.fail())
    {
      cout<<"Failed to open instream to "<<info_file<<endl; 
      exit(1);
    }

  //2. populate proteomes
  int num_orgs, num_orgs_ref, tag_length;
  instream>>num_orgs_ref>>num_orgs>>tag_length;
  
  num_orgs+=num_orgs_ref;

  //a. hold them all in memory
  //read in info file.
  if(FM=="M")
    {
      vector<Mproteome> proteomes;
      string eater,org_name,faa_path;
      int ordinal;
      instream>>ordinal>>org_name>>faa_path;
      while(!instream.eof())
	{
	  bool already_there=false;
	  for(int i=0; i<proteomes.size() && !already_there; i++)
	    {
	      if(ordinal==proteomes.at(i).ordinal && org_name==proteomes.at(i).org_name)
		{
		  proteomes.at(i).files.push_back(faa_path);
		  //cout<<"Already there. ordinal="<<ordinal<<" org_name="<<org_name<<endl;
		  //cout<<"pushing back "<<faa_path<<endl;
		  already_there=true;
		}
	    }
	  if(!already_there)
	    {
	      //cout<<"Not already there."<<endl;
	      //cout<<"ordinal = "<<ordinal<<" org_name = "<<org_name<<endl;
	      //cout<<"adding file "<<faa_path<<endl;
	      Mproteome mp_temp;
	      mp_temp.ordinal=ordinal;
	      mp_temp.org_name=org_name;
	      mp_temp.files.push_back(faa_path);
	      mp_temp.total=0.0;
	      proteomes.push_back(mp_temp);
	    }
	  instream>>ordinal>>org_name>>faa_path;
	}

      cout<<"Size of proteomes (M) = "<<proteomes.size()<<endl;
      
      //read in all the proteins.
      for(int i=0; i<proteomes.size(); i++)
	{
	  for(int j=0; j<proteomes.at(i).files.size(); j++)
	    {
	      ifstream instream;
	      instream.open( proteomes.at(i).files.at(j).c_str());
	      if(instream.fail())
		{
		  cout<<"Failed to open instream to "<<proteomes.at(i).files.at(j)<<endl;
		  exit(1);
		}

	      cout<<"Successfully opened "<<proteomes.at(i).files.at(j)<<endl;

	      string eater;
	      instream>>eater;
	      while(!instream.eof())
		{
		  if(eater[0]=='>')
		    {
		      char buffer[2000];
		      instream.getline(buffer,2000);
		      string info_line = eater+buffer;
		      
		      string protein;
		      instream>>eater;
		      while(eater[0]!='>' && !instream.eof())
			{
			  protein+=eater;
			  instream>>eater;
			}
		      proteomes.at(i).proteins.push_back(protein);	
		      //cout<<"info_line = "<<info_line<<endl<<protein<<endl;
		    }
		}
	      instream.close();
	      // cout<<"Size of "<<i<<" proteins = "<<proteomes.at(i).proteins.size()<<endl;
	    }

	  for(int j=0; j<20; j++)
	    {
	      proteomes.at(i).bit_scores[j]=0.0;
	    }

	  for(int j=0; j<proteomes.at(i).proteins.size(); j++)
	    {
	      for(int k=0; k<proteomes.at(i).proteins.at(j).length(); k++)
		{
		  int index=aas.get_aa(proteomes.at(i).proteins.at(j)[k]);
		  if(index>=0 && index<20)
		    {
		      proteomes.at(i).bit_scores[index]+=1.0;
		      proteomes.at(i).total+=1.0;
		    }
		}
	    }
	}
      
      ofstream outstream2;
      string outfile2 = path+"/"+bfn+"/bin_corrections_histograms.txt"; 
      outstream2.open(outfile2.c_str());
      if(outstream2.fail())
	{
	  cout<<"Failed to open outstream2 to "<<outfile2<<endl;
	  exit(1);
	}

      for(int i=0; i<proteomes.size(); i++)
	{
	  outstream2<<i<<" "<<proteomes.at(i).total<<endl;
	  //cout<<"Doing bitting for proteome i="<<i<<endl;
	  for(int j=0; j<20; j++)
	    {
	      //cout<<j<<" "<<proteomes.at(i).bit_scores[j]<<" ";
	      //outstream2<<aas.get_aa(j)<<" "<<proteomes.at(i).bit_scores[j]<<" ";
	      proteomes.at(i).bit_scores[j]/=proteomes.at(i).total; //now frequencies.	  
	      //cout<<proteomes.at(i).bit_scores[j]<<" ";
	      //outstream2<<proteomes.at(i).bit_scores[j]<<" ";
	      proteomes.at(i).bit_scores[j]=abs(log(proteomes.at(i).bit_scores[j]));
	      //outstream2<<proteomes.at(i).bit_scores[j]<<" ";
	      proteomes.at(i).bit_scores[j]/=2.0;
	      outstream2<<proteomes.at(i).bit_scores[j]<<endl;
	      //cout<<proteomes.at(i).bit_scores[j]<<endl;
	    }
	  outstream2<<endl;
	}
      
      ofstream outstream;
      string outfile = path+"/"+bfn+"/bin_corrections.txt";
      outstream.open(outfile.c_str());//, ios::out | ios::app);
      if(outstream.fail())
	{
	  cout<<"Failed to open outstream to "<<outfile<<endl;
	  exit(1);
	}

      ofstream outstream3;
      string outfile3 = path+"/"+bfn+"/bin_corrections_unlikelies.txt"; 
      outstream3.open(outfile3.c_str(), ios::out | ios::app);
      if(outstream3.fail())
	{
	  cout<<"Failed to open outstream3 to "<<outfile3<<endl;
	  exit(1);
	}

      outstream<<num_bins<<endl<<endl;

      //      int i_starter=32;
      //int j_starter=1687;
      //bool start_writing=false;

      if(PS=="P")
	{
	  cout<<"Using omp...presumably..."<<endl;

	  vector<v> stored_counts;

#pragma omp parallel for
	  for(int i=0; i<proteomes.size()-1; i++)
	    {
	      for(int j=i+1; j<proteomes.size(); j++)
		{
		  cout<<i<<" to "<<j<<endl;
		  
		  vector<int> array;		  
		  for(int f=0; f<num_bins; f++)
		    {
		      array.push_back(0);
		    }

		  int max(0),min(0);
		  string most_unlikely="";
		  int uscore=0;
		  do_counts(proteomes.at(i),proteomes.at(j),tag_length,array,max,min,most_unlikely,uscore,SOS_CUTOFF);
		  do_counts(proteomes.at(j),proteomes.at(i),tag_length,array,max,min,most_unlikely,uscore,SOS_CUTOFF);
		  //avg result;
		  //result = calc_running_avgs(array,max,min);
		  
		  v vtemp;
		  for(int f=0; f<array.size(); f++)
		    {
		      vtemp.counts.push_back(array[f]);
		      vtemp.id1=i;
		      vtemp.id2=j;
		      vtemp.uscore=uscore;
		      vtemp.most_unlikely=most_unlikely;
		    }
#pragma omp critical
		  stored_counts.push_back(vtemp);
		}		  
	    }

	  cout<<"Size of stored_counts = "<<stored_counts.size()<<endl;
	  
	  for(int r=0; r<proteomes.size()-1; r++)
	    {
	      for(int q=r+1; q<proteomes.size(); q++)
		{
		  for(int i=0; i<stored_counts.size(); i++)
		    {
		      if(stored_counts.at(i).id1==r && stored_counts.at(i).id2==q)
			{
			  cout<<stored_counts.at(i).id1<<" "<<stored_counts.at(i).id2<<endl;
			  outstream<<stored_counts.at(i).id1<<" "<<stored_counts.at(i).id2<<endl;
			  for(int j=0; j<stored_counts.at(i).counts.size(); j++)
			    {
			      outstream<<stored_counts.at(i).counts.at(j)<<endl;
			    }
			  outstream<<endl;
			}
		    }		  
		}
	    }
	}
      else if(PS=="S")
	{
	  for(int i=0; i<proteomes.size()-1; i++)
	    {
	      for(int j=i+1; j<proteomes.size(); j++)
		{
		  cout<<i<<" to "<<j<<endl;
		  int array[100];
		  for(int f=0; f<100; f++)
		    {
		      array[f]=0;
		    }
		  
		  int max(0),min(0);
		  string most_unlikely="";
		  int uscore=0;
		  //do_counts(proteomes.at(i),proteomes.at(j), tag_length, array,max,min,most_unlikely,uscore);
		  //do_counts(proteomes.at(j),proteomes.at(i), tag_length, array,max,min,most_unlikely,uscore);
		  //avg result;
		  //result = calc_running_avgs(array,max,min);
		  
		  outstream<<i<<" to "<<j<<endl;
		  //outstream<<result.index<<endl<<result.avg<<endl;
		  for(int f=0; f<100; f++)
		    {
		      outstream<<array[f]<<endl;
		    }

		  outstream<<endl;
		  outstream3<<i<<" to "<<j<<endl;
		  outstream3<<uscore<<endl<<most_unlikely<<endl<<endl;
		}		  
	    }	  
	}
    }      

  //b. read them in from file
  else if(FM=="F")
    {
      vector<Fproteome> proteomes;
    }
  else
    {
      cout<<"Invalid option for File/Memory (FM)."<<endl;
      exit(1);
    }

  //3. do all comparisons and write out.  
}

//void do_counts(Mproteome p1, Mproteome p2, int tag_length, int array[],int& max, int& min,string& most_unlikely,int& uscore)
void do_counts(Mproteome p1, Mproteome p2, int tag_length, vector<int>& array,int& max, int& min,string& most_unlikely,int& uscore,int SOS_CUTOFF)
{
  int sum=0;
  for(int i=0; i<array.size(); i++)
    {
      sum+=array[i];
    }
  bool second_do=false;
  if(sum>0)
    {
      second_do=true;
    }

  amino_acids aas;
  for(int i=0; i<p1.proteins.size(); i++)
    {
      for(int j=0; j<p1.proteins.at(i).length(); j++)
	{
	  string s="";
	  double score=0.0;
	  for(int k=0; k<tag_length && (j+k)<p1.proteins.at(i).length(); k++)
	    {
	      //	      bool composition_ok=check_composition(p1.proteins.at(i).substr(k,tag_length),tag_length,aas,SOS_CUTOFF);
	      //	      if(composition_ok)
	      //		{
		  s+=p1.proteins.at(i)[j+k];
		  int index=aas.get_aa(p1.proteins.at(i)[j+k]);
		  
		  if(index>=0 && index<20)
		    {
		      score+=p1.bit_scores[index]+p2.bit_scores[index];
		      array[(int)(score+.5)]+=1;
		      if( (score+.5)<min || min==0)
			{
			  min=score+.5;
			  //cout<<"min score = "<<min<<" sequence = "<<s<<endl;
			}
		      if( (score+.5)>max)
			{
			  max=score+.5;
			  //cout<<"max score = "<<max<<" sequence = "<<s<<endl;
			  uscore=max;
			  most_unlikely=s;
			}
		    }
		  else
		    {
		      k=tag_length;
		      j+=k-1;
		    }
		  //		}
		  //	      else
		  //		{
		  //		  k=tag_length;
		  //		  j+=k-1;
		  //		}
	    }
	}
    }
}

bool check_composition(string s,int tag_length,amino_acids& aas,int SOSCUTOFF)
{
  if(s.length()<tag_length)
    {
      return false;
    }

  vector<int> array;
  for(int i=0; i<20; i++) 
    {
      array.push_back(0);
    }

  for(int i=0; i<s.length(); i++)
    {
      int index=aas.get_aa(s[i]);
      if(index>=0 && index<20)
	{
	  array[index]++;
	}
      else
	{
	  return false;
	}
    }

  int result=0;
  for(int i=0; i<array.size(); i++)
    {
      result+=array[i]*array[i];
    }

  if(result>SOSCUTOFF)
    {
      return false;
    }

  return true;
}

avg calc_running_avgs(int array[],int max, int min)
{
  /*
  avg chosen_avg;
  chosen_avg.avg=0;
  chosen_avg.index=0;

  cout<<"In calc_running_avgs"<<endl;
  for(int i=0; i<1000; i++)
    {
      cout<<i<<" "<<array[i]<<endl;
    }

  vector<avg> rolling_avg;

  //go forwards
  double denom=0;
  double sum=0;
  double cur_avg=0;
  for(int i=min; i<max; i++)
    {
      if(array[i]>0)
	{
	  sum+=array[i];
	  denom+=1;
	  cur_avg=sum/denom;
	  avg atemp;
	  atemp.avg=log(cur_avg);
	  atemp.index=i;
	  rolling_avg.push_back(atemp);
	}
    }

  for(int i=rolling_avg.size()-1; i>=0; i--)
    {
      double diff=log(array[rolling_avg.at(i).index])-rolling_avg.at(i).avg;
      //double diff=log(counts.at(rolling_avg.at(i).index))-rolling_avg.at(i).avg;
      //cout<<i<<" "<<rolling_avg.at(i).index<<" "<<rolling_avg.at(i).avg<<" "<<diff<<endl;
      if(abs(diff)<0.08)
	{
	  chosen_avg=rolling_avg.at(i);
	  //cout<<"Choosing the avg. This better only appear once."<<endl;
	  break;
	}
    }
  //cout<<"chosen avg = "<<chosen_avg.index<<" "<<chosen_avg.avg<<endl;
  return chosen_avg;
  */
}

int get_num_bins(string path,string bfn)
{
  ifstream instream;
  instream.open( (path+"/"+bfn+"/DATA/data.txt").c_str());
  if(instream.fail())
    {
      cout<<"Failed to open data to get num bins."<<endl;
      exit(1);
    }

  int result;
  instream>>result;
  return result;
}

int calc_sos_cutoff(string path,string bfn)
{
  int result=0;
  
  ifstream instream;
  string infile = path+"/"+bfn+"/"+bfn+"_info.txt";
  instream.open(infile.c_str());
  if(instream.fail())
    {
      cout<<"Failed to open path to "<<infile<<" in calc_sos_cutoff(string,string)."<<endl;
      exit(1);
    }

  instream>>result>>result>>result;
  result=result*6.5;
  instream.close();
  return result;
}