/* program:  二つのアミノ酸配列のアラインメントと列一致度の計算
   file:     align.c
   date:     2003/7/23 by O.W
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define MAXLENG 3000    /* 配列長最大値 */
#define ASIZE 23        /* アミノ酸の種類（種類は 22 個，他の分類が 1 個）*/
#define INFTY -70000000 /* 無限値（ありえないスコア） NEW */
#define SHIFT  10000    /* スコアの値をシフトするための係数 NEW */
#define GOPEN -100000   /* ギャップを開くときのペナルティ * SHIFT NEW */
#define GEXT  -10000    /* ギャップを続けるときのペナルティ * SHIFT NEW */

/* スコア表 BLOSUM-62 */
int TBL[ASIZE+1][ASIZE+1] = {
  { -9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9 },
  { -9, 4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-2,-1, 0 },
  { -9,-1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-1, 0,-1 },
  { -9,-2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3, 3, 0,-1 },
  { -9,-2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1, 0,-1,-4,-3,-3, 4, 1,-1 },
  { -9, 0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-3,-3,-2 },
  { -9,-1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2, 1, 0,-3,-1, 0,-1,-2,-1,-2, 0, 3,-1 },
  { -9,-1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3, 1,-2,-3,-1, 0,-1,-3,-2,-2, 1, 4,-1 },
  { -9, 0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2, 0,-2,-2,-3,-3,-1,-2,-1 },
  { -9,-2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2, 2,-3, 0, 0,-1 },
  { -9,-1,-3,-3,-3,-1,-3,-3,-4,-3, 4, 2,-3, 1, 0,-3,-2,-1,-3,-1, 3,-3,-3,-1 },
  { -9,-1,-2,-3,-4,-1,-2,-3,-4,-3, 2, 4,-2, 2, 0,-3,-2,-1,-2,-1, 1,-4,-3,-1 },
  { -9,-1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2, 5,-1,-3,-1, 0,-1,-3,-2,-2, 0, 1,-1 },
  { -9,-1,-1,-2,-3,-1, 0,-2,-3,-2, 1, 2,-1, 5, 0,-2,-1,-1,-1,-1, 1,-3,-1,-1 },
  { -9,-2,-3,-3,-3,-2,-3,-3,-3,-1, 0, 0,-3, 0, 6,-4,-2,-2, 1, 3,-1,-3,-3,-1 },
  { -9,-1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4, 7,-1,-1,-4,-3,-2,-2,-1,-2 },
  { -9, 1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2, 0,-1,-2,-1, 4, 1,-3,-2,-2, 0, 0, 0 },
  { -9, 0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1, 1, 5,-2,-2, 0,-1,-1, 0 },
  { -9,-3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1, 1,-4,-3,-2,11, 2,-3,-4,-3,-2 },
  { -9,-2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1, 3,-3,-2,-2, 2, 7,-1,-3,-2,-1 },
  { -9, 0,-3,-3,-3,-1,-2,-2,-3,-3, 3, 1,-2, 1,-1,-2,-2, 0,-3,-1, 4,-3,-2,-1 },
  { -9,-2,-1, 3, 4,-3, 0, 1,-1, 0,-3,-4, 0,-3,-3,-2, 0,-1,-4,-3,-3, 4, 1,-1 },
  { -9,-1, 0, 0, 1,-3, 3, 4,-2, 0,-3,-3, 1,-1,-3,-1, 0,-1,-3,-2,-2, 1, 4,-1 },
  { -9, 0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2, 0, 0,-2,-1,-1,-1,-1,-1 }
};
/* 補足：
  1. BLAST http://www.ncbi.nlm.nih.gov/blast/html/sub_matrix.html より
  2. アルファベット（アミノ酸記号）とインデックスの対応は以下の表の通り
*/
/* アルファベット表（X は「その他」用）*/
int C2N[26] = { /* char --> index */
   1,21, 5, 4, 7,14, 8, 9,10, 0,12,11,13, 3, 0,15, 6, 2,16,17, 0,20,18,23,19,22
/* A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z */
};
int N2C[24] = { /* index --> char */
   '*', 'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L',
   'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'B', 'Z', 'X'
};

/* 大域変数 */
int SEQ1[MAXLENG+1], SEQ2[MAXLENG+1]; /* アミノ酸配列 */
int STB[MAXLENG+1][MAXLENG+1][3];     /* スコア計算のための表 */

void maketb(int, int);
void out(int, int);

int match(int c1, int c2) { return (c1 != 23) & (c1 == c2); }
int max(int x, int y)     { return x > y ? x : y; }
int max3(int x, int y, int z) { return max(max(x, y), z); }
int min(int x, int y)     { return x < y ? x : y; }


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

  int leng;
  int m, n, i, j;

  FILE *fp;
  char instr[MAXLENG];

  /* データ列の読み込み */
  fp = fopen("proteins-smallexample.txt", "r");
  if(!fp){
    printf("File open error!\n");
    exit(1);
  }

  fgets(instr, MAXLENG + 1, fp);
  leng = strlen(instr);
  for(i = 0; i < leng; i++) {
    if(instr[i] == '%')
      break;
    SEQ1[i+1] = C2N[instr[i] - 'A'];
  }
  m = i;
  fgets(instr, MAXLENG + 1, fp);
  leng = strlen(instr);
  for(i = 0; i < leng; i++) {
    if(instr[i] == '%')
      break;
    SEQ2[i+1] = C2N[instr[i] - 'A'];
  }
  n = i;
  fclose(fp);

  /* TBL の変換（SHIFT 倍する）NEW */
  for(i = 0; i <= ASIZE; i++)
    for(j = 0; j <= ASIZE; j++)
      TBL[i][j] = TBL[i][j] * SHIFT;

  /* alignment, seq. id.  の計算 */
  maketb(m, n);
  out(m, n);

  return 0;
}

void maketb(int m, int n) {

  int i, j, ch, ch1, ch2, t;

  /* 探索打ち切りのための境界条件 NEW */
  /* 解説：
     1. ギャップ途中で一方の列が終わることは許されない ==> 最後は必ず比較
     2. 通常の時に一方の列が終わったらペナルティ 0 で打ち切ってよい
  */
  ch = SEQ2[n];
  for(i = 1; i <= m; i++)
    STB[i][n][0] = STB[i][n][1] = TBL[SEQ1[i]][ch] + match(SEQ1[i], ch);
  ch = SEQ1[m];
  for(j = 1; j <= n; j++)
    STB[m][j][0] = STB[m][j][2] = TBL[ch][SEQ2[j]] + match(SEQ2[j], ch);

  /* STB の計算 */
  for(i = m - 1; i >= 1; i--) {
    ch1 = SEQ1[i];
    for(j = n - 1; j >= 1; j--) {
      ch2 = SEQ2[j];
      t = TBL[ch1][ch2];
      t = t + match(ch1, ch2);  // NEW
      STB[i][j][0] = max3(STB[i+1][j+1][0] + t,
			  STB[ i ][j+1][1] + GOPEN,
			  STB[i+1][ j ][2] + GOPEN);

      STB[i][j][1] = max(STB[i+1][j+1][0] + t,
			 STB[ i ][j+1][1] + GEXT);

      STB[i][j][2] = max(STB[i+1][j+1][0] + t,
			 STB[i+1][ j ][2] + GEXT);
    }
  }

  return;
}

int putout(int g, int i, int j, int k, int s, int s1, int s2, int s3,
	   char xout[], char yout[]) {
  if (s == s1) {
    xout[k] = N2C[SEQ1[i]];  yout[k] = N2C[SEQ2[j]];
    return 0;
  }

  if (s == s2) {
    xout[k] = '-';  yout[k] = N2C[SEQ2[j]];
    return 1;
  }

  if (s == s3) {
    xout[k] = N2C[SEQ1[i]];  yout[k] = '-';
    return 2;
  }

  printf("ERROR AT OUT: %d/ %d, %d; %d, %d, %d, %d\n",
	 g, i, j, s, s1, s2, s3);
  exit(-1);
}

void out(int m, int n) {
  int i, j, ii, jj, k, g;
  int s, s1, s2, t, idnum;
  int maxi, maxj, maxk;
  char xout[MAXLENG+1], yout[MAXLENG+1];

  /* 最適スタート地点を探す */
  s1 = INFTY;
  for(j = 1; j <= n; j++) {
    s = STB[1][j][0];
    if(s > s1) {
      s1 = s;  maxj = j;
    }
  }

  s2 = INFTY;
  for(i = 1; i <= m; i++) {
    s = STB[i][1][0];
    if(s > s2) {
      s2 = s;  maxi = i;
    }
  }

  if(s1 > s2) {
    s = s1;  maxi = 1;  maxk = maxj;
  } else {
    s = s2;  maxj = 1;  maxk = maxi;
  }

  /* 出力する文字列の作成（先頭部分）*/
  i = j = 1;
  for(k = 1; k < maxk; k++) {
    if(maxi == 1) {
      xout[k] = ' ';  yout[k] = N2C[SEQ2[j]];  j++;
    } else {
      xout[k] = N2C[SEQ1[i]];  yout[k] = ' ';  i++;
    }
  }

  /* 出力する文字列の作成（本体）*/
  /* 番兵の設定 NEW */
  for(i = 1; i <= m; i++) {
    STB[i][n+1][0] = 0;
    STB[i][n+1][1] = INFTY;
  }
  for(j = 1; j <= n; j++) {
    STB[m+1][j][0] = 0;
    STB[m+1][j][2] = INFTY;
  }
  STB[m+1][n+1][0] = 0;                          

  i = maxi;  j = maxj;  k = maxk - 1;
  g = 0;

  while (i <= m && j <= n) {
    t = TBL[SEQ1[i]][SEQ2[j]] + match(SEQ1[i], SEQ2[j]);
    k++;

    switch (g) {
    case 0:
      g = putout(g, i, j, k,
		 STB[ i ][ j ][g],
		 STB[i+1][j+1][0] + t,
		 STB[ i ][j+1][1] + GOPEN,
		 STB[i+1][ j ][2] + GOPEN,
		 xout, yout);
      break;

    case 1:
      g = putout(g, i, j, k,
		 STB[ i ][ j ][g],
		 STB[i+1][j+1][0] + t,
		 STB[ i ][j+1][1] + GEXT,
		 INFTY,
		 xout, yout);
      break;

    case 2:
      g = putout(g, i, j, k,
		 STB[ i ][ j ][g],
		 STB[i+1][j+1][0] + t,
		 INFTY,
		 STB[i+1][ j ][2] + GEXT,
		 xout, yout);
    }

    switch (g) {
    case 0: i++; j++; break;
    case 1: j++; break;
    case 2: i++;
    }
  }

  /* idnum を mixed score から取り出す NEW */
  if(s < 0) idnum = s + (int)(-(float)s/(float)SHIFT + 0.5) * SHIFT;
  else      idnum = s % SHIFT;

  /* 結果の出力 */
  printf("score = %d, idnum = %d, seq. id. = %d\n",
	 s, idnum, idnum * 100 / min(m, n));
  printf("SEQ1 (%4d): ", m);
  for(ii = 1; ii <= k; ii++)
    printf("%c", xout[ii]);
  for(;i <= m; i++)
    printf("%c", N2C[SEQ1[i]]);
  printf("\n");
  printf("SEQ2 (%4d): ", n);
  for(jj = 1; jj <= k; jj++)
    printf("%c", yout[jj]);
  for(;j <= n; j++)
    printf("%c", N2C[SEQ2[j]]);
  printf("\n");

  return;
}