/* A software package to get pairwise alignments
	and generate projected regions */
#include <stdio.h>
#include <ctype.h>

#define max(x,y)  ((x) >= (y) ? (x) : (y))
#define min(x,y)  ((x) <= (y) ? (x) : (y))

#define NUMBER 100		 /* maximum sequence piece number */
#define TRUE	1
#define FALSE	0
#define MAXBLOCKS 300

struct seq_struct {
	char *name;		/* sequence name */
	int from;		/* start position */
	int to;			/* end position */
	int len;		/* length of the sequence */
	int bfrom;		/* projected base start */
	int bto;		/* projected base end */
	int block[MAXBLOCKS][4];/* aligned blocks */
	int bn;			/* block number */
};
struct seq_struct seq[NUMBER+1];
int seq_no;
int Width, *From, *To;
char **Sname;
extern int     *Seq_Script[];  /* aligned sequence script */
extern int     Entry_No[];     /* Number of script entries */

/* get_all - get all alignments */
int get_all(char *Rname, int width, char *ISname[], int *IFrom, int *ITo)
{
	FILE *fp, *lav_fp, *ckopen();
	char fname[30];
	char *sname1, *sname2, *get_seqname();
	int from1, to1, from2, to2, i, bn;

	Width = width;
	From = IFrom;
	To = ITo;
	Sname = ISname;
	fp = ckopen(Rname,"r");
	seq_no = 0;
	while (fscanf(fp,"%s",&fname) != EOF) {
	   lav_fp = ckopen(fname,"r");
	   lav_read_init(lav_fp);
	   sname1 = get_seqname(&from1, &to1, 0);
	   sname2 = get_seqname(&from2, &to2, 0);;
	   while (get_alignment(lav_fp)) {
		bn = seq[seq_no].bn;
		seq[seq_no].name = sname2;
		seq[seq_no].from = seq[seq_no].from + from2 - 1;
		seq[seq_no].to = seq[seq_no].to + from2 - 1;
		seq[seq_no].len = seq[seq_no].block[bn-1][3];

		/*
		printf("*******************************\n");
		printf("seq_no: %d\n",seq_no);
		printf("name: %s\n",seq[seq_no].name);
		printf("from: %d\n",seq[seq_no].from);
		printf("to: %d\n",seq[seq_no].to);
		printf("len: %d\n",seq[seq_no].len);
		printf("bfrom: %d\n",seq[seq_no].bfrom);
		printf("bto: %d\n",seq[seq_no].bto);
		printf("bn: %d\n",bn);
		for (i=0; i<bn; ++i)
			printf("%d %d %d %d\n",seq[seq_no].block[i][0],
			   seq[seq_no].block[i][1], seq[seq_no].block[i][2],
			   seq[seq_no].block[i][3]);
		printf("*******************************\n");
		*/

		seq_no++;
	   }
	   fclose(lav_fp);
	}
	fclose(fp);
}

/* get_alignment - get one alignment */
int get_alignment(FILE *lav_fp)
{
	char lflag, c, buf[100], *s, *strtok();
	int bn;

	lflag = 0;
	while (fgets(buf, 100, lav_fp)) {
		if (buf[0] == 'a') {
			lflag = 1;
			break;
		}
	}
	if (lflag) {
	   fgets(buf, 100, lav_fp);
	   fgets(buf, 100, lav_fp);
	   sscanf(buf, "%s %d %d", &c, &seq[seq_no].bfrom, &seq[seq_no].from);
	   fgets(buf, 100, lav_fp);
	   sscanf(buf, "%s %d %d", &c, &seq[seq_no].bto, &seq[seq_no].to);
	   bn = 0;
	   while (fgets(buf, 100, lav_fp)) {
		if (sscanf(buf, "%s %d %d %d %d",&c, &seq[seq_no].block[bn][0],
			   &seq[seq_no].block[bn][1], &seq[seq_no].block[bn][2],
			   &seq[seq_no].block[bn][3]) == 5) {
			seq[seq_no].block[bn][1] = seq[seq_no].block[bn][1]
						   - seq[seq_no].from + 1;
			seq[seq_no].block[bn][3] = seq[seq_no].block[bn][3]
						   - seq[seq_no].from + 1;
			bn++;
		} else break;
	   }
	   seq[seq_no].bn = bn;
	   return 1;
	} else return 0;
}

/* LR - compute the region defined by pairwise regions */
void LR(int M1, int N1, int L1, int M2, int N2, int L2, int *LB, int *RB, int *lp, int *rp)
{
	int i, j, m1, m2, t, i1, j1, k;
	int SL1, SL2;
	char *ckalloc();
	int ii, pi, pl, pr;
	char *LDEF, *RDEF;
	int *P_index1, *P_index2;
	int sn, sn1, sn2;
	int *LL, *RR;
	int bfrom, bto, len;
	int lmost, rmost;

	P_index1 = (int *) ckalloc((N1+1)*sizeof(int));
	P_index2 = (int *) ckalloc((N2+1)*sizeof(int));
	LDEF = (char *) ckalloc(N1+1);
	RDEF = (char *) ckalloc(N1+1);
	lmost = N1;
	rmost = 0;

	for (i=0; i<=N1; ++i) {
		LB[i] = 0;
		RB[i] = N2;
		LDEF[i] = RDEF[i] = 0;
	}
	/* Projection omitted */
	for (m1=0; m1<=M1; ++m1) {
		SL1 = Padded_index(L1+m1, P_index1);
		for (m2=0; m2<=M2; ++m2) {
		   SL2 = Padded_index(L2+m2, P_index2);
		   if (L1+m1 == 1) { /* L1+m1 is the principle sequence */
		      /* get the alignment between L1+m1 and L2+m2 */
		      /* transfer the alignment into boundary lines */
		      /* define projected region */
		      sn = get_sn(L2+m2);
		      if (P_index1[seq[sn].bfrom] < lmost) lmost = P_index1[seq[sn].bfrom];
		      if (P_index1[seq[sn].bto] > rmost) rmost = P_index1[seq[sn].bto];
		      for (i=0; i<seq[sn].bn; ++i) {
			 if (i>0 && seq[sn].block[i][0] > seq[sn].block[i-1][2]+1) {	/* horizontal line */
			    k = P_index2[seq[sn].block[i-1][3]];
			    for (j=seq[sn].block[i-1][2]+1; j<seq[sn].block[i][0]; ++j) {
				j1 = P_index1[j];
				LB[j1] = max(LB[j1], k);
				RB[j1] = min(RB[j1], k);
				LDEF[j1] = RDEF[j1] = 1;
			    }
			 }

			 if (i>0 && seq[sn].block[i][1] > seq[sn].block[i-1][3]+1) {	/* vertical line */
			    k = P_index2[seq[sn].block[i][1]];
			    j1 = P_index1[seq[sn].block[i-1][2]];
			    /* update the upper bound */
			    RB[j1] = max(RB[j1], k);
			 }

			 /* print out diagonal lines */
			 for (j=seq[sn].block[i][0], i1=seq[sn].block[i][1];
			      j<=seq[sn].block[i][2]; ++j, ++i1) {
			    k = P_index2[i1];
			    j1 = P_index1[j];
			    LB[j1] = max(LB[j1], k);
			    RB[j1] = min(RB[j1], k);
			    LDEF[j1] = RDEF[j1] = 1;
			 }
		      }
		   } else {
		      sn1 = get_sn(L1+m1);
		      sn2 = get_sn(L2+m2);

		      /* get the alignment between 1 and L2+m2 */
		      /* transfer the alignment into boundary lines */
		      bfrom = seq[sn2].bfrom;
		      bto = seq[sn2].bto;
		      j = bto - bfrom + 1;
		      LL = (int *) ckalloc(j*sizeof(int)) - bfrom;
		      RR = (int *) ckalloc(j*sizeof(int)) - bfrom;
		      for (i=0; i<seq[sn2].bn; ++i) {
			 if (i>0 && seq[sn2].block[i][0] > seq[sn2].block[i-1][2]+1) {	/* horizontal line */
			    k = P_index2[seq[sn2].block[i-1][3]];
			    for (j=seq[sn2].block[i-1][2]+1; j<seq[sn2].block[i][0]; ++j) {
				LL[j] = k;
				RR[j] = k;
			    }
			 }

			 if (i>0 && seq[sn2].block[i][1] > seq[sn2].block[i-1][3]+1) {	/* vertical line */
			    k = P_index2[seq[sn2].block[i][1]];
			    /* update the upper bound */
			    j = seq[sn2].block[i-1][2];
			    RR[j] = k;
			 }

			 /* print out diagonal lines */
			 for (j=seq[sn2].block[i][0], i1=seq[sn2].block[i][1];
			      j<=seq[sn2].block[i][2]; ++j, ++i1) {
			    k = P_index2[i1];
			    LL[j] = k;
			    RR[j] = k;
			 }
		      }

		      /* get the alignment between L1+m1 and 1 */
		      /* transfer the alignment into boundary lines */
		      /* define projected region */
		      for (i=0; i<seq[sn1].bn; ++i) {
			 if (i>0 && seq[sn1].block[i][1] > seq[sn1].block[i-1][3]+1) {	/* horizontal line */
			    k = seq[sn1].block[i-1][2];
			    for (j=seq[sn1].block[i-1][3]+1; j<seq[sn1].block[i][1]; ++j) {
				j1 = P_index1[j];
				if (k>=bfrom && k<=bto) {
				   LB[j1] = max(LB[j1],LL[k]);
				   RB[j1] = min(RB[j1],RR[k]);
				   LDEF[j1] = RDEF[j1] = 1;
				} else if (k<bfrom) {
				   RB[j1] = min(RB[j1],LL[bfrom]-1);
				   RDEF[j1] = 1;
				} else {
				   LB[j1] = max(LB[j1],RR[bto]);
				   LDEF[j1] = 1;
				}
			    }
			 }

			 if (i>0 && seq[sn1].block[i][0] > seq[sn1].block[i-1][2]+1) {	/* vertical line */
			    k = seq[sn1].block[i][0];
			    j1 = P_index1[seq[sn1].block[i-1][3]];
			    /* update the upper bound */
			    if (k>=bfrom && k<=bto) {
			       RB[j1] = RR[k];
			       RDEF[j1] = 1;
			    }
			 }

			 /* print out diagonal lines */
			 for (j=seq[sn1].block[i][1], k=seq[sn1].block[i][0];
			      j<=seq[sn1].block[i][3]; ++j, ++k) {
			    j1 = P_index1[j];
			    if (k>=bfrom && k<=bto) {
				   LB[j1] = max(LB[j1],LL[k]);
				   RB[j1] = min(RB[j1],RR[k]);
				   LDEF[j1] = RDEF[j1] = 1;
			    } else if (k<bfrom) {
				   RB[j1] = min(RB[j1],LL[bfrom]-1);
				   RDEF[j1] = 1;
			    } else {
				   LB[j1] = max(LB[j1],RR[bto]);
				   LDEF[j1] = 1;
			    }
			 }
		      }
		      free(LL+bfrom); free(RR+bfrom);
		   }

		   /*
		   if (RD[L1+m1][L2+m2]==1) {
		      SL2 = Padded_index(L2+m2, P_index2);
		      LL = PLB[L1+m1][L2+m2];
		      RR = PRB[L1+m1][L2+m2];
		      for (i=0; i<=SL1; ++i) {
			   ii = P_index1[i];
			   if (i==0) pi=-1;
			   else pi = P_index1[i-1];
			   for (j=pi+1; j<=ii; ++j) {
			      LB[j] = max(LB[j], P_index2[LL[i]]);
			      RB[j] = min(RB[j], P_index2[RR[i]]);
			   }
		      }
		   }
		   */
		}
	}

	/*
	printf("N1=%d,N2=%d\n",N1,N2);
	for (i=0; i<=N1; ++i) printf("i=%d, %d %d\n",i,LB[i],RB[i]);
	printf("end\n");
	fflush(stdout);
	*/

	/* Normalization */
	pl = 0;
	for (i=0; i<=N1; ++i)
		if (LDEF[i]) pl = LB[i];
		else LB[i] = pl;
	pr = N2;
	for (i=N1; i>=0; --i)
		if (RDEF[i]) pr = RB[i];
		else RB[i] = pr;
	for (i=N1-1; i>0; --i) LB[i] = min(LB[i],LB[i+1]);
	for (i=1; i<N1; ++i) RB[i] = max(RB[i-1],RB[i]);
	for (i=1; i<N1; ++i) 
		if (LB[i] > RB[i]) {
			t = (LB[i] + RB[i])/2;
			LB[i] = t;
			RB[i] = t;
		}
	lmost = max(0, lmost-Width);
	rmost = min(N1, rmost+Width);

	/*
	fprintf(stderr, "lmost=%d, rmost=%d\n",lmost,rmost);
	*/

	if (lmost > 0) pl = LB[lmost-1];
	else pl = 0;
	if (rmost < N1) pr = RB[rmost+1];
	else pr = N2;
	for (i=rmost; i>=lmost; --i) {
		ii = max(0,i-Width);
		if (LB[i]-Width > LB[ii])
			LB[i] = max(pl, LB[ii]);
		else LB[i] = max(pl,LB[i]-Width);
	}
	for (i=lmost; i<=rmost; ++i) {
		ii = min(N1,i+Width);
		if (RB[i]+Width < RB[ii])
			RB[i] = min(pr, RB[ii]);
		else RB[i] = min(pr,RB[i]+Width);
	}

	/*
	for (i=0; i<=N1; ++i) printf("i=%d, %d %d\n",i,LB[i],RB[i]);
	fflush(stdout);
	*/

	for (i=0; i<N1; ++i)
		if (RB[i] < LB[i+1]-1) {
			RB[i] = LB[i+1]-1;
		}
	LB[0] = 0;
	RB[N1] = N2;

	/*
	printf("N1=%d,N2=%d\n",N1,N2);
	for (i=0; i<=N1; ++i) printf("i=%d, %d %d\n",i,LB[i],RB[i]);
	fflush(stdout);
	*/

	/*
	for (i=1; i<=N1; ++i) {
		if (LB[i] < LB[i-1] || RB[i] < RB[i-1] || LB[i] > RB[i]) {
			printf("i-1=%d, %d %d\n",i-1,LB[i-1],RB[i-1]);
			printf("i=%d, %d %d\n",i,LB[i],RB[i]);
			exit(1);
		}
	}
	if (LB[0] != 0) {
			printf("i=%d, %d %d\n",0,LB[0],RB[0]);
			exit(1);
	}
	if (RB[N1] != N2) {
			printf("i=%d, %d %d\n",N1,LB[N1],RB[N1]);
			printf("N2=%d\n",N2);
			exit(1);
	}
	*/

	free(P_index1);
	free(P_index2);
	free(LDEF);
	free(RDEF);
	if (lmost > rmost) {
		lmost = 0;
		rmost = N1;
	}
	*lp = lmost;
	*rp = rmost;
}


/* Padded_index - compute the new indexes in the alignment */
int Padded_index(int m, int *P_index)
{
	int ci, i, j, k, c, *sp;

	P_index[0] = 0;

	ci = 1;
	j = 1;
	sp = Seq_Script[m];
	for (k=0; k<Entry_No[m]; ++k) {
		if ((c = sp[k]) > 0) {	/* non-gap */
			for (i=ci; i<ci+c; ++i) P_index[i] = j++;
			ci = ci+c;
		} else {	/* gap */
			j = j-c;
		}
	}
	return i;
}

/* get_sn - get the sequence index in seq structure */
int get_sn(int k)
{
	int i, from, to;
	char *sname;

	sname = Sname[k];
	from = From[k];
	to = To[k];
	for (i=0; i<seq_no; ++i)
	   if (strsame(sname,seq[i].name) && (from == seq[i].from)
	       && (to == seq[i].to)) return i;
	fatalf("The alignment file of seqeunce %d is not given!",k);
}