/* sample program for invoking g_any.c */
#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(argc, argv)
int argc;
char *argv[];
{
	double atof();
	char *ckalloc();
	char *seq1, *seq2;
	int M, N;		/* lengths of sequences */
	long ms, Q, R, V[128][128];
	int *LB, *RB;
	long *S;		/* conversion operations */
	int i, j, b_argc;

	if (argc < 7)
		fatalf("%s seq1 seq2 boundary ms(<0) g(>=0) h(>0)",argv[0]);
	M = get_seq(argv[1], &seq1);
	N = get_seq(argv[2], &seq2);
	j = (M+1)*sizeof(int);
	LB = (int *) ckalloc(j);
	RB = (int *) ckalloc(j);
	b_argc = 4;
	if (strsame(argv[3],"rec")) set_rec_boundary(M, N, LB, RB);
	else if (strsame(argv[3],"band")) {
		if (argc != 9) 
		    fatalf("%s seq1 seq2 band low up ms(<0) g(>=0) h(>0)",argv[0]);
		set_band_boundary(M, N, atoi(argv[4]), atoi(argv[5]), LB, RB);
		b_argc = 6;
	} else if (strsame(argv[3],"ran")) {
		if (argc != 8)
		    fatalf("%s seq1 seq2 ran # ms(<0) g(>=0) h(>0)",argv[0]);
		set_ran_boundary(M, N, atoi(argv[4]), LB, RB);
		b_argc = 5;
	} else get_boundary(argv[3], M, LB, RB);
	/*
	for (i=0; i<=M; ++i) printf("%5d, %5d\n",LB[i],RB[i]);
	*/
	ms = DIGIT * atof(argv[b_argc]);
	if (ms >= 0)
		fatal("The mismatch weight should be negative.");
	Q = DIGIT * atof(argv[b_argc+1]);
	if (Q < 0)
		fatal("The gap-open penalty should be nonnegative.");
	R = DIGIT * atof(argv[b_argc+2]);
	if (R < 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;

	S = (long *)ckalloc((M+N)*sizeof(long));

	printf("seq1: %s   length=%d\n",argv[1],M);
	printf("seq2: %s   length=%d\n",argv[2],N);
	printf("q=%.1f, r=%.1f\n",(float)Q/DIGIT,(float)R/DIGIT);
	printf("\nOptimal global score = %.1f\n",
		((float)ALIGN(seq1-1,seq2-1,M,N,LB,RB,V,Q,R,S))/DIGIT);
	DISPLAY(seq1-1,seq2-1,M,N,S,1,1);
	exit(0);
}

/* get_seq - read a sequence */
int get_seq(file_name, seqptr)
char *file_name, **seqptr;
{
	FILE *qp, *ckopen();
	char *p, *fgets(), *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;
	while((c=getc(qp)) != EOF) 
		if (c != '\n') *p++ = c;
	*p = '\0';
	fclose(qp);
	return p - *seqptr;
}
	
/* get_boundary - read the boundary of the region */
get_boundary(file_name, rc, lb, rb)
char *file_name;
int rc, *lb, *rb;
{
	FILE *fp, *ckopen();
	int i;

	fp = ckopen(file_name, "r");
	for (i=0; i<=rc; ++i, ++lb, ++rb)
		fscanf(fp,"%i %i",lb,rb);
	fclose(fp);
	return;
}

/* set_rec_boundary - setup the boundary for the rectangular case */
set_rec_boundary(rc, cc, lb, rb)
int rc, cc, *lb, *rb;
{
	int i;

	for (i=0; i<=rc; ++i) {
		*lb++ = 0;
		*rb++ = cc;
	}
	return;
}

/* set_band_boundary - setup the boundary for the band */
set_band_boundary(rc, cc, low, up, lb, rb)
int rc, cc, low, up, *lb, *rb;
{
	int i;

	low = min(max(-rc, low),min(0,cc-rc));
	up = max(min(cc, up), max(cc-rc, 0));
	for (i=0; i<=rc; ++i) {
		*lb++ = min(cc,max(0,i+low));
		*rb++ = max(0,min(cc,i+up));
	}
	return;
}

/* set_ran_boundary - setup the boundary for testing */
set_ran_boundary(rc, cc, n, lb, rb)
int rc, cc, n, *lb, *rb;
{
	int i, clb, crb, cn;

	if (n == 0) fatal("intial value should be greater than 0\n");
	if (n >= cc) fatal ("intial value is too large!\n");
	cn = cc/(rc/n+1);
	clb = 0;
	crb = cc-rc/n;
	if (crb < 0) crb = 0;
	for (i=0; i<rc; ++i) {
		*lb++ = clb;
		*rb++ = crb;
		if (i % n == n-1) {
			clb += cn;
			clb = min(crb,min(cc,clb));
			crb += 1;
			crb = min(cc,crb);
		}
	}
	*lb = clb;
	*rb = cc;
	return;
}

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

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

/* fatalf - format message, print it, and die */
fatalf(msg, val)
char *msg, *val;
{
	fprintf(stderr, msg, val);
	putc('\n', stderr);
	exit(1);
}
	
/* ckopen - open file; check for success */
FILE *ckopen(name, mode)
char *name, *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(amount)
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(s, t)
char *s, *t;
{
	return(strcmp(s, t) == 0);
}