#include <stdio.h>
#include <pcre.h>
#include <string.h>
#include <math.h>
#include <search.h>

#define max(x,y) (x>y?x:y)
#define min(x,y) (x<y?x:y)

char* get_parameter(int argc, char** argv, char* param);
int   findSites(char* r, char* name, char* seq, int* positions, int* np, int* orientations, int* no, int singlestrand); 
void  chomp(char* s); 
char* complement(char* s); 
char* substr(char *string, int start, int length); 
char* nextSequence(FILE* fp, char** name, int* size, int* nextSequence_started, int* nextSequence_ended, char* nextSequence_currentLine); 


int main(int argc, char** argv) {

  int   i, len, l=0, j;
  FILE* fp1, *fp2;
  FILE *fpc, *fpo, *fpp;
  char* buff;
  
  char* seq1 = 0;
  char* seq2 = 0;

  char* name1;
  char* name2;
  
  int   size1;
  int   size2;
  
  int*  positions1;
  int*  orientations1;
  int*  positions2;
  int*  orientations2;

  int   np1;
  int   no1;
  int   np2;
  int   no2;

  char* fasta1;
  char* fasta2;
  char* re;

  char* outfile;
  char* posfile;
  char* orifile;

  int   singlestrand = 0;

  int   nextSequence_started1 = 0;
  int   nextSequence_ended1   = 0;
  char* nextSequence_currentLine1; 
 
  
  int   nextSequence_started2 = 0;
  int   nextSequence_ended2   = 0;
  char* nextSequence_currentLine2;


  if (argc < 4) {
    printf("Usage : do_fastcompare_conserved_set -fasta1 FILE -fasta2 FILE -kmer STR -cons FILE -pos FILE -ori FILE");
    exit(0);
  }

  fasta1      = get_parameter(argc, argv, "-fasta1");
  fasta2      = get_parameter(argc, argv, "-fasta2");
  re          = get_parameter(argc, argv, "-kmer");
  outfile     = get_parameter(argc, argv, "-cons");
  posfile     = get_parameter(argc, argv, "-pos");
  orifile     = get_parameter(argc, argv, "-ori");

  if (strcmp(get_parameter(argc, argv, "-singlestrand"), "") != 0)
    singlestrand      = atoi(get_parameter(argc, argv, "-singlestrand"));
  
  nextSequence_currentLine1 = (char*)malloc(50000 * sizeof(char));
  nextSequence_currentLine2 = (char*)malloc(50000 * sizeof(char));

  positions1    = (int*)malloc(5000 * sizeof(int));
  orientations1 = (int*)malloc(5000 * sizeof(int));

  positions2    = (int*)malloc(5000 * sizeof(int));
  orientations2 = (int*)malloc(5000 * sizeof(int));

  fp1 = fopen(fasta1, "r");
  if (!fp1) {
    printf("cannot open %s ..\n", fasta1);
    exit(0);
  }

  fp2 = fopen(fasta2, "r");
  if (!fp2) {
    printf("cannot open %s ..\n", fasta2);
    exit(0);
  }

  //
  //  START SCANNING SPECIES 1 AND 2
  //

  if (strcmp(outfile, "") != 0) {
    fpc = fopen(outfile, "w");
    if (!fpc) {
      printf("Could not open %s ..\n" ,outfile);
    }
  }

  if (strcmp(posfile, "") != 0) {
    fpp = fopen(posfile, "w");
    if (!fpp) {
      printf("Could not open %s ..\n" , posfile);
    }
  }

  if (strcmp(orifile, "") != 0) {
    fpo = fopen(orifile, "w");
    if (!fpo) {
      printf("Could not open %s ..\n" , orifile);
    }
  }
  
  while (seq1 = nextSequence(fp1, &name1, &size1, &nextSequence_started1, &nextSequence_ended1, nextSequence_currentLine1))  {   
    
    seq2 = nextSequence(fp2, &name2, &size2, &nextSequence_started2, &nextSequence_ended2, nextSequence_currentLine2);
    
    findSites(re, name1, seq1, positions1, &np1, orientations1, &no1, singlestrand);

    findSites(re, name2, seq2, positions2, &np2, orientations2, &no2, singlestrand);

    if ((np1 > 0) && (np2 > 0)) {
      
      if (strcmp(outfile, "") != 0) {
	fprintf(fpc, "%s\n", name1);
      }

      if (strcmp(posfile, "") != 0) {

	for (i=0; i<np1; i++) {
	  fprintf(fpp, "%s\n", positions1[i]);
	}

      }

      if (strcmp(orifile, "") != 0) {

	for (i=0; i<no1; i++) {
	  fprintf(fpo, "%s\n", orientations1[i]);
	}

      }

    }
    
    free(name2);
    free(seq2);
    free(name1);
    free(seq1);
  
  }

  if (strcmp(outfile, "") != 0)
    fclose(fpc);

  if (strcmp(posfile, "") != 0)
    fclose(fpp);
  
  if (strcmp(orifile, "") != 0)
    fclose(fpo);
      
  return 0;


}
  
int findSites(char* r, char* name, char* seq, int* positions, int* np, int* orientations, int* no, int singlestrand) 
{
  
  int i, j; //, j, n;
  int   score1, score2;
  pcre* re;
  const char* error;
  int   erroffset;
  int   rc;
  int   ovector[30];
  char* cr = 0;
  char*  substring_start;
  int  substring_length;
  int startoffset = 0;
  int  pos;
  char* seq_pos;


  //
  //  store the position where a RE has been found, to avoid counting twice the palindromes
  //
  seq_pos = (char*)calloc(strlen(seq), sizeof(char));

  *np = 0;
  *no = 0;
  

  
  //printf("searching for %s in %s\n", r, seq);


  //
  // ONE STRAND
  //
  
  
  re = pcre_compile(r,
		    0,
		    &error,
		    &erroffset,
		    NULL);
  
  while ( (rc = pcre_exec(re,
		 NULL,
		 seq,
		 strlen(seq),
		 startoffset,
		 0,
		 ovector,
			  30) ) > 0) {
  

    substring_start     = seq         + ovector[0]; 
    substring_length    = ovector[1]  - ovector[0]; 
    pos                 = strlen(seq) - ovector[0] - 1;
    // save the position
    positions[ *np ]    = pos;
    // save the orientation
    orientations[ *no ] = 1;
    startoffset         = ovector[0]  + substring_length;
    seq_pos[ pos ]      = 1;

    

    (*np)++;
    (*no)++;

  }


  //return 0;

  if (singlestrand == 0) {

    // OTHER STRAND
    
    cr = complement(r);
    
    //printf("Complement of %s is %s\n", r, cr);
    
    startoffset = 0;
    
    re = pcre_compile(cr,
		      0,
		      &error,
		      &erroffset,
		      NULL);
    
    while ( (rc = pcre_exec(re,
			    NULL,
			    seq,
			    strlen(seq),
			    startoffset,
			    0,
			    ovector,
			    30)) > 0) {
      substring_start     = seq + ovector[0]; 
      substring_length    = ovector[1] - ovector[0]; 
      pos                 = strlen(seq) - ovector[0] - 1;

      //  if not a palindromic version
      if (seq_pos[ pos ] != 1) {
	
	// save the position
	positions[ *np ] = pos;
	
	
	
	(*np)++;
      }  
      
      // save the orientation in all cases
      orientations[ *no ] = -1;
      (*no)++;
      
      startoffset = ovector[0] + substring_length;
      
      
    }

  }

    
  if (cr != 0) 
    free(cr);
  free(seq_pos);
  pcre_free(re);
  return 0;
}


//
// take a char* and replace every nt by its complement
//
char* complement(char* s) 
{

  int l = strlen(s);
  char* c;
  char  d;
  int i;
  
  c = (char*)calloc(l+1, sizeof(char));

  for (i=l-1; i>=0; i--) {
    if (s[i] == 'A')
      d = 'T';
    else if (s[i] == 'T')
      d = 'A';
    else if (s[i] == 'G')
      d = 'C';
    else if (s[i] == 'C') 
      d = 'G';
    else if (s[i] == '[')
      d = ']';
    else if (s[i] == ']')
      d = '[';
    else
      d = s[i];
    
    strncat(c, &d, 1);
  }
	 
  return c;
  
}

void chomp(char* s) 
{

  int j;

  j = strlen(s);
  while (j >= 0) {
    if (s[j] == '\n') {
       s[j] = '\0';
       return;
    }
    j--;
  }
}



//
//  returns the next sequence in the file fp
//
char* nextSequence(FILE* fp, char** name, int* size, int* nextSequence_started, int* nextSequence_ended, char* nextSequence_currentLine) 
{
  
  int len = 20000;
  char* seq;
  

  // not yet started, get the first line
  if (*nextSequence_started == 0) {
    fgets(nextSequence_currentLine, len, fp);  
    chomp(nextSequence_currentLine);      
    *nextSequence_started = 1;
 
    
 }
  
  //printf("%s\n", nextSequence_currentLine);

  // if file is finished, return 0
  if (*nextSequence_ended == 1) {
    return 0;
  }

  // if started, line should be filled with ">..."
  if (nextSequence_currentLine[0] == '>') {
    

    *name = strdup(nextSequence_currentLine + 1);
    
    // create a new line
    seq = (char*)calloc(100000, sizeof(char) ); 

    while (1) {
      
      fgets(nextSequence_currentLine, len, fp);    
      
      if (feof(fp)) {
	*nextSequence_ended  = 1;
	return seq;
      }

      chomp(nextSequence_currentLine);      
      if (strlen(nextSequence_currentLine) == 0) { 
	continue;
      }
      
      if (nextSequence_currentLine[0] == '>') {	
	return seq;	
      } else {
	strcat(seq, nextSequence_currentLine);
      }
      
    }

  } else return 0;
  
}


char* get_parameter(int argc, char** argv, char* param)
{
  int i = 0;
  while ((i < argc) && (strcmp(param, argv[i])))
    i++;
  if (i<argc)
    return (argv[i+1]);
  else
    return "";
}


char *substr(char *string, int start, int length) 
{
  char *substring;
  
  substring = (char*)calloc(length+1, sizeof(char));

  memcpy(substring, string+start, length*sizeof(char));
  substring[length] = '\0';

  return substring;
}