#if 0
	cc -o gen2 gen2.c -lm
	exit
#endif

/*
 * Genome comparison program.
 * Ken Shirriff				shirriff@sprite.Berkeley.EDU
 * Usage: gen file1 file2 > output.pgm
 * Output is a PGM image that can be displayed with xloadimage or xv.
 * You may have to play with the gamma to make the output visible
 * I.e. xloadimage -gamma 2 output.pgm
 * Flags are:
 *	-all	Normally only the parts of the sequence near each other
 *		are compared; -all compares the whole sequences (slow).
 *	-num n	Specifies how many bases are in each compared substring.
 *		If n is too small, there's a lot of noise.  If n is too big,
 *		there won't be any matches at all.
 *	-scale s	Controls how big the resulting image is.
 *	-brighten f	Will brighten the resulting image by a factor of f.
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifdef sun
#define bcmp(a,b,c) memcmp(a,b,c)
#define bzero(a,b) memset(a,0,b)
#endif

int all=0;
int n=5;
int raw=1;
int grid=0;
#define GSPC 100
float brighten = 1;
float gam = 3.;

#define MAX 10000
int sizex=500, sizey=500;
int scale=20;
char *rawbuf;

main(argc,argv)
int argc;
char **argv;
{
    char seq0[MAX], seq1[MAX];
    int nseq0,nseq1;
    while (argc>1 &&argv[1][0]=='-') {
#define USAGE "gen [-xgraph] [-all] [-grid] [-num n] [-scale n] [-brighten n] [-gamma g] file1 file2\n"
	if (strcmp(argv[1],"-all")==0) {
	    all = 1;
	    argc--;
	    argv++;
	} else if (strcmp(argv[1],"-grid")==0) {
	    grid = 1;
	    argc--;
	    argv++;
	} else if (strcmp(argv[1],"-xgraph")==0) {
	    raw = 0;
	    argc--;
	    argv++;
	} else if (argc>2 && strcmp(argv[1],"-scale")==0) {
	    scale = atoi(argv[2]);
	    argc -= 2;
	    argv += 2;
	} else if (argc>2 && strcmp(argv[1],"-gamma")==0) {
	    gam = atof(argv[2]);
	    argc -= 2;
	    argv += 2;
	} else if (argc>2 && strcmp(argv[1],"-brighten")==0) {
	    brighten = atof(argv[2]);
	    argc -= 2;
	    argv += 2;
	} else if (argc>2 && strcmp(argv[1],"-num")==0) {
	    n = atoi(argv[2]);
	    argc -= 2;
	    argv += 2;
	} else {
	    fprintf(stderr,USAGE);
	    exit(-1);
	}

    }
    if (argc<3) {
	fprintf(stderr,USAGE);
	exit(-1);
    }

    readseq(seq0,&nseq0,argv[1]);
    readseq(seq1,&nseq1,argv[2]);
    sizex = nseq0/scale+1;
    sizey = nseq1/scale+1;

    if (raw==0) {
	printf("PixelMarkers: 1\n");
	printf("NoLines: 1\n\n");
    } else {
	rawbuf = (char *)malloc(sizex*sizey);
	if (rawbuf==0) {
		fprintf(stderr,"Malloc failed: size %d x %d\n", sizex,
		sizey);
	}
	bzero(rawbuf,sizex*sizey);
    }

    cmpseq(seq0,nseq0,seq1,nseq1);
    if (raw) {
	dumpbuf();
    }
}

readseq(seq,nseq,name)
char *seq;
int *nseq;
char *name;
{
    FILE *f;
    char buf[80];
    f = fopen(name,"r");
    if (f==NULL) {
	fprintf(stderr,"Couldn't open %s\n", name);
	exit(-1);
    }
    *nseq = 0;
    while (1) {
	if (fgets(buf,79,f)==NULL) break;
	if (strncmp(buf,"ORIGIN",6)==0) {
	    break;
	}
    }
    while (1) {
	if (fgets(buf,79,f)==NULL) break;
	if (buf[0]=='/') break;
	doline(seq,nseq,buf);
    }
    fclose(f);
}

doline(seq,nseq,buf)
char *seq,*buf;
int *nseq;
{
    int i;
    for (i=0;i<80;i++) {
	if (isalpha(buf[i])) {
	    (*nseq)++;
	    seq[*nseq] = buf[i];
	} else if (buf[i]=='\0') {
	    break;
	}
    }
}

cmpseq(seq0,nseq0,seq1,nseq1)
char *seq0,*seq1;
int nseq0,nseq1;
{
    int i,j;
    int j0,j1;
    for (i=0;i<nseq0-n;i++) {
    if (all) {
        j0 = 0;
	j1 = nseq1-n;
    } else {
	j0 = i-2000;
	j1 = i+2000;
    }
	if (j0<0) j0=0;
	if (j1>nseq1-n) j1=nseq1-n;
	for (j=j0;j<j1;j++) {
	    if (bcmp(seq0+i,seq1+j,n)==0) {
		if (raw==0) {
		    printf("%d %d\n",i,j);
		} else {
		    rawbuf[(i/scale)+sizex*(j/scale)]++;
		}
	    }
	}
	if ((i%50)==0) {
	    fprintf(stderr,"Done %d of %d\n",i,nseq0);
	}
    }

    if (grid && raw) {
	for (i=0;i<nseq0;i+=GSPC) {
	    for (j=0;j<nseq1;j++) {
		    rawbuf[(i/scale)+sizex*(j/scale)] = 1;
	    }
	}
	for (j=0;j<nseq1;j+=GSPC) {
	    for (i=0;i<nseq0;i++) {
		    rawbuf[(i/scale)+sizex*(j/scale)] = 1;
	    }
	}
    }
}

dumpbuf()
{
    int i,v;
    int max=1;
    printf("P5\n%d %d\n%d\n",sizex,sizey,255);
    for (i=0;i<sizex*sizey;i++) {
	if (gam != 1) {
		v = pow(rawbuf[i]*1./scale,1./gam)*255.*brighten;
	} else {
		v = rawbuf[i]*255./scale*brighten;
	}
	if (v>255) {
	    putchar(255);
	} else {
	    putchar(v);
	}
    }
}