/*  'p64MP.c' program in C performs local alignments of protein sequence */
/*   uses MacMP interface by Viktor K. Decyk to run on 2 processors */
/*   2001 J. M. Ward*/

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

#include "MacMP.h"

#define SIZE 4752    /* number of sequences */
#define LENGTH 5200     /* max length of sequence */

int nproc;

void dynamatch (const char **st1, const char **st2, int *score);

int sl1, sl2, sl3;
int idtask, nargs;
double fact1, fact2;
clock_t start_program, end_program;


int main (void)
{
int i, j, num1, num2, score1, score2, nproc; 
char line [LENGTH + 1];
char *seqs[SIZE];

typedef void (*aProcPtr) ();
FILE  *ifp, *ofp;
start_program = clock ();
MP_Init (&nproc);
MP_Setstack(200000);
printf("Processors: %d \n", nproc);

printf ("Pairwise Alignment Score Maxtrix Program \n");

if ((ifp = fopen("hydro.dat", "r")) == NULL) {printf ("File could not be opened\n"); }
ofp = fopen("matrixar.a", "w"); 

i = 0;
while (!feof(ifp)) {                      

  fgets (line, sizeof(line), ifp);
  if (!(line[0] == '>')) {
    line[strlen(line) - 1] = '\0';          /* gets rid of the new line */
    seqs[i] = malloc((strlen(line) + 1));
    strcpy(seqs[i],line); 
	i++;
  }
}
printf ("Number or records = %d \n", (i -1));
fclose(ifp); 

nargs = 3;

  for (j = 0; j < SIZE; j++) {     
     sl1 = strlen(seqs[j]);
	
     for (i = 0; i <= j; i = i+2) {	 
	   sl2 = strlen(seqs[i]);
	   if (sl1 > sl2) 
         fact1 = 1 + (log10 ((double) sl1/(double) sl2));
       else 
          fact1 = 1 + (log10 ((double) sl2/(double) sl1));
		  
	   MP_Taskstart(&idtask, (aProcPtr)&dynamatch, &nargs, &seqs[j], &seqs[i], &score1);
	 
      
	   if (i + 1 <= j) {
	     sl3 = strlen(seqs[i+1]);
         if (sl1 > sl3) 
           fact2 = 1 + (log10 ((double) sl1/(double) sl3));
         else 
            fact2 = 1 + (log10 ((double) sl3/(double) sl1));
			
	     dynamatch(&seqs[j], &seqs[i+1], &score2);
	     num2 = score2 / fact2; 
       }
	   
       MP_Taskwait(&idtask); 
       num1 = score1 / fact1; 
	   fprintf(ofp, "%6d ", num1);  
	   if (i + 1 <= j) {fprintf(ofp, "%6d ", num2);}    
	 }
  fprintf(ofp, "\n");
 /* printf("Row: %d \n", j);*/

  }
  
MP_End();
fclose(ofp); 
end_program = clock ();
printf("\ntime in min: %7.2f\n", (end_program-start_program)/3600.);
return (0);

}

void dynamatch (const char **st1, const char **st2, int *maxscore)
{
  const char *str1, *str2;
  int i,n,j, ls1, ls2, next1, next2, tempmax; 
  int gap = 3;
  int ci ['Z'];
 
  int current[LENGTH], previous[LENGTH];
 
int phat_T75_B23[24][24] = 
{5, -6, -2, -5,   1, -3, -5,  1, -3,  0, -1,  -7, -1, -1, -3,  2,  0, -4, -3,  1, -3, -4,  0, -10, 
-6,  9, -3, -7,  -8, -2, -6, -5, -4, -6, -6,  -1, -6, -7, -7, -6, -6, -7, -6, -7, -4, -4, -6, -10, 
-2, -3, 11,  2,  -2,  2,  0, -1,  4, -3, -3,  -2, -2, -1, -4,  1, -1, -5,  2, -3,  8,  1, -1, -10, 
-5, -7,  2, 12,  -7,  0,  6, -2, -1, -5, -5,  -5, -5, -5, -5, -4, -5, -7, -4, -5,  5,  3, -4, -10, 
 1, -8, -2, -7,   7, -5, -7, -2, -7, -3, -2, -10, -2,  0, -8,  1, -1, -4, -1, -2, -4, -6, -2, -10, 
-3, -2,  2,  0,  -5,  9,  1, -2,  2, -3, -3,  -1, -1, -2, -3, -1, -3,  1,  0, -3,  1,  6, -2, -10, 
-5, -6,  0,  6,  -7,  1, 12, -3, -1, -5, -5,  -4, -5, -5, -5, -3, -5, -7, -2, -5,  2,  6, -4, -10, 
 1, -5, -1, -2,  -2, -2, -3,  9, -4, -2, -2,  -5, -1, -2, -3,  1, -1, -5, -3, -2, -1, -2, -1, -10, 
-3, -4,  4, -1,  -7,  2, -1, -4, 11, -5, -4,  -5, -4, -2, -6, -2, -4, -3,  3, -5,  2,  1, -3, -10, 
 0, -6, -3, -5,  -3, -3, -5, -2, -5,  5,  2,  -7,  3,  0, -4, -2, -1, -4, -3,  3, -4, -4,  0, -10, 
-1, -6, -3, -5,  -2, -3, -5, -2, -4,  2,  4,  -7,  2,  1, -5, -2, -1, -3, -2,  1, -3, -4,  0, -10, 
-7, -1, -2, -5, -10, -1, -4, -5, -5, -7, -7,   5, -6, -7, -4, -5, -6, -8, -4, -8, -3, -3, -6, -10, 
-1, -6, -2, -5,  -2, -1, -5, -1, -4,  3,  2,  -6,  6,  0, -5, -2,  0, -4, -2,  1, -3, -3,  0, -10, 
-1, -7, -1, -5,   0, -2, -5, -2, -2,  0,  1,  -7,  0,  6, -5, -2, -2,  0,  4, -1, -3, -3,  0, -10, 
-3, -7, -4, -5,  -8, -3, -5, -3, -6, -4, -5,  -4, -5, -5, 13, -3, -4, -6, -5, -4, -4, -4, -4, -10, 
 2, -6,  1, -4,   1, -1, -3,  1, -2, -2, -2,  -5, -2, -2, -3,  6,  1, -5, -2, -2,  0, -2, -1, -10, 
 0, -6, -1, -5,  -1, -3, -5, -1, -4, -1, -1,  -6,  0, -2, -4,  1,  3, -7, -3,  0, -2, -4, -1, -10, 
-4, -7, -5, -7,  -4,  1, -7, -5, -3, -4, -3,  -8, -4,  0, -6, -5, -7, 11,  1, -4, -5, -3, -3, -10, 
-3, -6,  2, -4,  -1,  0, -2, -3,  3, -3, -2,  -4, -2,  4, -5, -2, -3,  1, 11, -3,  0, -1, -1, -10, 
 1, -7, -3, -5,  -2, -3, -5, -2, -5,  3,  1,  -8,  1, -1, -4, -2,  0, -4, -3,  4, -4, -4,  0, -10, 
-3, -4,  8,  5,  -4,  1,  2, -1,  2, -4, -3,  -3, -3, -3, -4,  0, -2, -5,  0, -4,  7,  1, -2, -10, 
-4, -4,  1,  3,  -6,  6,  6, -2,  1, -4, -4,  -3, -3, -3, -4, -2, -4, -3, -1, -4,  1,  6, -4, -10, 
 0, -6, -1, -4,  -2, -2, -4, -1, -3,  0,  0,  -6,  0,  0, -4, -1, -1, -3, -1,  0, -2, -4, -1, -10, 
-10, -10, -10,  -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10, -10,  1};

  ci['A'] = 0;
  ci['R'] = 1;
  ci['N'] = 2;
  ci['D'] = 3;
  ci['C'] = 4;
  ci['Q'] = 5;
  ci['E'] = 6;
  ci['G'] = 7;
  ci['H'] = 8;
  ci['I'] = 9;
  ci['L'] = 10;
  ci['K'] = 11;
  ci['M'] = 12;
  ci['F'] = 13;
  ci['P'] = 14;
  ci['S'] = 15;
  ci['T'] = 16;
  ci['W'] = 17;
  ci['Y'] = 18;
  ci['V'] = 19;
  ci['B'] = 20;
  ci['Z'] = 21;
  ci['X'] = 22;
  ci['*'] = 23;
  
  str1 = *st1;
  str2 = *st2;

  ls1 = strlen(str1);
  ls2 = strlen(str2);
  
  for (j=0; j < LENGTH; ++j) {previous[j] = 0;} 
  
  tempmax = 0;
  for (i=1; i<= ls1; ++i) {
  current[0] = 0;
    for (n=1; n <= ls2; ++n) { 
	  next1 = phat_T75_B23 [ci[str1[i-1]]] [ci[str2[n-1]]];
	  next2 = next1 + previous[n-1];
	  if (current[n-1] - gap > next2) {next2 = current[n-1] - gap;}
	  if (previous[n] - gap > next2) {next2 = previous[n] - gap;}
	  current[n] = next2;
	  if (next2 > tempmax) {tempmax = next2;}	
    }
	for (j=0; j < ls2; ++j) previous[j] = current[j];
  }
  *maxscore = tempmax;
  return;
}