/* ARIADNE V.1.3

Copyright 

Richard Mott 2000

Wellcome Trust Centre For Human Genetics
Univeristy of Oxford
Roosevelt Drive
Oxford OX3 7AD
UK

The software package ARIADNE is distributed in the hope that it will be
useful, but in order that the University as a charitable foundation
protects its assets for the benefit of its educational and research
purposes, the University makes clear that no condition is made or to
be implied, nor is any warranty given or to be implied, as to the
accuracy of ARIADNE, or that it will be suitable for
any particular purpose or for use under any specific conditions, or that
the content or use of ARIADNE will not constitute
or result in infringement of third-party rights.  Furthermore, the
University disclaims all responsibility for the use which is made of 
ARIADNE.


*/

#include<stdio.h>
#include<math.h>
#include"seq_util.h"
#include"ariadne.h"
#include"gapstat.h"

#define DIAG 1
#define LEFT 2
#define ABOVE 3

typedef struct origin {
  int n, m;
} ORIGIN;

static int **S, **P, **Q;
static char **path_S, **path_P, **path_Q, **path;
static ORIGIN **OriginP, **OriginQ, **Origin;
static char *align1, *align2, *align0;
static int len1, len2, ready=0, maxlen;

void invert_s(  char *seq );
void PrintThisAlignment( char *name1, char *align1, char *name2, char *align2, char *align0, int start1, int start2 );

void
BestAffineSWSetup( SEQUENCE *seq, PROFILE *pro, int self, int gap_start, int gap_extend ) {

  int *s1, *s2, *p1, *p2, *q1, *q2, *tmp, x, *mat;
  int n, m;

  gap_start = gap_start > 0 ? gap_start : -gap_start;
  gap_extend = gap_extend > 0 ? gap_extend : -gap_extend;

  gap_start += gap_extend;


  len1 = pro->length;
  len2 = seq->len;
  maxlen = ( len1 > len2 ? len1 : len2 );

  align1 = (char*)calloc( maxlen+1, sizeof(char));
  align2 = (char*)calloc( maxlen+1, sizeof(char));
  align0 = (char*)calloc( maxlen+1, sizeof(char));

  S = (int**)calloc(len1+1,sizeof(int*))+1;
  P = (int**)calloc(len1+1,sizeof(int*))+1;
  Q = (int**)calloc(len1+1,sizeof(int*))+1;

  OriginP = (ORIGIN **)calloc(len1+1,sizeof(ORIGIN*))+1;
  OriginQ = (ORIGIN **)calloc(len1+1,sizeof(ORIGIN*))+1;
  Origin = (ORIGIN **)calloc(len1+1,sizeof(ORIGIN*))+1;

  path_S = (char**)calloc(len1+1,sizeof(char*))+1;
  path_P = (char**)calloc(len1+1,sizeof(char*))+1;
  path_Q = (char**)calloc(len1+1,sizeof(char*))+1;

  for(n=-1;n<len1;n++) {
    S[n] = (int*)calloc(len2+1,sizeof(int))+1;
    P[n] = (int*)calloc(len2+1,sizeof(int))+1;
    Q[n] = (int*)calloc(len2+1,sizeof(int))+1;

    path_S[n] = (char*)calloc(len2+1,sizeof(char))+1;
    path_P[n] = (char*)calloc(len2+1,sizeof(char))+1;
    path_Q[n] = (char*)calloc(len2+1,sizeof(char))+1;

    OriginP[n] = (ORIGIN*)calloc(len2+1,sizeof(ORIGIN))+1;
    OriginQ[n] = (ORIGIN*)calloc(len2+1,sizeof(ORIGIN))+1;
    Origin[n] = (ORIGIN*)calloc(len2+1,sizeof(ORIGIN))+1;
  }
  
  for(n=0;n<len1;n++){
    Origin[n][0].n = OriginP[n][0].n = OriginQ[n][0].n = n;
    Origin[n][-1].n = OriginP[n][-1].n = OriginQ[n][-1].n = n;
  }
  for(m=0;m<len2;m++){
    Origin[0][m].m = OriginP[0][m].m = OriginQ[0][m].m = m;
    Origin[-1][m].m = OriginP[-1][m].m = OriginQ[-1][m].m = m;
  }

  for(n=0;n<len1;n++){
    s1 = S[n-1]; s2 = S[n];
    p1 = P[n-1]; p2 = P[n];
    q1 = Q[n-1]; q2 = Q[n];

    mat = pro->profile[n];

    for(m=0;m<len2;m++) {
      
      p2[m] = 0.0;
      x = p1[m] - gap_extend;
      if ( p2[m] < x ) {
        p2[m] = x;
        path_P[n][m] = ABOVE;
	OriginP[n][m] = OriginP[n-1][m];
      }
      x = s1[m] - gap_start;
      if ( p2[m] < x ) {
        p2[m] = x;
        path_P[n][m] = DIAG;
	OriginP[n][m] = Origin[n-1][m];
      }
      if ( p2[m] == 0.0 ) {
	OriginP[n][m].n = n;
	OriginP[n][m].m = m;
      }

      q2[m] = 0.0;
      x = q2[m-1] - gap_extend;
      if ( q2[m] < x ) {
        q2[m] = x;
        path_Q[n][m] = LEFT;
	OriginQ[n][m] = OriginQ[n][m-1];
      }
      x = s2[m-1] - gap_start;
      if ( q2[m] < x ) {
        q2[m] = x;
        path_Q[n][m] = DIAG;
	OriginQ[n][m] = Origin[n][m-1];
      }      
      if ( q2[m] == 0.0 ) {
	OriginQ[n][m].n = n;
	OriginQ[n][m].m = m;
      }
      
      s2[m] = 0.0;
      x = s1[m-1] + mat[seq->s[m]];
      if ( s2[m] < x ) {
        s2[m] = x;
        path_S[n][m] = DIAG;
	Origin[n][m] = Origin[n-1][m-1];
      }
      if ( s2[m] < q2[m] ) {
        s2[m] = q2[m];
        path_S[n][m] = LEFT;
	Origin[n][m] = OriginQ[n][m];
      }
      if ( s2[m] < p2[m] ) {
        s2[m] = p2[m];
        path_S[n][m] = ABOVE;
	Origin[n][m] = OriginP[n][m];
      }
      if ( s2[m] == 0.0 ) {
	Origin[n][m].n = n;
	Origin[n][m].m = m;
      }
    }
  }
  ready = 1;
}


int BestAffineSWNext( SEQUENCE *seq, PROFILE *pro, int self, int gap_start, int gap_extend, int stats, double lambda0, double K0, double H, double alpha ) {
  
  int n, m, max_m, max_n, pat;
  int  max = -1.0, bscore;
  int pos, ms;
  char c;
  ORIGIN origin;
  double K, L, pval, theta, logkappa;

  if ( ready ) {
    
    if ( self ) {
      for(n=0;n<len1;n++){
	for(m=0;m<n;m++) {
	  if ( max < S[n][m] && Origin[n][m].n != -1 && Origin[n][m].n != 0 && Origin[n][m].m != 0) {
	    max = S[n][m];
	    max_m = m;
	    max_n = n;
	  }
	}
      }
    }
    else {
      for(n=0;n<len1;n++){
	for(m=0;m<len2;m++) {
	  if ( max < S[n][m] && Origin[n][m].n != -1 ) {
	    max = S[n][m];
	    max_m = m;
	    max_n = n;
	  }
	}
      }
    }      
    if ( max > 0.0 ) {
	
      /* backtrack */
      
      m = max_m;
      n = max_n;
      bscore = 0.0;
      pos = 0;
      while( n>=0 && m>=0 && (pat=path_S[n][m]) ) {
	if ( pat == DIAG ) {
	  ms=pro->profile[n][seq->s[m]];
	  align1[pos] = pro->seq[n];
	  align2[pos] = seq->s[m];
	  if ( toupper(align1[pos]) == toupper(align2[pos]) && ms > 0 ) 
	    align0[pos] = '|';
	  else if ( ms > 0 ) 
	    align0[pos] = ':';
	  else
	    align0[pos] = ' ';
	  pos ++;
	  bscore += ms;
	  n--;
	  m--;
	}
	else if ( pat == LEFT ) {
	  while( m>= 0 && path_Q[n][m] == LEFT ) {
	    align1[pos] = '-';
	    align2[pos] = seq->s[m];
	    align0[pos] = ' ';
	    pos++;
	    bscore -= gap_extend;
	    m--;
	  }

	  align1[pos] = '-';
	  align2[pos] = seq->s[m];
	  align0[pos] = ' ';
	  pos++;
	  bscore -= gap_start;
	  m--;
	}
	else if ( pat == ABOVE ) {
	  while( n>= 0 && path_P[n][m] == ABOVE) {
	    align2[pos] = '-';
	    align1[pos] = pro->seq[n];
	    align0[pos] = ' ';
	    pos++;

	    bscore -= gap_extend;
	    n--;
	  }
	  align2[pos] = '-';
	  align1[pos] = pro->seq[n];
	  align0[pos] = ' ';
	  pos++;

	  bscore -= gap_start;
	  n--;
	}
      }
      
      align1[pos] = align2[pos] = align0[pos] = 0;
      invert_s( align1 ); invert_s( align2 ); invert_s(align0 );

/*       printf ("\n%s\n%s\n%s\n\n", align1, align0, align2 ); */

      if ( stats ) {
	pval = EmpiricalGEM( lambda0, K0, H, alpha, len1, len2, &theta, &logkappa, max );
	K = K0*exp(logkappa);
	L = lambda0*theta;
      }
      else 
	pval = 1.0;

      fprintf(stdout, "> %s len %d from %d to %d  vs  %s len %d from %d to %d   score %d  eval %e identity %.2f%% K %e L %e H %e alpha %e\n", pro->name, pro->length, start1, stop1, protein->name, protein->len, start2, stop2, score, evalue, percent, K, L, H, alpha );

      PrintThisAlignment( pro->name, align1, seq->name, align2, align0, n+1, m+1  );
      
      /* max-out */
      
      origin = Origin[max_n][max_m];
      for(n=0;n<len1;n++){
	for(m=0;m<len2;m++) {
	  if ( Origin[n][m].n == origin.n && Origin[n][m].m == origin.m ) {
	    Origin[n][m].n = -1;
	    Origin[n][m].m = -1;
	  }
	}
      }
    }
  }
  return max;
}
  
 

void BestAffineSWCleanup() {

  int n;
  
  if ( ready ) {

    free(align1); free(align2); free(align0);

    for(n=-1;n<len1;n++) {
      free(S[n]-1);
      free(P[n]-1);
      free(Q[n]-1);

      free(path_S[n]-1);
      free(path_P[n]-1);
      free(path_Q[n]-1);

      free(OriginP[n]-1);
      free(OriginQ[n]-1);
      free(Origin[n]-1);
    }

    free(S-1);
    free(P-1);
    free(Q-1);

    free(path_S-1);
    free(path_P-1);
    free(path_Q-1);

    free(OriginP-1);
    free(OriginQ-1);
    free(Origin-1);

    ready = 0;
  }
  
}

void PrintThisAlignment( char *name1, char *align1, char *name2, char *align2, char *align0, int start1, int start2 ) {

  int i, pos1=start1, pos2 = start2;
  int len1, len2;
  int len = strlen(align1);

  start1++;
  start2++;

  len1 = start1;
  len2 = start2;
  while ( *align1 ) {
    for(i=0;i<50&&align1[i];i++)
      len1 += ( align1[i] != '-' );

    for(i=0;i<50&&align2[i];i++)
      len2 += ( align2[i] != '-' );

    printf( "%6d %-.50s %6d %s\n", start1, align1, len1, name1 ? name1 : "seq1" );
    printf( "%6s %-.50s\n", "", align0 );
    printf( "%6d %-.50s %6d %s\n\n", start2, align2, len2, name2 ? name2 : "seq2" );
    
    start1 = len1+1;
    start2 = len2+1;
    align1+=i;    
    align2+=i;    
    align0+=i;
  }
}

void invert_s(  char *seq ){
  char *s = seq;
  char *t = &seq[strlen(seq)-1];
  char u;

  while ( t >= s ) {
    u = *s;
    *s = *t;
    *t = u;
    s++;
    t--;
  }
}
