/* yama - Yet Another Multiple Aligner */
/* Version Creation Date: 10/30/92 */
/* Revision Histories:
	March 93 : lav format
	April 93 : linked with lib.c
	May 93   : left/right end gaps
	April 94 : rewritten in Standard C
*/

#include <stdio.h>
#include <ctype.h>

#define max(x,y)  ((x) >= (y) ? (x) : (y))
#define min(x,y)  ((x) <= (y) ? (x) : (y))
#define DEFAULT_LL 50		/* output line length */
#define SIGMA 128
#define SUB(a,b)  PD(i,j,(a),(b)) = PD(i,j,(b),(a)) = \
		PD(j,i,(a),(b)) = PD(j,i,(b),(a))
#define SUBI(a,b)  PD(i,j,(a),(b)) = PD(j,i,(b),(a))
#define SUBJ(a,b)  PD(j,i,(a),(b)) = PD(i,j,(b),(a))
#define DAG(a)    PD(i,j,(a),(a)) = PD(j,i,(a),(a))
#define DT(x)	PT[i][j]x = PT[j][i]x
#define DTI(x)	PT[i][j]x
#define DTJ(x)	PT[j][i]x
#define NUMBER	10		/* Maximum number of sequences */
#define SYMA	'A'
#define SYMC	'C'
#define SYMG	'G'
#define SYMT	'T'
#define DASH	'-'
#define LSTAR	'*'
#define RSTAR	'@'
#define MININT	-9999999
#define TRUE	1
#define FALSE	0
/*
#define	USAGE   "Usage: yama [-e <endpoints file>] [-c <score file>] [-t] [L=??] <sequence order>"
*/

#define	USAGE   "Usage: yama [-c <score file>] [-t] [L=??] <sequence order>"

/*
score file format :
	seq1 seq2 60 2 8 -10
	seq2 seq3 50 2 10 -10
	default 60 2 10 -10

	seq2 3 2 4 3
	default 0 0 0 0

Default : to output LAV-format text
	-t : normal text output
	 L : text length

input sequence order format:
	"((seq1[20,50] seq2[30,45]) seq3[20,50])"
*/

/* int	PD(NUMBER+1,NUMBER+1,SIGMA,SIGMA); */	/* symbol distance */
int	*PD_storage;		/* symbol distance */
int	K;			/* number of sequences */
int	LL;			/* line length */
#define PD_SIZE(K) (sizeof(int) *  SIGMA * SIGMA * ((K)+1) * ((K)+1))
#define PD_ALLOC(K) do{PD_storage = (int*)ckalloc(PD_SIZE(K));}while(0)
#define PD(a,b,c,d) PD_storage[(d) + SIGMA*((c) + SIGMA*((b) + (K+1)*(a)))]
#define PD2(a,b) (PD_storage + SIGMA*SIGMA*((b) + (K+1)*(a)))


int	PT[NUMBER+1][NUMBER+1][2][2][4][4];	/* Altschul gap counts */
char	*S [NUMBER+1];		/* sequences */
int	L [NUMBER+1];		/* lengths of sequences */
int	From[NUMBER+1];		/* Starting position */
int	To[NUMBER+1];		/* End position */
int	SD[NUMBER+1][NUMBER+1];	/* sequence distance */
char 	*Sname[NUMBER+1];	/* sequence name */
char	Ename[30];
/*
int	gflag = 0;
*/
char	tflag = 1;
char	cflag = 0;
char	eflag = 0;

static int sl;
static int *sapp;				/* Current script append ptr */
static int  last;				/* Last script op appended */

						/* Append "Delete k" op */
#define DEL(k)				\
{ sl += k;				\
  if (last < 0)				\
    last = sapp[-1] -= (k);		\
  else					\
    last = *sapp++ = -(k);		\
}
						/* Append "Insert k" op */
#define INS(k)				\
{ sl += k;				\
  if (last > 0)				\
    last = sapp[-1] += (k);		\
  else					\
    last = *sapp++ = (k);		\
}

#define REP { ++sl; last = *sapp++ = 0; }	/* Append "Replace" op */

/* - -> * for some cases */
#define ET(ii,jj)  ((ii == DASH)? jj : ii)

/* compute the value for gap table T */
/* 0 : DASH;  1 : letter;  2 : LSTAR ; 3 : RSTAR */
#define GV(ii)  ((ii==DASH)? 0 : ((ii==LSTAR)? 2 : ((ii==RSTAR)? 3 : 1)))

#define ETGV(kk,ii,jj) {	\
	kk = ET(ii,jj);		\
	kk = GV(kk);		\
}

int 	*C1, *C2;
int 	*HH, *VV, *DD;
int 	*RH, *RV, *RD;
char	**E1, **E2;	/* End gaps indicator */
char	**HE1, **VE2;	/* End gaps indicator */


typedef struct tnode *Tptr;
typedef struct tnode {
	Tptr left;
	Tptr right;
	int L;		/* leftmost sequence */
	int R;		/* rightmost sequence */
	char **AL;	/* return alignment */
	int N;		/* alignment length */
} Tnode;
char order[256];

void main(), data(), StoAL(), ALStoAL(), MALStoMAL(), MDisplay(), 
	LDisplay(), order_print(), print_score();
Tptr build_tree();
int traverse_tree();

void main(int argc, char *argv[])
{
	int e, m, v;
	char cname[30];
	char **AL;
	int  *ALS;
	int  SL;
	char *ckalloc();
	int i,j;
	int c;
	Tptr tp;

	/*
	gflag = 0;
	*/
	LL = DEFAULT_LL;
	if (argc <= 1) fatal(USAGE);
	for ( ++argv; --argc; argv++ )
	   if ( (*argv)[0] == '-' ) {
		/*
		if ( (*argv)[1] == 'g' ) {
			if ( (*argv)[2] == '2' ) gflag = 2;
			else gflag = 1;
		} else if ( (*argv)[1] == 't' ) tflag = 0;
		*/
		if ( (*argv)[1] == 't' ) tflag = 0;
		else if ( (*argv)[1] == 'c' ) { /* score file is provided */
			++argv;
			--argc;
			sscanf(*argv,"%s",cname);
			cflag = 1;
			score_name(cname);
		} else if ( (*argv)[1] == 'e' ) { /* end points are provided */
			++argv;
			--argc;
			sscanf(*argv,"%s",Ename);
			eflag = 1;
		}
		else fatal(USAGE);
	   } else if ((*argv)[0] == 'L' && (*argv)[1] == '=')
			LL = atoi(*argv + 2);
	   else {
		strcpy(order,*argv);
		tp = build_tree(order);
	   }	
	/*
	for (i=1; i<=K; ++i) printf("Seq %d: Length: %d %s\n",i,L[i],S[i]+1);
	*/

	PD_ALLOC(K);
	if (tflag) {
		printf("#:lav\n\n");
		printf("d {\n  \"YAMA output with parameters:\n");
	}
	order_print(order);
	/*
	if (cflag) printf("  score file: %s\n",cname);
	else printf("  default score file used\n");
	*/
	data();    /* get input scoring scheme */
	c = traverse_tree(tp);
	if (tflag) {
		printf("  \"\n");
		printf("}\n");
		LDisplay(tp->AL,K-1,tp->N,c);
	} else {
		printf("\n  alignment score = %d\n\n",c);
		MDisplay(tp->AL,K-1,tp->N);
	}
#ifdef CHECK
	printf("check_score=%d\n",align_score(tp->AL,K-1,tp->N));
#endif
	free(tp->AL[0]);
	free(tp->AL);
	exit(0);
}

/* order_print - print out the sequence order */
void order_print(char *order)
{
	char *op, buf[80], flag;
	int LL = 60, IL = 54, i;
	int line_no=0;

	flag = TRUE;
	while (flag) {
		while (*order == ' ') ++order;
		i = 0;
		op = order;
		if (++line_no > 1) LL = IL;
		while ((*op++ != '\0') && (i++ < LL)); 
		op--;
		if (*op == '\0') flag = FALSE;
		if (i > 0) {
			while (*op != ']' && *op !=')' && *op != '\0' && i-- > 0) op--;
			if (i <= 0) {
				flag = FALSE;
				if (line_no == 1)
					printf("  order=%s\n",order);
				else
					printf("        %s\n",order);
				printf("Last line is too long!\n");
			} else {
				strncpy(buf,order,i);
				buf[i] = '\0';
				if (line_no == 1)
					printf("  order=%s\n",buf);
				else
					printf("        %s\n",buf);
			}
		}
		order = ++op;
	}
}

/* build_tree - return the root of the aligning tree */
Tptr build_tree(char *order)
{
	Tptr Stack[NUMBER], p;
	int i,j,k;
	int sp;
	char temp[20];
	char *ckalloc();
	char *tname[NUMBER+1];	/* sequence name */
	int tfrom[NUMBER], tto[NUMBER], tK;
	char *get_seqname();
	FILE *fp, *ckopen();
	int strsame();
	char fname[40];

	if (eflag) { /* read end points */
	   fp = ckopen(Ename,"r");
	   lav_read_init(fp);
	   tK = 0;
	   while ((tname[tK] = get_seqname(&tfrom[tK], &tto[tK],0))) 
		   if (++tK > NUMBER) fatal("The maximum number of sequences is %d",NUMBER);
	   fclose(fp);
	}
	sp = 0;
	K = 0;
	i = 0;
	while (order[i] != '\0') {
		switch (order[i]) {
		case ')':
			 if (sp < 2) fatal("wrong input format!");
			 p = (Tptr)ckalloc(sizeof(Tnode));
			 p->right = Stack[--sp];
			 p->left = Stack[--sp];
			 p->L = (p->left)->L;
			 p->R = (p->right)->R;
			 Stack[sp++] = p;
		case '(':
		case ']':
		case ' ':
		case '\t':
			 ++i;
			 break;
		default:
			if (++K > NUMBER) fatal("The maximum number of sequences is %d",NUMBER);
			j = strcspn(order+i,"( )\0\t");
			strncpy(fname,order+i,j); 
			/*
			printf("fname\n",K,fname[K]);
			*/
			i += j;
			L[K] = get_seq(fname,&S[K],&Sname[K],&From[K],&To[K]);
			S[K]--;
			/*
			   printf("K=%d, Sname=%s, From=%d, To=%d\n",K,Sname[K],From[K],To[K]);
			   fflush(stdout);
			*/
			p = (Tptr)ckalloc(sizeof(Tnode));
			p->left = p->right = NULL;
			p->L = p->R = K;
			Stack[sp++] = p;
			break;
		}
	}
	if (sp != 1) fatal("wrong format!");
	return Stack[0];
}

/* traverse_tree - traverse the aligning tree */
int traverse_tree(Tptr tp)
{
	int *ALS;
	int SL, i, l, r, n;
	int c;
	Tptr p,q;
	char *ckalloc();

	p = tp->left;
	q = tp->right;
	if (p->L == p->R) { /* one to many */
		p->AL = (char **) ckalloc((L[p->L]+1)*sizeof(char *));
		p->AL[0] = (char *) ckalloc(L[p->L]+1);
		for (i=1; i<=L[p->L]; ++i) p->AL[i] = (p->AL[i-1])+1;
		StoAL(S[p->L],L[p->L],p->AL);
		p->N = L[p->L];
	} else traverse_tree(p);
	if (q->L == q->R) { /* many to one */
		q->AL = (char **) ckalloc((L[q->L]+1)*sizeof(char *));
		q->AL[0] = (char *) ckalloc(L[q->L]+1);
		for (i=1; i<=L[q->L]; ++i) q->AL[i] = (q->AL[i-1])+1;
		StoAL(S[q->L],L[q->L],q->AL);
		q->N = L[q->L];
	} else traverse_tree(q);
	/*
	printf("(%d %d)",tp->L,tp->R);
	*/
	ALS = (int *)ckalloc((p->N + q->N +1)*sizeof(int));
	/*
	Mscore(p->AL,p->R - p->L, p->N, q->AL, q->R - q->L, q->N,p->L,q->L);
	*/
	c = MALIGN(p->AL,p->R - p->L, p->N, q->AL, q->R - q->L, q->N, ALS, &SL,p->L,q->L);
	n = tp->R - tp->L + 1;
	tp->AL = (char **)ckalloc((SL+1)*sizeof(char *));
	tp->AL[0] = (char *)ckalloc(n*(SL+1));
	for (i=1; i<=SL; ++i) tp->AL[i] = tp->AL[i-1] + n;
	MALStoMAL(p->AL,p->R - p->L, p->N, q->AL, q->R - q->L, q->N, ALS, SL, tp->AL);
	tp->N = SL;
	free(ALS);
	free(p->AL[0]);
	free(p->AL);
	free(q->AL[0]);
	free(q->AL);
	return c;
}

#define COST 		\
{	DAG(DASH) = 0;	\
	DAG(SYMA) = m;	\
	DAG(SYMC) = m;	\
	DAG(SYMG) = m;	\
	DAG(SYMT) = m;	\
	SUB(DASH,SYMA) = -GE; 	\
	SUB(DASH,SYMC) = -GE; 	\
	SUB(DASH,SYMG) = -GE; 	\
	SUB(DASH,SYMT) = -GE; 	\
	SUB(SYMA,SYMC) = v;	\
	SUB(SYMA,SYMG) = v;	\
	SUB(SYMA,SYMT) = v;	\
	SUB(SYMC,SYMG) = v;	\
	SUB(SYMC,SYMT) = v;	\
	SUB(SYMG,SYMT) = v;	\
			\
	DT([0][0][0][0]) = 0;  /* -- : -- */	\
	DT([0][0][0][1]) = G;  /* -- : -x */	\
	DT([0][0][1][0]) = G;  /* -x : -- */	\
	DT([0][0][1][1]) = 0;  /* -x : -x */	\
	DT([0][1][0][0]) = 0;  /* -- : x- */	\
	DT([0][1][0][1]) = 0;  /* -- : xx */	\
	DT([0][1][1][0]) = G;  /* -x : x- */	\
	DT([0][1][1][1]) = 0;  /* -x : xx */	\
	DT([1][0][0][0]) = 0;  /* x- : -- */	\
	DT([1][0][0][1]) = G;  /* x- : -x */	\
	DT([1][0][1][0]) = 0;  /* xx : -- */	\
	DT([1][0][1][1]) = 0;  /* xx : -x */	\
	DT([1][1][0][0]) = 0;  /* x- : x- */	\
	DT([1][1][0][1]) = G;  /* x- : xx */	\
	DT([1][1][1][0]) = G;  /* xx : x- */	\
	DT([1][1][1][1]) = 0;  /* xx : xx */	\
						\
			\
	/* Left end gaps */	\
	DAG(LSTAR) = 0;	\
	SUBI(LSTAR,SYMA) = -LGGEI; 	\
	SUBI(LSTAR,SYMC) = -LGGEI; 	\
	SUBI(LSTAR,SYMG) = -LGGEI; 	\
	SUBI(LSTAR,SYMT) = -LGGEI; 	\
	SUBJ(LSTAR,SYMA) = -LGGEJ; 	\
	SUBJ(LSTAR,SYMC) = -LGGEJ; 	\
	SUBJ(LSTAR,SYMG) = -LGGEJ; 	\
	SUBJ(LSTAR,SYMT) = -LGGEJ; 	\
	SUB(LSTAR,DASH) = 0;	\
						\
	DT([0][0][0][2]) = 0;  /* -- : -* */	\
	DTI([0][0][1][2]) = LGGJ; /* -x : -* */	\
	DTJ([0][0][1][2]) = LGGI; /* -x : -* */	\
	DT([0][0][2][0]) = 0;  /* -* : -- */	\
	DTI([0][0][2][1]) = LGGI; /* -* : -x */	\
	DTJ([0][0][2][1]) = LGGJ; /* -* : -x */	\
	DT([0][0][2][2]) = 0;  /* -* : -* */	\
	DT([0][1][0][2]) = 0;  /* -- : x* */	\
	DTI([0][1][1][2]) = LGGJ; /* -x : x* */	\
	DTJ([0][1][1][2]) = LGGI; /* -x : x* */	\
	DT([0][1][2][0]) = 0;  /* -* : x- */	\
	DT([0][1][2][1]) = 0;  /* -* : xx */	\
	DT([0][1][2][2]) = 0;  /* -* : x* */	\
	DT([1][0][0][2]) = 0;  /* x- : -* */	\
	DT([1][0][1][2]) = 0;  /* xx : -* */	\
	DT([1][0][2][0]) = 0;  /* x* : -- */	\
	DTI([1][0][2][1]) = LGGI; /* x* : -x */	\
	DTJ([1][0][2][1]) = LGGJ; /* x* : -x */	\
	DT([1][0][2][2]) = 0;  /* x* : -* */	\
	DT([1][1][0][2]) = 0;  /* x- : x* */	\
	DTI([1][1][1][2]) = LGGJ; /* xx : x* */	\
	DTJ([1][1][1][2]) = LGGI; /* xx : x* */	\
	DT([1][1][2][0]) = 0;  /* x* : x- */	\
	DTI([1][1][2][1]) = LGGI; /* x* : xx */	\
	DTJ([1][1][2][1]) = LGGJ; /* x* : xx */	\
	DT([1][1][2][2]) = 0;  /* x* : x* */	\
						\
	/* Right end gaps */			\
	DAG(RSTAR) = 0;	\
	SUBI(RSTAR,SYMA) = -RGGEI; 	\
	SUBI(RSTAR,SYMC) = -RGGEI; 	\
	SUBI(RSTAR,SYMG) = -RGGEI; 	\
	SUBI(RSTAR,SYMT) = -RGGEI; 	\
	SUBJ(RSTAR,SYMA) = -RGGEJ; 	\
	SUBJ(RSTAR,SYMC) = -RGGEJ; 	\
	SUBJ(RSTAR,SYMG) = -RGGEJ; 	\
	SUBJ(RSTAR,SYMT) = -RGGEJ; 	\
	SUB(RSTAR,DASH) = 0;	\
						\
	DT([0][0][0][3]) = 0;  /* -- : -* */	\
	DTI([0][0][1][3]) = RGGJ; /* -x : -* */	\
	DTJ([0][0][1][3]) = RGGI; /* -x : -* */	\
	DT([0][0][3][0]) = 0;  /* -* : -- */	\
	DTI([0][0][3][1]) = RGGI; /* -* : -x */	\
	DTJ([0][0][3][1]) = RGGJ; /* -* : -x */	\
	DT([0][0][3][3]) = 0;  /* -* : -* */	\
	DT([0][1][0][3]) = 0;  /* -- : x* */	\
	DTI([0][1][1][3]) = RGGJ; /* -x : x* */	\
	DTJ([0][1][1][3]) = RGGI; /* -x : x* */	\
	DT([0][1][3][0]) = 0;  /* -* : x- */	\
	DT([0][1][3][1]) = 0;  /* -* : xx */	\
	DT([0][1][3][3]) = 0;  /* -* : x* */	\
	DT([1][0][0][3]) = 0;  /* x- : -* */	\
	DT([1][0][1][3]) = 0;  /* xx : -* */	\
	DT([1][0][3][0]) = 0;  /* x* : -- */	\
	DTI([1][0][3][1]) = RGGI; /* x* : -x */	\
	DTJ([1][0][3][1]) = RGGJ; /* x* : -x */	\
	DT([1][0][3][3]) = 0;  /* x* : -* */	\
	DT([1][1][0][3]) = 0;  /* x- : x* */	\
	DTI([1][1][1][3]) = RGGJ; /* xx : x* */	\
	DTJ([1][1][1][3]) = RGGI; /* xx : x* */	\
	DT([1][1][3][0]) = 0;  /* x* : x- */	\
	DTI([1][1][3][1]) = RGGI; /* x* : xx */	\
	DTJ([1][1][3][1]) = RGGJ; /* x* : xx */	\
	DT([1][1][3][3]) = 0;  /* x* : x* */	\
				\
	SUB(LSTAR,RSTAR) = 0;		\
	DT([0][0][2][3]) = 0;  /* -* : -@ */	\
	DT([0][0][3][2]) = 0;  \
	DT([0][1][2][3]) = 0;  \
	DT([0][1][3][2]) = 0;  \
	DT([1][0][2][3]) = 0;  \
	DT([1][0][3][2]) = 0;  \
	DT([1][1][2][3]) = 0;  \
	DT([1][1][3][2]) = 0;  \
}

/* data :   read input scoring scheme */
void data ()
{
	auto		int	i, j;
	auto		int	m, v;
	auto		int	G, GE;	/* gap cost */
	auto		int	LGGI, LGGJ, LGGEI, LGGEJ;  /* end gap cost */
	auto		int	RGGI, RGGJ, RGGEI, RGGEJ;  /* end gap cost */
	auto		char	first_time;

	/*
	printf("\n                              G  GE  match  mismatch\n");
	*/
	for (i=1; i<K; ++i)
	   for (j=i+1; j<=K; ++j) {
		if (!get_DNA_scores(Sname[i],Sname[j],&m,&v,&G,&GE))
			fatal("wrong cost file!");
		/*
		printf("  %-12s %-12s %3d %3d    %3d       %3d\n",Sname[i],Sname[j],m,v,G,GE);
		*/
		if (!get_DNA_end_scores(Sname[i],&LGGI,&LGGEI,&RGGI,&RGGEI)) {
			LGGI = G;
			LGGEI = GE;
			RGGI = G;
			RGGEI = GE;
		}
		/*
		printf("i=%d, LGGI=%d, LGGEI=%d, RGGI=%d, RGGEI=%d\n",i,LGGI,LGGEI,RGGI,RGGEI);
		*/
		if (!get_DNA_end_scores(Sname[j],&LGGJ,&LGGEJ,&RGGJ,&RGGEJ)) {
			LGGJ = G;
			LGGEJ = GE;
			RGGJ = G;
			RGGEJ = GE;
		}
		/*
		printf("j=%d, LGGJ=%d, LGGEJ=%d, RGGJ=%d, RGGEJ=%d\n",j,LGGJ,LGGEJ,RGGJ,RGGEJ);
		*/
		COST
	   }
	 print_score();
	   /*
	   first_time = 1;
	   for (i=1; i<=K; ++i)
	     if (get_DNA_end_scores(Sname[i],&LGGI,&LGGEI,&RGGI,&RGGEI)) {
		if (first_time) {
	   		printf("\n  Penalties for the end gaps:\n");
	   		printf("                LG  LGE  RG  RGE\n");
			first_time = 0;
		}
		printf("  %-12s %3d  %3d %3d  %3d\n",Sname[i],LGGI,LGGEI,RGGI,RGGEI);
	     }
	   */
}

/* StoAL - transfer a sequence (S1) to a one-sequence alignment */
void StoAL(char *S1, int L1, char **AL)
{
	int i,j,k,op;

	AL[0][0] = DASH;
	for (i=1; i<=L1; ++i) {
		AL[i][0] = S1[i];
	}
}

/* ALStoAL - generate a pairwise alignment (AL) according to the script (ALS) */
void ALStoAL(char *S1, char *S2, int L1, int L2, int ALS[], int SL, char **AL)
{
	int i,j,k,op;

	AL[0][0] = DASH;
	AL[0][1] = DASH;

 	i = j = k = op = 0;
  	while (i < L1 || j < L2)
    	  { if (op == 0 && *ALS == 0)
        	{ op = *ALS++;
          	  AL[++k][0] = S1[++i];
          	  AL[k][1]  = S2[++j];
        	}
      	   else
        	{ if (op == 0)
            		op = *ALS++;
          	  if (op > 0)
            		{ AL[++k][0] = DASH;
              		  AL[k][1] = S2[++j];
              		  op--;
            		}
          	  else
            		{ AL[++k][0] = S1[++i];
              		  AL[k][1] = DASH;
              		  op++;
            		}
        	}
	}
}

/* MALStoMAL - generate a multiple alignment according to the script */
void MALStoMAL(char **AL1, int M1, int N1, char **AL2, int M2, int N2, int ALS[], int SL, char **AL)
{
	int i,j,k,op;
	char *alk, *alt;
	int I,J;
	int MM;

	MM = M1+M2+1;

	for (I=0; I<=MM; ++I) AL[0][I] = DASH;

 	i = j = k = op = 0;
  	while (i < N1 || j < N2)
    	  { if (op == 0 && *ALS == 0)
        	{ op = *ALS++;
		  ++i; ++j; ++k;
		  alk = AL[k];
		  alt = AL1[i];
		  for (I=0; I<=M1; ++I) alk[I]=alt[I];
		  alt = AL2[j];
		  for (I=0; I<=M2; ++I) alk[M1+1+I]=alt[I];
        	}
      	   else
        	{ if (op == 0)
            		op = *ALS++;
          	  if (op > 0)
			{ ++j; ++k;
		 	  alk = AL[k];
		  	  for (I=0; I<=M1; ++I) alk[I]=DASH;
		  	  alt = AL2[j];
		  	  for (I=0; I<=M2; ++I) alk[M1+1+I]=alt[I];
              		  op--;
            		}
          	  else
			{ ++i; ++k;
		 	  alk = AL[k];
		  	  alt = AL1[i];
		  	  for (I=0; I<=M1; ++I) alk[I]=alt[I];
		  	  for (I=0; I<=M2; ++I) alk[M1+1+I]=DASH;
              		  op++;
            		}
        	}
	}
}

/* MDisplay - display a multiple alignment */
void MDisplay(char **AL, int M, int N)
{
	int i,j,k;
	int POS[NUMBER];

	for (k=0; k<=M; ++k) POS[k] = From[k+1];
	for (i=1; i<=N; i=i+LL) {
	   printf("      ");
	   for (j=i+1; j<=N+1 && j<i+LL; ++j)
		if ((j-i)%10 == 0)
			putchar('|');
		else
			putchar(' ');
	   putchar('\n');
	   printf("%5d ",POS[0]);
	   for (j=i; j<=N && j<i+LL; ++j) {
		printf("%c",AL[j][0]);
		if (AL[j][0] != DASH) POS[0]++;
	   }
	   printf(" %s\n",Sname[1]);
	   for (k=1; k<=M; ++k) {
		printf("%5d ",POS[k]);
		for (j=i; j<=N && j<i+LL; ++j) {
			if (AL[j][k] != DASH) {
			   if (AL[j][k] == AL[j][0]) printf("%c",'.');
			   else printf("%c",AL[j][k]);
			   POS[k]++;
			} else printf("%c",DASH);
		}
	   	printf(" %s\n",Sname[k+1]);
	   }
	   printf("\n");
	}
}

/* LDisplay - display a multiple alignment in "lax" format */
void LDisplay(char **AL, int M, int N, int c)
{
	int i, j1, j2, k;
	int POS[NUMBER];
	char *A1, *A2;
	char same;

	printf("s {\n");
	for (k=0; k<=M; ++k)
		printf("  \"%s\" %d %d\n",Sname[k+1], From[k+1],To[k+1]);
	printf("}\n");
	printf("a {\n");
	printf("  s %d\n",c);
	printf("  b");
	for (k=0; k<=M; ++k) printf(" %d",1);
	printf("\n  e");
	for (k=0; k<=M; ++k) printf(" %d",To[k+1]-From[k+1]+1);
	printf("\n");
	j1 = 1;
	A1 = AL[1];
	j2 = j1+1;
	for (k=0; k<=M; ++k) POS[k] = 0;
	while (j2 <= N) {
		A2 = AL[j2];
		same = 1;
		for (k=0; k<=M; ++k)
			if ((A1[k] == DASH) != (A2[k] == DASH)) same = 0;
		if (same == 0) {
			printf("  l");
			for (k=0; k<=M; ++k)
			   printf(" %d",POS[k]+1);
			for (k=0; k<=M; ++k)
			   if (A1[k] == DASH) printf(" %d",POS[k]);
			   else {
				POS[k] = POS[k] + (j2-j1);
				printf(" %d",POS[k]);
			   }
			printf(" %.1f\n",0.0);
			j1 = j2;
			A1 = AL[j1];
		}
		++j2;
	}
	/* Output the last block */
	printf("  l");
	for (k=0; k<=M; ++k)
	   printf(" %d",POS[k]+1);
	for (k=0; k<=M; ++k)
	   if (A1[k] == DASH) printf(" %d",POS[k]);
	   else {
		POS[k] = POS[k] + (j2-j1);
		printf(" %d",POS[k]);
	   }
	printf(" %.1f\n",0.0);
	printf("}\n");
}

char **AL1, **AL2;
int M1,N1,M2,N2;
int L1,L2;

/* MALIGN - initial routine for calling malign */
int MALIGN(char **IAL1, int IM1, int IN1, char **IAL2, int IM2, int IN2, int *ALS, int *SL, int IL1, int IL2)
{
	char *pal, *cal, ci, cj, gci, gcj;
	int i,j,I,J;
	int t;
	char *ckalloc();
	char *e1, *e2;
	char *he1, *ve2;
	int c;

	AL1 = IAL1;
	AL2 = IAL2;
	M1 = IM1;
	N1 = IN1;
	M2 = IM2;
	N2 = IN2;
	L1 = IL1;
	L2 = IL2;
	j = (N1+1) * sizeof(int);
	C1 = (int *) ckalloc(j);
	j = (N2+1) * sizeof(int);
	C2 = (int *) ckalloc(j);
	HH = (int *) ckalloc(j);
	VV = (int *) ckalloc(j);
	DD = (int *) ckalloc(j);
	RH = (int *) ckalloc(j);
	RV = (int *) ckalloc(j);
	RD = (int *) ckalloc(j);

	sl = 0;
	sapp = ALS;
	last = 0;

	E1 = (char **) ckalloc((N1+1)*sizeof(char *));
	E1[0] = (char *) ckalloc((M1+1)*(N1+1));
	HE1 = (char **) ckalloc((N1+1)*sizeof(char *));
	HE1[0] = (char *) ckalloc((M1+1)*(N1+1));
	for (i=1; i<=N1; ++i) {
		E1[i] = E1[i-1]+(M1+1);
		HE1[i] = HE1[i-1]+(M1+1);
	}
	E2 = (char **) ckalloc((N2+1)*sizeof(char *));
	E2[0] = (char *) ckalloc((M2+1)*(N2+1));
	VE2 = (char **) ckalloc((N2+1)*sizeof(char *));
	VE2[0] = (char *) ckalloc((M2+1)*(N2+1));
	for (i=1; i<=N2; ++i) {
		E2[i] = E2[i-1]+(M2+1);
		VE2[i] = VE2[i-1]+(M2+1);
	}

	/* Compute E1 and E2 */
	for (j=0; j<=N1; ++j) {
		e1 = E1[j];
		he1 = HE1[j];
		for (i=0; i<=M1; ++i) he1[i] = e1[i] = DASH;
	}
	for (i=0; i<=M1; ++i) {
		j = 0;
		while (j<=N1 && AL1[j][i]==DASH) {
			HE1[j][i] = E1[j][i] = LSTAR;
			++j;
		}
		j = N1;
		while (j>=1 && AL1[j][i]==DASH) {
			HE1[j][i] = E1[j][i] = RSTAR;
			--j;
		}
		HE1[j][i] = RSTAR;
	}

	for (j=0; j<=N2; ++j) {
		e2 = E2[j];
		ve2 = VE2[j];
		for (i=0; i<=M2; ++i) ve2[i] = e2[i] = DASH;
	}
	for (i=0; i<=M2; ++i) {
		j = 0;
		while (j<=N2 && AL2[j][i]==DASH) {
			VE2[j][i] = E2[j][i] = LSTAR;
			++j;
		}
		j = N2;
		while (j>=1 && AL2[j][i]==DASH) {
			VE2[j][i] = E2[j][i] = RSTAR;
			--j;
		}
		VE2[j][i] = RSTAR;
	}


	/* Compute projection cost within AL1 and AL2 */
	C1[0] = 0;
	pal = AL1[0];
	for (i=1; i<=N1; ++i) {
		t = 0;
		e1 = E1[i];
		cal = AL1[i];
		for (I=0; I<M1; ++I) {
			ci = ET(cal[I],e1[I]);
			gci = GV(ci);
			for (J=I+1; J<=M1; ++J) {
				cj = ET(cal[J],e1[J]);
				gcj = GV(cj);
				t += PD(L1+I,L1+J,ci,cj);
				t -= PT[L1+I][L1+J][pal[I]!=DASH][pal[J]!=DASH][gci][gcj];
			}
		}
		C1[i] = t;
		pal = cal;
	}

	C2[0] = 0;
	pal = AL2[0];
	for (j=1; j<=N2; ++j) {
		t = 0;
		e2 = E2[j];
		cal = AL2[j];
		for (I=0; I<M2; ++I) {
			ci = ET(cal[I],e2[I]);
			gci = GV(ci);
			for (J=I+1; J<=M2; ++J) {
				cj = ET(cal[J],e2[J]);
				gcj = GV(cj);
				t += PD(L2+I,L2+J,ci,cj);
				t -= PT[L2+I][L2+J][pal[I]!=DASH][pal[J]!=DASH][gci][gcj];
			}
		}
		C2[j] = t;
		pal = cal;
	}

	c = malign(0,0,0,0,0,N1,N2,0,0,0);

	*SL = sl;
	free(C1); free(C2); free(VV); free(DD); free(HH);
	free(RH); free(RV); free(RD);
	free(E1[0]); free(E1);
	free(E2[0]); free(E2);
	return c;
}

/* malign - return an optimal alignment of two alignments */
int malign(int i1, int j1, int hh, int vv, int dd, int i2, int j2, int rh, int rv, int rd)
{
	int i, j, I, J;
	int h, v, d, t, c;
	int hg, vg, dg;
	char *cal1, *pal1, *cal2, *pal2;
	char ci, cj, pi, pj, ei, ej;
	int midi, midj;
	int midc, score;
	char *e1, *e2, *he1, *ve2;

	if (j2-j1 <= 0) {
		if (i2 > i1) DEL(i2-i1)
		return;
	}
	if (i2-i1 <= 0) {
		INS(j2-j1)
		return;
	}

	/* rounded up; for the sake of the two-row case */
	midi = (i1+i2+1)/2;

	/* Forward computation */
	/* Initialize the first row */
	HH[j1] = hh;
	VV[j1] = vv;
	DD[j1] = dd;
	cal1 = AL1[i1];
	pal2 = AL2[j1];
	he1 = HE1[i1];
	for (j=j1+1; j<=j2; ++j) {
		cal2 = AL2[j];
		e2 = E2[j];
		hg = HH[j-1];
		vg = VV[j-1];
		dg = DD[j-1];
		c = C2[j];
		for (I=0; I<=M2; ++I) {
			cj = ET(cal2[I],e2[I]);
			for (J=0; J<=M1; ++J) c += PD(L1+J,L2+I,he1[J],cj);
		}
		for (I=0; I<=M1; ++I) {
		   ci = cal1[I];
		   ei = GV(he1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = GV(cal2[J]);
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ei][cj];
			vg -= PT[L1+I][L2+J][ci!=DASH][0][ei][cj];
			dg -= PT[L1+I][L2+J][ci!=DASH][pj!=DASH][ei][cj];
		   }
		}
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
		   	ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		c += max3(hg,vg,dg);
		HH[j] = c;
		VV[j] = DD[j] = c+MININT;
		pal2 = cal2;
	   /*
	   printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i1,j,DD[j],HH[j],VV[j]);
	   */
	}

	/* Compute dynamic programming matrix */
	pal1 = AL1[i1];
	for (i=i1+1; i<=midi; ++i) {
	   cal1 = AL1[i];
	   e1 = E1[i];
	   he1 = HE1[i];

	   /* Initialize the first column */
	   ve2 = VE2[j1];
	   c = C1[i];
	   for (I=0; I<=M1; ++I) {
		ci = ET(cal1[I],e1[I]);
		for (J=0; J<=M2; ++J) c += PD(L1+I,L2+J,ci,ve2[J]);
	   }
	   hg = HH[j1];
	   vg = VV[j1];
	   dg = DD[j1];
	   cal2 = AL2[j1];
	   for (I=0; I<=M1; ++I) {
	      ci = GV(cal1[I]);
	      pi = pal1[I];
	      for (J=0; J<=M2; ++J) {
		   cj = cal2[J];
		   ej = GV(ve2[J]);
		   hg -= PT[L1+I][L2+J][0][cj!=DASH][ci][ej];
		   vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][ej];
		   dg -= PT[L1+I][L2+J][pi!=DASH][cj!=DASH][ci][ej];
	      }
	   }
	   /*
	   for (I=0; I<M1; ++I) {
	      ETGV(ci,cal1[I],e1[I])
	      pi = pal1[I];
	      for (J=I+1; J<=M1; ++J) {
	      	   ETGV(cj,cal1[J],e1[J])
		   pj = pal1[J];
		   hg -= PT[L1+I][L1+J][0][0][ci][cj];
		   vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
	      }
	   }
	   */
	   c += max3(hg,vg,dg);
	   v = c;
	   h = d = v+MININT;

	   /*
	   printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i,j1,d,h,v);
	   */

	   /* Compute the remaining columns */
	   pal2 = AL2[j1];
	   for (j=j1+1; j<=j2; ++j) {
		cal2 = AL2[j];
		e2 = E2[j];
		ve2 = VE2[j];

		/* Compute DD */
		c = C1[i] + C2[j];
		for (I=0; I<=M1; ++I) {
		   ci = ET(cal1[I],e1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = ET(cal2[J],e2[J]);
			c += PD(L1+I,L2+J,ci,cj);
		   }
		}
		hg = HH[j-1];
		vg = VV[j-1];
		dg = DD[j-1];
		for (I=0; I<=M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ci][cj];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][cj];
			dg -= PT[L1+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		/*
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
			ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		c += max3(hg,vg,dg);
		hg = HH[j-1] = h;
		vg = VV[j-1] = v;
		dg = DD[j-1] = d;
		d = c;

		/* Compute HH */
		c = C2[j];
		for (I=0; I<=M2; ++I) {
			cj = ET(cal2[I],e2[I]);
			for (J=0; J<=M1; ++J) c += PD(L1+J,L2+I,he1[J],cj);
		}
		for (I=0; I<=M1; ++I) {
		   ci = cal1[I];
		   ei = GV(he1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = GV(cal2[J]);
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ei][cj];
			vg -= PT[L1+I][L2+J][ci!=DASH][0][ei][cj];
			dg -= PT[L1+I][L2+J][ci!=DASH][pj!=DASH][ei][cj];
		   }
		}
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		c += max3(hg,vg,dg);
		h = c;

		/* Compute VV */
	   	c = C1[i];
	   	for (I=0; I<=M1; ++I) {
			ci = ET(cal1[I],e1[I]);
			for (J=0; J<=M2; ++J) c += PD(L1+I,L2+J,ci,ve2[J]);
	   	}
		hg = HH[j];
		vg = VV[j];
		dg = DD[j];
		for (I=0; I<=M1; ++I) {
		   ci = GV(cal1[I]);
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			cj = cal2[J];
			ej = GV(ve2[J]);
			hg -= PT[L1+I][L2+J][0][cj!=DASH][ci][ej];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][ej];
			dg -= PT[L1+I][L2+J][pi!=DASH][cj!=DASH][ci][ej];
		   }
		}
		/*
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
			ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		c += max3(hg,vg,dg);
		v = c;
		/*
		printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i,j,d,h,v);
		*/
		pal2 = cal2;
	   }
	   pal1 = cal1;
	   HH[j2] = h;
	   VV[j2] = v;
	   DD[j2] = d;
	}
	/*
	printf("malign: d=%d, v=%d, h=%d\n",d,v,h);
	printf("score=%d\n",max3(d,v,h));
	fflush(stdout);
	*/

	/* Backward computation */
	/* Initialize the last row */
	h = RH[j2] = rh;
	RV[j2] = rv;
	RD[j2] = rd;
	cal2 = AL2[j2];
	cal1 = AL1[i2];
	he1 = HE1[i2];
	for (j=j2-1; j>=j1; --j) {
		pal2 = AL2[j];
		c = C2[j+1];
		e2 = E2[j+1];
		for (I=0; I<=M2; ++I) {
			cj = ET(cal2[I],e2[I]);
			for (J=0; J<=M1; ++J) c += PD(L2+I,L1+J,cj,he1[J]);
		}
		d = v = h += c;
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			h -= PT[L2+I][L2+J][pi!=DASH][pal2[J]!=DASH][ci][cj];
		   }
		}
		*/
		for (I=0; I<=M2; ++I) {
		   ci = GV(cal2[I]);
		   pi = pal2[I];
		   for (J=0; J<=M1; ++J) 
		   	h -= PT[L1+J][L2+I][0][pi!=DASH][GV(he1[J])][ci];
		}
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			d -= PT[L2+I][L2+J][pi!=DASH][pal2[J]!=DASH][ci][cj];
		   }
		}
		*/
		for (J=0; J<=M2; ++J) {
		   cj = GV(cal2[J]);
		   pj = pal2[J];
		   for (I=0; I<=M1; ++I)
			d -= PT[L1+I][L2+J][cal1[I]!=DASH][pj!=DASH][GV(he1[I])][cj];
		}
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			v -= PT[L2+I][L2+J][0][0][ci][cj];
		   }
		}
		*/
		for (J=0; J<=M2; ++J) {
		   cj = GV(cal2[J]);
		   for (I=0; I<=M1; ++I)
			v -= PT[L1+I][L2+J][cal1[I]!=DASH][0][GV(he1[I])][cj];
		}
		RH[j] = h;
		RD[j] = d;
		RV[j] = v;
		cal2 = pal2;

		/*
		printf("i=%d, j=%d, d= %d, v=%d, h=%d\n",i2,j,d,v,h);
		*/
	}

	/* Compute dynamic programming matrix */
	cal1 = AL1[i2];
	for (i=i2-1; i>=midi; --i) {
	   pal1 = AL1[i];
	   cal2 = AL2[j2];

	   /* Initialize the last column */
	   c = C1[i+1];
	   e1 = E1[i+1];
	   ve2 = VE2[j2];
	   for (I=0; I<=M1; ++I) {
		ci = ET(cal1[I],e1[I]);
		for (J=0; J<=M2; ++J) c += PD(L1+I,L2+J,ci,ve2[J]);
	   }
	   d = h = v = RV[j2] + c;
	   /*
	   for (I=0; I<M1; ++I) {
		ETGV(ci,cal1[I],e1[I])
		pi = pal1[I];
		for (J=I+1; J<=M1; ++J) {
		   ETGV(cj,cal1[J],e1[J])
		   v -= PT[L1+I][L1+J][pi!=DASH][pal1[J]!=DASH][ci][cj];
		}
	   }
	   */
	   for (I=0; I<=M1; ++I) {
		ci = GV(cal1[I]);
		pi = pal1[I];
		for (J=0; J<=M2; ++J)
			v -= PT[L1+I][L2+J][pi!=DASH][0][ci][GV(ve2[J])];
	   }
	   /*
	   for (I=0; I<M1; ++I) {
		ETGV(ci,cal1[I],e1[I])
	   	pi = pal1[I];
	   	for (J=I+1; J<=M1; ++J) {
		   ETGV(cj,cal1[J],e1[J])
		   d -= PT[L1+I][L1+J][pi!=DASH][pal1[J]!=DASH][ci][cj];
		}
	   }
	   */
	   for (I=0; I<=M1; ++I) {
	   	ci = GV(cal1[I]);
	   	pi = pal1[I];
	   	for (J=0; J<=M2; ++J)
		   d -= PT[L1+I][L2+J][pi!=DASH][cal2[J]!=DASH][ci][GV(ve2[J])];
	   }
	   /*
	   for (I=0; I<M1; ++I) {
		ETGV(ci,cal1[I],e1[I])
	   	for (J=I+1; J<=M1; ++J) {
		   ETGV(cj,cal1[J],e1[J])
		   h -= PT[L1+I][L1+J][0][0][ci][cj];
		}
	   }
	   */
	   for (I=0; I<=M1; ++I) {
	   	ci = GV(cal1[I]);
	   	for (J=0; J<=M2; ++J)
		   h -= PT[L1+I][L2+J][0][cal2[J]!=DASH][ci][GV(ve2[J])];
	   }

	   /*
	   printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i,j2,d,h,v);
	   */

	   /* Compute the remaining columns */
	   e1 = E1[i+1];
	   he1 = HE1[i];
	   for (j=j2-1; j>=j1; --j) {
		pal2 = AL2[j];
		e2 = E2[j+1];
		ve2 = VE2[j];

		/* Compute RD, RH, RV affected by the diagonal edge */
		c = C1[i+1] + C2[j+1];
		for (I=0; I<=M1; ++I) {
		   ci = ET(cal1[I],e1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = ET(cal2[J],e2[J]);
			c += PD(L1+I,L2+J,ci,cj);
		   }
		}
		hg = vg = dg = RD[j+1];
		RH[j+1] = h;
		RV[j+1] = v;
		RD[j+1] = d;
		for (I=0; I<=M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ci][cj];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][cj];
			dg -= PT[L1+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		/*
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
			ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		h = c+hg;
		v = c+vg;
		d = c+dg;

		/* Compute RD, RH, RV affected by the horizontal edge */
		c = C2[j+1];
		for (I=0; I<=M2; ++I) {
			cj = ET(cal2[I],e2[I]);
			for (J=0; J<=M1; ++J) c += PD(L1+J,L2+I,he1[J],cj);
		}
		hg = vg = dg = RH[j+1];
		for (I=0; I<=M1; ++I) {
		   pi = pal1[I];
		   ei = GV(he1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = GV(cal2[J]);
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ei][cj];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ei][cj];
			dg -= PT[L1+I][L2+J][pi!=DASH][pj!=DASH][ei][cj];
		   }
		}
		/*
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		if (c+hg > h) h = c+hg;
		if (c+vg > v) v = c+vg;
		if (c+dg > d) d = c+dg;

		/* Compute RD, RH, RV affected by the vertical edge */
	   	c = C1[i+1];
	   	for (I=0; I<=M1; ++I) {
			ci = ET(cal1[I],e1[I]);;
			for (J=0; J<=M2; ++J) c += PD(L1+I,L2+J,ci,ve2[J]);
	   	}
		hg = dg = vg = RV[j];
		for (I=0; I<=M1; ++I) {
		   ci = GV(cal1[I]);
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			ej = GV(ve2[J]);
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ci][ej];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][ej];
			dg -= PT[L1+I][L2+J][pi!=DASH][pj!=DASH][ci][ej];
		   }
		}
		/*
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
			ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		*/
		if (c+hg > h) h = c+hg;
		if (c+vg > v) v = c+vg;
		if (c+dg > d) d = c+dg;

		/*
		printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i,j,d,h,v);
		*/

		cal2 = pal2;
	   }
	   cal1 = pal1;
	   RH[j1] = h;
	   RV[j1] = v;
	   RD[j1] = d;
	}
	/*
	printf("malign: d=%d, v=%d, h=%d\n",d,v,h);
	printf("score=%d\n",max3(d,v,h));
	*/

	/* Determine optimal midpoint */
	midc = max(VV[j1]+RV[j1],DD[j1]+RD[j1]);
	midj = j1;
	for (j=j1+1; j<=j2; ++j) {
	   if ((t = max(VV[j]+RV[j],DD[j]+RD[j])) > midc) {
	 	midc = t;
		midj = j;
	   }
	}

/*
	printf("i1=%d, i2=%d, j1=%d, j2=%d, midscore=%d, midj=%d\n",i1,i2,j1,j2,midc,midj);
*/

	if (i2-i1 <= 1) { /* two-row case */
		if (midc == (VV[midj]+RV[midj])) { /* vertical edge */
			if (midj > j1) INS(midj-j1)
			DEL(1)
			if (midj < j2) INS(j2-midj)
		} else { /* diagonal edge */
			if (midj-1 > j1) INS(midj-1-j1)
			REP
			if (midj < j2) INS(j2-midj)
		}
	} else {
		/* Divide into two subproblems */
		/* save overlapped variables */
		h = HH[midj];
		v = VV[midj];
		d = DD[midj];
		if (midc == (VV[midj]+RV[midj])) { /* vertical edge */
			malign(i1,j1,hh,vv,dd,midi,midj,MININT,RV[midj],MININT);
			malign(midi,midj,MININT,v,MININT,i2,j2,rh,rv,rd);
		} else { /* diagonal edge */
			malign(i1,j1,hh,vv,dd,midi,midj,MININT,MININT,RD[midj]);
			malign(midi,midj,MININT,MININT,d,i2,j2,rh,rv,rd);
		}
	}
	fflush(stdout);

	return midc;
}

/* Mscore - return the optimum score of two alignments alignment */
int Mscore(char **AL1, int M1, int N1, char **AL2, int M2, int N2, int L1, int L2)
{
	int i, j, I, J;
	int h, v, d, t, c;
	int hg, vg, dg;
	char *cal, *pal;
	char *cal1, *pal1, *cal2, *pal2;
	char ci, cj, pi, pj, ei, ej;
	char *ckalloc();
	char *e1, *e2;
	char *he1, *ve2;

	j = (N1+1) * sizeof(int);
	C1 = (int *) ckalloc(j);
	j = (N2+1) * sizeof(int);
	C2 = (int *) ckalloc(j);
	HH = (int *) ckalloc(j);
	VV = (int *) ckalloc(j);
	DD = (int *) ckalloc(j);

	E1 = (char **) ckalloc((N1+1)*sizeof(char *));
	E1[0] = (char *) ckalloc((M1+1)*(N1+1));
	HE1 = (char **) ckalloc((N1+1)*sizeof(char *));
	HE1[0] = (char *) ckalloc((M1+1)*(N1+1));
	for (i=1; i<=N1; ++i) {
		E1[i] = E1[i-1]+(M1+1);
		HE1[i] = HE1[i-1]+(M1+1);
	}
	E2 = (char **) ckalloc((N2+1)*sizeof(char *));
	E2[0] = (char *) ckalloc((M2+1)*(N2+1));
	VE2 = (char **) ckalloc((N2+1)*sizeof(char *));
	VE2[0] = (char *) ckalloc((M2+1)*(N2+1));
	for (i=1; i<=N2; ++i) {
		E2[i] = E2[i-1]+(M2+1);
		VE2[i] = VE2[i-1]+(M2+1);
	}

	/* Compute E1 and E2 */
	for (j=0; j<=N1; ++j) {
		e1 = E1[j];
		he1 = HE1[j];
		for (i=0; i<=M1; ++i) he1[i] = e1[i] = DASH;
	}
	for (i=0; i<=M1; ++i) {
		j = 0;
		while (j<=N1 && AL1[j][i]==DASH) {
			HE1[j][i] = E1[j][i] = LSTAR;
			++j;
		}
		j = N1;
		while (j>=1 && AL1[j][i]==DASH) {
			HE1[j][i] = E1[j][i] = RSTAR;
			--j;
		}
		HE1[j][i] = RSTAR;
	}

	for (j=0; j<=N2; ++j) {
		e2 = E2[j];
		ve2 = VE2[j];
		for (i=0; i<=M2; ++i) ve2[i] = e2[i] = DASH;
	}
	for (i=0; i<=M2; ++i) {
		j = 0;
		while (j<=N2 && AL2[j][i]==DASH) {
			VE2[j][i] = E2[j][i] = LSTAR;
			++j;
		}
		j = N2;
		while (j>=1 && AL2[j][i]==DASH) {
			VE2[j][i] = E2[j][i] = RSTAR;
			--j;
		}
		VE2[j][i] = RSTAR;
	}

	/* Compute projection cost within AL1 and AL2 */
	for (i=0; i<=N1; ++i) {
		t = 0;
		cal = AL1[i];
		e1 = E1[i];
		for (I=0; I<M1; ++I) {
			ci = ET(cal[I],e1[I]);
			for (J=I+1; J<=M1; ++J) {
				cj = ET(cal[J],e1[J]);
				t += PD(L1+I,L1+J,ci,cj);
			}
		}
		C1[i] = t;
	}
	for (j=0; j<=N2; ++j) {
		t = 0;
		cal = AL2[j];
		e2 = E2[j];
		for (I=0; I<M2; ++I) {
			ci = ET(cal[I],e2[I]);
			for (J=I+1; J<=M2; ++J) {
				cj = ET(cal[J],e2[J]);
				t += PD(L2+I,L2+J,ci,cj);
			}
		}
		C2[j] = t;
	}

	/* Initialize the first row */
	HH[0] = h = VV[0] = DD[0] = 0;
	e1 = E1[0];
	he1 = HE1[0];
	pal = AL2[0];
	for (j=1; j<=N2; ++j) {
		cal = AL2[j];
		e2 = E2[j];
		c = C2[j];
		for (I=0; I<=M2; ++I) {
			ci = ET(cal[I],e2[I]);
			for (J=0; J<=M1; ++J) {
				c += PD(L2+I,L1+J,he1[J],ci);
			}
		}
		h += c;
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal[I],e2[I])
		   pi = pal[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal[J],e2[J])
			h -= PT[L2+I][L2+J][pi!=DASH][pal[J]!=DASH][ci][cj];
		   }
		}
		for (I=0; I<=M2; ++I) {
		   ci = GV(cal[I]);
		   pi = pal[I];
		   for (J=0; J<=M1; ++J)
			h -= PT[L2+I][L1+J][0][pi!=DASH][GV(he1[J])][ci];
		}
		pal = cal;
		HH[j] = h;
		VV[j] = DD[j] = h+MININT;
		/*
		printf("i=0, j=%d, h=%d\n",j,h);
		*/
	}

	/* Compute dynamic programming matrix */
	pal1 = AL1[0];
	for (i=1; i<=N1; ++i) {
	   cal1 = AL1[i];
	   e1 = E1[i];
	   he1 = HE1[i];

	   /* Initialize the first column */
	   e2 = E2[0];
	   ve2 = VE2[0];
	   v = C1[i];
	   for (I=0; I<=M1; ++I) {
		ci = ET(cal1[I],e1[I]);
		for (J=0; J<=M2; ++J) v += PD(L1+I,L2+J,ci,ve2[J]);
	   }
	   v += VV[0];
	   for (I=0; I<M1; ++I) {
		ETGV(ci,cal1[I],e1[I])
		pi = pal1[I];
		for (J=I+1; J<=M1; ++J) {
		   ETGV(cj,cal1[J],e1[J])
		   v -= PT[L1+I][L1+J][pi!=DASH][pal1[J]!=DASH][ci][cj];
		}
	   }
	   for (I=0; I<=M1; ++I) {
		ci = GV(cal1[I]);
		pi = pal1[I];
		for (J=0; J<=M2; ++J)
		   v -= PT[L1+I][L2+J][pi!=DASH][0][ci][GV(ve2[J])];
	   }
	   h = d = v+MININT;

	   /*
	   printf("i=%d, j=0, v=%d\n",i,v);
	   */
	   
	   /* Compute the remaining columns */
	   pal2 = AL2[0];
	   for (j=1; j<=N2; ++j) {
		cal2 = AL2[j];
		e2 = E2[j];
		ve2 = VE2[j];

		/* Compute DD */
		c = C1[i] + C2[j];
		for (I=0; I<=M1; ++I) {
		   ci = ET(cal1[I],e1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = ET(cal2[J],e2[J]);
			c += PD(L1+I,L2+J,ci,cj);
		   }
		}
		hg = HH[j-1];
		vg = VV[j-1];
		dg = DD[j-1];
		for (I=0; I<=M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ci][cj];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][cj];
			dg -= PT[L1+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
		   	ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
		   	ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		c += max3(hg,vg,dg);
		hg = HH[j-1] = h;
		vg = VV[j-1] = v;
		dg = DD[j-1] = d;
		d = c;

		/* Compute HH */
		c = C2[j];
		for (I=0; I<=M2; ++I) {
			cj = ET(cal2[I],e2[I]);
			for (J=0; J<=M1; ++J) c += PD(L2+I,L1+J,he1[J],cj);
		}
		for (I=0; I<=M1; ++I) {
		   ci = cal1[I];
		   ei = GV(he1[I]);
		   for (J=0; J<=M2; ++J) {
			cj = GV(cal2[J]);
			pj = pal2[J];
			hg -= PT[L1+I][L2+J][0][pj!=DASH][ei][cj];
			vg -= PT[L1+I][L2+J][ci!=DASH][0][ei][cj];
			dg -= PT[L1+I][L2+J][ci!=DASH][pj!=DASH][ei][cj];
		   }
		}
		for (I=0; I<M2; ++I) {
		   ETGV(ci,cal2[I],e2[I])
		   pi = pal2[I];
		   for (J=I+1; J<=M2; ++J) {
			ETGV(cj,cal2[J],e2[J])
			pj = pal2[J];
			hg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
			vg -= PT[L2+I][L2+J][0][0][ci][cj];
			dg -= PT[L2+I][L2+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		c += max3(hg,vg,dg);
		h = c;

		/* Compute VV */
	   	c = C1[i];
	   	for (I=0; I<=M1; ++I) {
			ci = ET(cal1[I],e1[I]);
			for (J=0; J<=M2; ++J) c += PD(L1+I,L2+J,ci,ve2[J]);
	   	}
		hg = HH[j];
		vg = VV[j];
		dg = DD[j];
		for (I=0; I<=M1; ++I) {
		   ci = GV(cal1[I]);
		   pi = pal1[I];
		   for (J=0; J<=M2; ++J) {
			cj = cal2[J];
			ej = GV(ve2[J]);
			hg -= PT[L1+I][L2+J][0][cj!=DASH][ci][ej];
			vg -= PT[L1+I][L2+J][pi!=DASH][0][ci][ej];
			dg -= PT[L1+I][L2+J][pi!=DASH][cj!=DASH][ci][ej];
		   }
		}
		for (I=0; I<M1; ++I) {
		   ETGV(ci,cal1[I],e1[I])
		   pi = pal1[I];
		   for (J=I+1; J<=M1; ++J) {
		   	ETGV(cj,cal1[J],e1[J])
			pj = pal1[J];
			hg -= PT[L1+I][L1+J][0][0][ci][cj];
			vg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
			dg -= PT[L1+I][L1+J][pi!=DASH][pj!=DASH][ci][cj];
		   }
		}
		c += max3(hg,vg,dg);
		v = c;
		/*
		printf("i=%d, j=%d, d=%d, h=%d, v=%d\n",i,j,d,h,v);
		*/
		pal2 = cal2;
	   }
	   pal1 = cal1;
	   HH[N2] = h;
	   VV[N2] = v;
	   DD[N2] = d;
	}
	free(C1); free(C2); free(VV); free(DD); free(HH);
	free(E1[0]); free(E1);
	free(E2[0]); free(E2);
	free(HE1[0]); free(HE1);
	free(VE2[0]); free(VE2);
	printf("d=%d, v=%d, h=%d\n",d,v,h);
	printf("score=%d\n",max3(d,v,h));
	fflush(stdout);
	return max3(d,v,h);
}

/* align_score - return the score of the alignment AL */
int align_score(char **AL, int M, int N)
{
	int I, J, i, j;
	int score;
	char *cal, *pal;
	char ci, cj, pi, ti, tj;
	char *e1;
	char *ckalloc();

	E1 = (char **) ckalloc((N+1)*sizeof(char *));
	E1[0] = (char *) ckalloc((M+1)*(N+1));
	for (i=1; i<=N; ++i) {
		E1[i] = E1[i-1]+(M+1);
	}
	for (j=0; j<=N; ++j) {
		e1 = E1[j];
		for (i=0; i<=M; ++i) e1[i] = DASH;
	}
	for (i=0; i<=M; ++i) {
		j = 0;
		while (j<=N && AL[j][i]==DASH) {
			E1[j][i] = LSTAR;
			++j;
		}
		j = N;
		while (j>=1 && AL[j][i]==DASH) {
			E1[j][i] = RSTAR;
			--j;
		}
	}

	score = 0;
	pal = AL[0];
	for (j=1; j<=N; ++j) {
		e1 = E1[j];
		cal = AL[j];
		for (I=0; I<M; ++I) {
			ci = cal[I];
			ti = ET(ci,e1[I]);
			ETGV(ci,cal[I],e1[I])
			pi = pal[I];
		   	for (J=I+1; J<=M; ++J) {
			   tj = ET(cal[J],e1[J]);
			   ETGV(cj,cal[J],e1[J])
			   score += PD(I+1,J+1,ti,tj);
			   score -= PT[I+1][J+1][pi!=DASH][pal[J]!=DASH][ci][cj];
		   	}
		}
		pal = cal;
	}
	free(E1[0]); free(E1);
	return score;
}

/* max3 :   three argument maximum */

int max3 (int x, int y, int z)
{
	if (x>y)
		if (x>z) return x;
		else return z;
	else
		if (y>z) return y;
		else return z;
}