/* dottyplot.c -- c version of dotplot routine 

 by d.g.gilbert May 1989
 
 28sep95 -- removed Macintosh dependencies  
 -- revise to do line segments along diagonals that have contiguous points
--------------------------------------*/

#define TESTPROG

#include <stdlib.h>

#define min(a,b)  (a<b)?a:b
#define max(a,b) (a>b)?a:b


/* 
makedotplot creates a dot-plot array from biosequence data
by comparing two sequences for matching bases.

Result is array of long integers.  
Each long integer is actually an x,y short value, as in

  typedef  struct {
   short  x;
   short  y;
   } xydot, *xydots;
  
Parameters include 
 - vseq, nv, vertical sequence and its length
 - hseq, nh, horizontal sequence and length
 - window size (commonly in the 10-20 base range),
 - nmatch, # matches required for a window
 - doself, vertical & horizontal sequence is same
 - allpts, plot all matching bases in a window (versus 1 dot / window)
 
*/
 
  
long* dottyplot( 
  char *vseq,  long nv,  /* vertical sequence */
  char *hseq,  long nh,  /* horizontal sequence */
  short window,  /* size of sliding window, #bases */
  short nmatch, /* number of that must match in window for a dot */
  char doself,   /* vert. & horizontal are same sequence */
  char allpts,   /* display all points in a window (or just one/window) */
  long *ndots /* return number of x,y dots in result */
  ) 
 {
  register char *matchp, *hp, *vp;
  register long i;
  char  *matchvec;
  long   matches, diag, xstart, ystart, yfinish, halfwind, 
       h, v, ndiag, j, lasti, n2;
  long  npic, maxpic;
  long  *xydots;
 
    /* initialize dot window */
  matchvec = malloc( nv+window+10);
  n2= nh+nv-1;
   halfwind = window >> 1;
  
    /* initialize xy dots */
  maxpic = 100;
  npic  = 0;
  xydots = malloc( maxpic * sizeof(long));
  
       /* Compare for each register shift (diagonal) */
  if (doself) ndiag= nv; else ndiag = n2;
  for (diag=1; diag<=ndiag; diag++) {
    ystart=  max(0, nv-diag);
   yfinish= min( nv-1, n2-diag);
   xstart=  ystart - nv + diag;
   
     /* Compare the two sequences along this diagonal
      and fill in matchVec while clearing pointVec */
   matchp= matchvec;
   hp= &(hseq[xstart]); 
   vp= &(vseq[ystart]);  
   for (i= ystart; i<= yfinish; i++) {
     *(matchp++)= (*(hp++) == *(vp++));
    } 
   
   /*!!! ^^^^^ revise this & data store to find/keep line segments !!! */
   /* if matchp[i] && matchp[i-1] ... */
   
         /* clear a window's width at the end */
    for (i= 0; i<window; i++) *(matchp++)= 0;  
   
      /* Count the number of matches in the first window */
   for (i=0, matchp= matchvec, matches=0; i<window; i++) {
    if (*(matchp++)) matches++;
    }
    
     /* step window along diagonal, transfer comparison 
      results to pointVec when nMatch is met. */
   h = xstart;  
   v = nv - ystart;
   if (!allpts) { h += halfwind; v -= halfwind; }
   for (i=lasti=0; i< yfinish-ystart-window; i++) {   
    if (matches >= nmatch) { 
    
        if (allpts) {
      for (j= max(i,lasti); j<i+window; j++) {
       if (matchvec[j]) {
        if (npic >= maxpic) {  
         maxpic += 100;  
         xydots= (long*) realloc( xydots, maxpic*sizeof(long)); 
         if (!xydots) goto done;
         } 
        xydots[npic++] = (h+j) + ((v-j)<<16); 
        }
       }
      lasti= i+window;
      }
      
     else { 
      /* MoveTo( h + i, v - i); Line(0, 0); */  
      if (npic >= maxpic) {  
       maxpic += 100;  
       xydots= (long*) realloc( xydots, maxpic*sizeof(long)); 
       if (!xydots) goto done;
       } 
       xydots[npic++] = (h+i) + ((v-i)<<16);
      }
     }
    matches += (matchvec[i + window] - matchvec[i]);
    }
   } /* diag */
   
done:   
    free( matchvec); 
     /* resize to data size */
  xydots= (long*) realloc( xydots, npic*sizeof(long)); 
  if (ndots) *ndots= npic;
  return xydots;
} /* dottyplot */
   
   
#ifdef TESTPROG

#include <stdio.h>

typedef  struct {
 short  x;
 short  y;
 } xydot, *xydots;


void main( int argc, char *argv[])
{
 char *vseq = NULL;  long nv = 0;  /* vertical sequence */
 char *hseq = NULL;  long nh = 0;  /* horizontal sequence */
 short window = 10;  /* size of sliding window, #bases */
 short nmatch = 5; /* number of that must match in window for a dot */
 int doself = 0;   /* vert. & horizontal are same sequence */
 int allpts = 0;   /* display all points in a window (or just one/window) */
 char * infile = NULL;
 int   i;
  FILE  * fp;
  xydot * dots; 
  long ndots; 
 char *atseq, *seq = NULL; 
 long  maxseq, nseq;
  char inseq;
  
  
 printf("Dottyplot test\n");
 for (i=0; i<argc; i++) {
    if (*argv[i] == '-')
     switch (argv[i][1]) {
     case 'a': allpts= 1; break;
     case 's': doself= 1; break;
     case 'w': i++; window= atoi(argv[i]); break;
     case 'm': i++; nmatch= atoi(argv[i]); break;
     }
    else {
     infile= argv[i];
     }
    }
    
 if (!infile) {
  printf("Usage:  dottyplot [-a] [-s] [-w 10] [-m 5] sequencefile \n");
  printf(" sequencefile in fasta format, with 1 or 2 sequences\n");
  printf("  -a[llpoints]  = plot all points\n");
  printf("  -s[elf]       = match sequence against self\n");
  printf("  -w[indow] #   = set window size\n");
  printf("  -m[atch]  #   = set minimum matches per window\n");
  }
  
 inseq= 0;
 seq= NULL;
 nseq= 0;
 fp= fopen(infile, "r");
 while (! feof(fp)) {
  char  ins[1024];
  char *cp;
  
  cp= fgets( ins, sizeof(ins), fp);
  if (cp) {
  
   if (*cp == '>') {
    /* new seq */
    inseq++;
    if (inseq == 2) { hseq= seq; nh= nseq; }
    else if (inseq > 2) goto dodots;
    maxseq= 500;
    nseq= 0;
    seq= malloc( maxseq);
    atseq= seq;
    }
    
   else if (inseq) while (*cp >= ' ') {
    if (*cp > ' ') {
     if (nseq > maxseq) {
      maxseq += 100;
      seq= realloc(seq, maxseq);
      atseq= seq + nseq;
      }
     *atseq++ = *cp;
     nseq++;
     }
    cp++;
    }
   else {
    
    }
  }
 }
 
dodots:
 fclose(fp);
 if (!nseq) { printf("No sequence.\n"); exit(-1); }
 vseq= seq; nv= nseq;
 
 if (doself)  
  dots= (xydot*)dottyplot( vseq, nv, vseq, nv, window, nmatch, doself, allpts, &ndots);  
 else  
  dots= (xydot*)dottyplot( vseq, nv, hseq, nh, window, nmatch, doself, allpts, &ndots);  
  
 
  if (dots) {
   printf("Dot plot (x,y) points\n");
   for (i= 0; i<ndots; i++) {
    printf("%d\t%d\n", dots[i].x, dots[i].y);
    }
   }
}

#endif /*TESTPROG*/