#include <stdio.h>
#define DIGIT 		10.0	/* all parameters are accurate to the tenth */
#define max(x,y)  ((x) >= (y) ? (x) : (y))
#define min(x,y)  ((x) <= (y) ? (x) : (y))

main(int argc, char *argv[])
{
	double atof();
	char *ckalloc();
	int K, Kcount;
	int *LB, *RB, cl, cr, pl, nr;
	int i, j, op;
	int fi, fj, fk, li, lj, lk;
	int ci, cj, ck;
	int bi, bj, ei, ej, t;
	double score, c;
	FILE *qp, *ckopen();
	int fscanf(), atoi();
	char *seq1, *seq2, sname1[50], sname2[50];
	int M, N;		/* lengths of sequences */
	int len1, len2;
	long ms, G, H, V[128][128];
	int *S, *sp;		/* conversion operations */
	int old_i, old_j, match_no;
	int pi, pj, ppi, ppj;
	int ac;


	if (argc == 5) {
		ac = 2;
		qp = ckopen(argv[1], "r");
	} else if (argc == 4) {	/* standar input; pipelined */
		ac = 1;	
		qp = stdin;
	} else fatalf("Usage: %s [frags_out] ms(<0) g(>=0) h(>0)",argv[0]);
	fscanf(qp,"%s %d\n",sname1,&M);
	fscanf(qp,"%s %d\n",sname2,&N);
	get_seq(sname1, &seq1);
	get_seq(sname2, &seq2);
	ms = DIGIT * atof(argv[ac]);
	if (ms >= 0)
		fatal("The mismatch weight should be negative.");
	G = DIGIT * atof(argv[ac+1]);
	if (G < 0)
		fatal("The gap-open penalty should be nonnegative.");
	H = DIGIT * atof(argv[ac+2]);
	if (H < 0)
		fatal("The gap-extend penalty should be positive.");

	/* set up match and mismatch weights */
	for (i = 0; i < 128; i++)
	   for (j = 0; j < 128; j++)
		if (i == j)
		   V[i][j] = DIGIT;
		else
		   V[i][j] = ms;

	fscanf(qp,"%d\n",&K);
	fscanf(qp,"%lf\n",&score);
	fscanf(qp,"%d %d %d\n",&fi,&fj,&fk);
	fscanf(qp,"%d %d %d\n",&li,&lj,&lk);
#ifdef STATS
	printf("ms=%.1f, g=%.1f, h=%.1f\n",(float)ms/DIGIT,(float)G/DIGIT,(float)H/DIGIT);
	printf("%s %d\n",sname1,M);
	printf("%s %d\n",sname2,N);
	printf("%.1f\n",score);
#endif
	printf("#:lav\n\n");
	printf("d {\n\"sim2 output with parameters:\n");
	printf("K = %d, ", K);
	printf("M = 1, I = V = %.1f,\nG = %.1f, H = %.1f\"\n}\n",
		((float)ms)/DIGIT, ((float)G)/DIGIT, ((float)H)/DIGIT);
	printf("s {\n  \"%s\" 1 %d\n  \"%s\" 1 %d\n}\n", sname1, M,
		sname2, N);

	Kcount = 0;
	while (fscanf(qp,"%d %d %d\n",&ci,&cj,&ck) != EOF) {
		if (ci == li) {
			if (fi == li) { /* only one fragment */
				printf("a {\n");
				printf("  s %.1f\n",(float)fk);
				printf("  b %d %d\n",fi+1,fj+1);
				printf("  e %d %d\n",li+lk,lj+lk);
				printf("  l %d %d %d %d %.1f\n",fi+1,fj+1,fi+fk,fj+fk,100.0);
				printf("}\n");
			} else {
				for (i = ppi; i <= ci; ++i) {
					LB[i] = min(ppj,LB[i]);
					RB[i] = max(cj,RB[i]);
				}
				len1 = ei - bi + fk + lk;
				len2 = ej - bj + fk + lk;
				S = (int *)ckalloc((len1+len2)*sizeof(int));
				for (i = 0; i < len1+len2; ++i) S[i]=0;
				c = ((float)ALIGN(seq1-1,seq2-1,bi,bj,ei,ej,LB,RB,V,G,H,S+fk))/DIGIT + (float)fk + (float)lk;
#ifdef STATS
				DISPLAY(seq1-1+fi,seq2-1+fj,len1,len2,S,fi+1,fj+1);
#endif
				printf("a {\n");
				printf("  s %.1f\n",c);
				printf("  b %d %d\n",fi+1,fj+1);
				printf("  e %d %d\n",li+lk,lj+lk);
				old_i = i = fi;
				old_j = j = fj;
				sp = S;
				match_no = 0;
				while (i < li+lk || j < lj+lk) {
					op = *sp++;
					if (op == 0) {
					   if (seq1[i++] == seq2[j++]) ++match_no;
					   if (i == li+lk || j == lj+lk || *sp != 0) {
					      printf("  l %d %d %d %d %.1f\n",old_i+1,old_j+1,i,j,100.0*((float)match_no)/((float)(i-old_i)));
					      match_no = 0;
					      old_i = i;
					      old_j = j;
					   }
					} else {
					   if (op > 0) {
						j = j+op;
						old_j = j;
					   } else {
						i = i-op;
						old_i = i;
					   }
					}
				}
				printf("}\n");
				free(S);
			}
			if (++Kcount == K) {
				if (fi != li) {
					free(LB+bi);
					free(RB+bi);
				}
				break;
			} else {
				if (fi != li) {
					free(LB+bi);
					free(RB+bi);
				}
				fscanf(qp,"%lf\n",&score);
				fscanf(qp,"%d %d %d\n",&fi,&fj,&fk);
				fscanf(qp,"%d %d %d\n",&li,&lj,&lk);
			}
		} else if (ci == fi) { /* fi != li */
			bi = fi + fk;
			bj = fj + fk;
			ei = li;
			ej = lj;
			t = (ei - bi + 1) * sizeof(int);
			LB = (int *) ckalloc(t) - bi;
			RB = (int *) ckalloc(t) - bi;
			for (i = bi; i <= ei; ++i) {
				LB[i] = ej;
				RB[i] = bj;
			}
			ppi = pi = bi;
			ppj = pj = bj;
		} else { /* intermediate fragments */
			for (i = ppi; i <= ci; ++i) {
				LB[i] = min(ppj,LB[i]);
				RB[i] = max(cj,RB[i]);
			}
			ppi = pi;
			ppj = pj;
			pi = ci+ck;
			pj = cj+ck;
		}
	}
	if (argc == 5) fclose(qp);
	exit(0);
}

/* get_seq - read a sequence */
int get_seq(char *file_name, char **seqptr)
{
	FILE *qp, *ckopen();
	char *p, *fgets(), *strchr(), *ckalloc();
	int c, n;
		
	qp = ckopen(file_name, "r");
	for (n=0; (c=getc(qp)) != EOF;)
		if (c != '\n') ++n;
	*seqptr = ckalloc(n+1);
	rewind(qp);
	p = *seqptr;
	for (p = *seqptr; ; ) {
		if (fgets(p, (int)n, qp) == NULL) {
			fclose(qp);
			*p = '\0';
			break;
		}
		if (p == *seqptr && !isupper(*p))
			continue;
		while (isupper(*p)) {
			if (strchr("ACGT", *p) == NULL)
				*p = 'A';
			++p;
		}
		if (*p != '\n' && *p != '\0') {
			fprintf(stderr, "Illegal character %c in file %s.\n",
				*p, file_name);
			exit(1);
		}
	}
	*p = '\0';
	return p - *seqptr;
}

/* lib.c - library of C procedures. */

/* fatal - print message and die */
fatal(char *msg)
{
	fprintf(stderr, "%s\n", msg);
	exit(1);
}

/* fatalf - format message, print it, and die */
fatalf(char *msg, char *val)
{
	fprintf(stderr, msg, val);
	putc('\n', stderr);
	exit(1);
}
	
/* ckopen - open file; check for success */
FILE *ckopen(char *name, char *mode)
{
	FILE *fopen(), *fp;

	if ((fp = fopen(name, mode)) == NULL)
		fatalf("Cannot open %s.", name);
	return(fp);
}

/* ckalloc - allocate space; check for success */
char *ckalloc(long amount)
{
	char *malloc(), *p;

	if ((p = malloc( (unsigned) amount)) == NULL)
		fatal("Ran out of memory.");
	return(p);
}

/* strsame - tell whether two strings are identical */
int strsame(char *s, char *t)
{
	return(strcmp(s, t) == 0);
}