/*Genome problem
run on super computer
Tokyo Institute of technology
http://www.gsic.titech.ac.jp/supercon/supercon2003e/shiba.c
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>

#define MAXLENG   3000
#define MAXDATA   1600
#define ASIZE 23
#define INFTY -70000000
#define SHIFT  10000
#define GOPEN -100000
#define GEXT  -10000
#define GSCO  25
#define PFILE "proteins.txt"
//#define PFILE "proteins-smallexample.txt"

int TBL[ASIZE+1][ASIZE+1] = {
  { -90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000,-90000 },
  { -90000, 40001,-10000,-20000,-20000,     0,-10000,-10000,     0,-20000,-10000,-10000,-10000,-10000,-20000,-10000, 10000,     0,-30000,-20000,     0,-20000,-10000,     0 },
  { -90000,-10000, 50001,     0,-20000,-30000, 10000,     0,-20000,     0,-30000,-20000, 20000,-10000,-30000,-20000,-10000,-10000,-30000,-20000,-30000,-10000,     0,-10000 },
  { -90000,-20000,     0, 60001, 10000,-30000,     0,     0,     0, 10000,-30000,-30000,     0,-20000,-30000,-20000, 10000,     0,-40000,-20000,-30000, 30000,     0,-10000 },
  { -90000,-20000,-20000, 10000, 60001,-30000,     0, 20000,-10000,-10000,-30000,-40000,-10000,-30000,-30000,-10000,     0,-10000,-40000,-30000,-30000, 40000, 10000,-10000 },
  { -90000,     0,-30000,-30000,-30000, 90001,-30000,-40000,-30000,-30000,-10000,-10000,-30000,-10000,-20000,-30000,-10000,-10000,-20000,-20000,-10000,-30000,-30000,-20000 },
  { -90000,-10000, 10000,     0,     0,-30000, 50001, 20000,-20000,     0,-30000,-20000, 10000,     0,-30000,-10000,     0,-10000,-20000,-10000,-20000,     0, 30000,-10000 },
  { -90000,-10000,     0,     0, 20000,-40000, 20000, 50001,-20000,     0,-30000,-30000, 10000,-20000,-30000,-10000,     0,-10000,-30000,-20000,-20000, 10000, 40000,-10000 },
  { -90000,     0,-20000,     0,-10000,-30000,-20000,-20000, 60001,-20000,-40000,-40000,-20000,-30000,-30000,-20000,     0,-20000,-20000,-30000,-30000,-10000,-20000,-10000 },
  { -90000,-20000,     0, 10000,-10000,-30000,     0,     0,-20000, 80001,-30000,-30000,-10000,-20000,-10000,-20000,-10000,-20000,-20000, 20000,-30000,     0,     0,-10000 },
  { -90000,-10000,-30000,-30000,-30000,-10000,-30000,-30000,-40000,-30000, 40001, 20000,-30000, 10000,     0,-30000,-20000,-10000,-30000,-10000, 30000,-30000,-30000,-10000 },
  { -90000,-10000,-20000,-30000,-40000,-10000,-20000,-30000,-40000,-30000, 20000, 40001,-20000, 20000,     0,-30000,-20000,-10000,-20000,-10000, 10000,-40000,-30000,-10000 },
  { -90000,-10000, 20000,     0,-10000,-30000, 10000, 10000,-20000,-10000,-30000,-20000, 50001,-10000,-30000,-10000,     0,-10000,-30000,-20000,-20000,     0, 10000,-10000 },
  { -90000,-10000,-10000,-20000,-30000,-10000,     0,-20000,-30000,-20000, 10000, 20000,-10000, 50001,     0,-20000,-10000,-10000,-10000,-10000, 10000,-30000,-10000,-10000 },
  { -90000,-20000,-30000,-30000,-30000,-20000,-30000,-30000,-30000,-10000,     0,     0,-30000,     0, 60001,-40000,-20000,-20000, 10000, 30000,-10000,-30000,-30000,-10000 },
  { -90000,-10000,-20000,-20000,-10000,-30000,-10000,-10000,-20000,-20000,-30000,-30000,-10000,-20000,-40000, 70001,-10000,-10000,-40000,-30000,-20000,-20000,-10000,-20000 },
  { -90000, 10000,-10000, 10000,     0,-10000,     0,     0,     0,-10000,-20000,-20000,     0,-10000,-20000,-10000, 40001, 10000,-30000,-20000,-20000,     0,     0,     0 },
  { -90000,     0,-10000,     0,-10000,-10000,-10000,-10000,-20000,-20000,-10000,-10000,-10000,-10000,-20000,-10000, 10000, 50001,-20000,-20000,     0,-10000,-10000,     0 },
  { -90000,-30000,-30000,-40000,-40000,-20000,-20000,-30000,-20000,-20000,-30000,-20000,-30000,-10000, 10000,-40000,-30000,-20000,110001, 20000,-30000,-40000,-30000,-20000 },
  { -90000,-20000,-20000,-20000,-30000,-20000,-10000,-20000,-30000, 20000,-10000,-10000,-20000,-10000, 30000,-30000,-20000,-20000, 20000, 70001,-10000,-30000,-20000,-10000 },
  { -90000,     0,-30000,-30000,-30000,-10000,-20000,-20000,-30000,-30000, 30000, 10000,-20000, 10000,-10000,-20000,-20000,     0,-30000,-10000, 40001,-30000,-20000,-10000 },
  { -90000,-20000,-10000, 30000, 40000,-30000,     0, 10000,-10000,     0,-30000,-40000,     0,-30000,-30000,-20000,     0,-10000,-40000,-30000,-30000, 40001, 10000,-10000 },
  { -90000,-10000,     0,     0, 10000,-30000, 30000, 40000,-20000,     0,-30000,-30000, 10000,-10000,-30000,-10000,     0,-10000,-30000,-20000,-20000, 10000, 40001,-10000 },
  { -90000,     0,-10000,-10000,-10000,-20000,-10000,-10000,-10000,-10000,-10000,-10000,-10000,-10000,-10000,-20000,     0,     0,-20000,-10000,-10000,-10000,-10000,-10000 }
};

int ALL[MAXDATA+1][MAXLENG+1];
int LENG[MAXDATA+1];
int SCO[MAXDATA+1][MAXDATA+1];
int MAXFAMILY=0;
int myrank;

int loadfile(int);
int score(int, int, int, int);
int find_family(int *,int,int,int *);
int ffopen();
int ffput(int*,int);
int ffclose();

int max(x, y)     { return x > y ? x : y; }
int max3(x, y, z) { return max(max(x, y), z); }
int min(x, y)     { return x < y ? x : y; }

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

  int num, i, j, k1, k2, count;
  int fam[MAXDATA+1];
  int famid[MAXDATA+1];
  int buf[MAXDATA+1];
  int t,p;
  MPI_Status status;

  MPI_Init( &argc, &argv );
  MPI_Comm_rank( MPI_COMM_WORLD, &myrank );
  MPI_Comm_size( MPI_COMM_WORLD, &p );

  if(myrank==0)ffopen();

  num=loadfile(myrank);

  if(myrank==0){
    k1=p-1;
    while(1){
      MPI_Recv(&buf,MAXDATA+1,MPI_INT,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status);
      if(status.MPI_TAG==1){
        MPI_Send(&k1,1, MPI_INT, status.MPI_SOURCE, 1, MPI_COMM_WORLD);
	k1++;
        if(k1>=num+p)break;
      }else{
        MPI_Get_count(&status,MPI_INT,&count);
        printf("root:receive:from %d:",status.MPI_SOURCE);
          for(j=0;j<count;j++)
            printf(" %d",buf[j]);
        printf("\n");
        if(MAXFAMILY<count){
          ffput(buf,count);
          MAXFAMILY=count;
        }
        MPI_Send(&MAXFAMILY,1, MPI_INT, status.MPI_SOURCE, 2, MPI_COMM_WORLD);
      }
    }
  }else{
    memset(SCO,0,sizeof(SCO));
    for(k1=myrank-1;k1<num;){
      printf("%d:NUM:%d\n",myrank,k1);
      i=0;
      for(k2=k1+1;k2<num;k2++){
        if(score(k1,k2,LENG[k1],LENG[k2])){
          fam[i]=k2;
          i++;
        }
      }
      if(i){
        famid[1]=k1;
        t=find_family(fam,i,1,famid);
        MAXFAMILY=max(t,MAXFAMILY);
      }
      MPI_Send(&k1,1, MPI_INT, 0, 1, MPI_COMM_WORLD);
      MPI_Recv(&k1,1,MPI_INT,0,1,MPI_COMM_WORLD,&status);
    }
  }

  MPI_Barrier( MPI_COMM_WORLD );

  if(myrank==0)ffclose();

  MPI_Finalize();

  return 0;
}

int find_family(int *fam,int k,int famnum,int *famid){
  int t,i,k1,k2,maxs=0;
/*  int j;*/
  int fam2[k];

  famnum++;

  *(famid+famnum)=*fam;

  if(k==1)return famnum;
  if(MAXFAMILY>=famnum+k)return famnum;

  for(k1=0;k1<k;k1++){
    i=0;
    for(k2=k1+1;k2<k;k2++){
      if(score(*(fam+k1),*(fam+k2),LENG[*(fam+k1)],LENG[*(fam+k2)])){
        fam2[i]=*(fam+k2);
        i++;
      }
    }
    if(i){
      *(famid+famnum)=*(fam+k1);
      t=find_family(fam2,i,famnum,famid);
      if(t>maxs){
        maxs=t;
        if(MAXFAMILY<t){
          MPI_Status status;
/*          printf("send %d:",myrank);
          for(j=1;j<=t;j++)
            printf(" %d",*(famid+j));
          printf("\n");*/
          MPI_Send(&famid[1],t, MPI_INT, 0, 2, MPI_COMM_WORLD);
          MPI_Recv(&MAXFAMILY,1,MPI_INT,0,2,MPI_COMM_WORLD,&status);
        }
      }
    }
  }

  return max(maxs,famnum);
}

int score(int k1,int k2,int m,int n) {

  int t,tmp,maxs,ans;
  int *pLSTB,*pLSTBs,*pLSTBe;
  int LSTB[m+1][3];
  int *pk1,*pk1s,*pk1e;
  int *pk2,*pk2s,*pk2e;

  if(SCO[k1][k2]==2)return 0;
  if(SCO[k1][k2]==3)return 1;

  pLSTBs=&LSTB[1][0];
  pLSTBe=&LSTB[m][0];
  pk1s=&ALL[k1][1];
  pk1e=&ALL[k1][m];
  pk2s=&ALL[k2][1];
  pk2e=&ALL[k2][n];

  for(pk1=pk1e,pLSTB=pLSTBe; pk1>=pk1s; pk1--,pLSTB-=3){
    *pLSTB=*(pLSTB+2) = TBL[*pk1][*pk2e];
  }
  maxs=*pLSTBs;

  for(pk2=pk2e-1;pk2>=pk2s;pk2--) {
    tmp=*(pLSTBe+1)=TBL[*pk1e][*pk2];
    for(pk1=pk1e-1,pLSTB=pLSTBe-3;pk1>=pk1s;pk1--,pLSTB-=3) {
      t=*(pLSTB+3)+TBL[*pk1][*pk2];
      *(pLSTB+3) = tmp;
      tmp = max3(t,*(pLSTB+4)+GOPEN,*(pLSTB+2)+GOPEN);
      *(pLSTB+1) = max(t,*(pLSTB+4)+GEXT);
      *(pLSTB+2) = max(t,*(pLSTB+2)+GEXT);
    }
    if(tmp>maxs)maxs=tmp;
  }
  for(pLSTB=pLSTBe;pLSTB>pLSTBs;pLSTB-=3){
    if(*pLSTB>maxs)maxs=*pLSTB;
  }
  if(maxs < 0) ans= maxs + (int)(-(float)maxs/(float)SHIFT + 0.5) * SHIFT;
  else         ans= maxs % SHIFT;
  ans=ans*100/min(m, n);

  SCO[k1][k2]=(ans>=GSCO)+2;

  return (ans>=GSCO);
}

int loadfile(int myrank){

  int C2N[26] = {
     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 leng, tleng, mleng;
  int num,i;

  FILE *fp;
  char instr[MAXLENG];

  fp = fopen(PFILE, "r");
  if(!fp) {
    printf("File open error!\n");
    return 0;
  }
  num = 0;  tleng = 0;
  while(1) {
    fgets(instr, MAXLENG + 1, fp);
    if(instr[0] == '@') break;
    leng = strlen(instr);
    for(i = 0; i < leng; i++) {
      if(instr[i] == '%') break;
      ALL[num][i+1] = C2N[instr[i] - 'A'];
    }
    LENG[num] = i;
    tleng = tleng + i;
    num++;
  }
  fclose(fp);
  mleng = tleng / num;

  if(myrank==0)printf("number of seq.s = %d,  average length = %d\n", num, mleng);

  return num;
}