/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */
/* From input file "denri.p" */


#include <p2c/p2c.h>


/* denri: compute density for delta Ri

  Tom Schneider
  NCI/FCRDC Bldg 469. Room 144
  P.O. Box B
  Frederick, MD  21702-1201
  (301) 846-5581 (-5532 for messages)
  toms@ncifcrf.gov
  http://www.lmmb.ncifcrf.gov/~toms/

  National Cancer Institute
  Laboratory of Mathematical Biology

 */

/* end of program */

/* begin module version */

#define version         1.09
/* of denri.p 1999 June 28
1.09: 1999 June 28: documentation brush up.
1.08: 1997 April 1: major use
origin 1997 March 31 */
#define updateversion   1.00
    /* defines lowest acceptable current parameter file */
/* end module version */

/* begin module describe.denri */
/*
name
   denri: compute density for delta Ri

synopsis
   denri(data: in, denrip: in, den: out, xyplop: out, output: out)

files
   data:  output of scan program

   den:  output of this program.  Use it as the xyin input to xyplo.

   denrip:  parameters to control the program.  The file must contain the
        following parameters, one per line:

      0. The version number of the program.  This allows the user to be
         warned if an old parameter file is used.

      1. numbertoprocess: number of data items to process (negative = all)

      2. wantfrom, wantto: position from-to range for computing

      3. wantminRi, wantmaxRi: bit from-to range for computing

   xyplop: control file for xyplo

   output: messages to the user

description

   Take the delta Ri values from the scan program and compute the density at
   each position and bit level.  Bit values from the scan data file are
   rounded and counts are placed in bins.  The output can then be plotted
   with xyplo.

examples

   denrip:

1.01      version of denri that this parameter file is designed for.
1000      number of data items to process (negative = all)
-20 20    wantfrom, wantto: position from-to range for computing
-70 20    wantminRi, wantmaxRi: bit from-to range for computing

documentation

see also
   scan.p, xyplo.p

author
   Thomas Dana Schneider

bugs

technical notes
   Four constants in module denri.const define the largest
   position range and bit range the program can handle.

*/
/* end module describe.denri */

/* begin module denri.const */

#define rangefrom       (-500)   /* lowest position allowed */
#define rangeto         500   /* highest position allowed */
#define minRi           (-70)   /* lowest Ri allowed */
#define maxRi           20   /* highest Ri allowed */


/* end module denri.const */

Static _TEXT data;   /* file used by this program */
Static _TEXT denrip;   /* file used by this program */
Static _TEXT den;   /* file used by this program */
Static _TEXT xyplop;   /* file used by this program */


Static jmp_buf _JL1;


/* begin module halt */
Static Void halt()
{
  /* stop the program.  the procedure performs a goto to the end of the
     program.  you must have a label:
        label 1;
     declared, and also the end of the program must have this label:
        1: end.
     examples are in the module libraries.
     this is the only goto in the delila system. */
  printf(" program halt.\n");
  longjmp(_JL1, 1);
}


/* end module halt version = 4.20; (@ of prgmod.p 1997 March 15 */

/* begin module skipblanks */
Static Void skipblanks(thefile)
_TEXT *thefile;
{
  /* skip over blanks until a non-blank, or end of line, is found */
  while ((P_peek(thefile->f) == ' ') & (!P_eoln(thefile->f)))
    getc(thefile->f);
}


Static Void skipnonblanks(thefile)
_TEXT *thefile;
{
  /* skip over nonblanks until a blank, or end of line, is found */
  while ((P_peek(thefile->f) != ' ') & (!P_eoln(thefile->f)))
    getc(thefile->f);
}


Static Void skipcolumn(thefile)
_TEXT *thefile;
{
  /* skip over a data column */
  skipblanks(thefile);
  skipnonblanks(thefile);
}


#define wid             7   /* character width of data output */
#define dec             5   /* decimals for data output */


/* Local variables for themain: */
struct LOC_themain {
  _TEXT *data, *den;
  long bins[rangeto - rangefrom + 1][maxRi - minRi + 1];
  /* bins to store the data */
  long c;   /* count number of data items to fill in */
  long d;   /* index to delta Ri values */
  double deltaRi;   /* delta Ri values */
  long numbertoprocess;   /* number of data items to process */
  double parameterversion;   /* parameter version number */
  long position;   /* position in the sequences */
  long wantfrom, wantto;   /* range to process */
  long wantminRi, wantmaxRi;   /* bit range to process */
} ;

Local Void readparameters(denrip, LINK)
_TEXT *denrip;
struct LOC_themain *LINK;
{
  /* read user defined parameters */
  if (*denrip->name != '\0') {
    if (denrip->f != NULL)
      denrip->f = freopen(denrip->name, "r", denrip->f);
    else
      denrip->f = fopen(denrip->name, "r");
  } else
    rewind(denrip->f);
  if (denrip->f == NULL)
    _EscIO2(FileNotFound, denrip->name);
  RESETBUF(denrip->f, Char);
  fscanf(denrip->f, "%lg%*[^\n]", &LINK->parameterversion);
  getc(denrip->f);
  if (LINK->parameterversion < updateversion) {
    printf("You have an old parameter file!\n");
    halt();
  }

  fscanf(denrip->f, "%ld%*[^\n]", &LINK->numbertoprocess);
  getc(denrip->f);
  fscanf(denrip->f, "%ld%ld%*[^\n]", &LINK->wantfrom, &LINK->wantto);
  getc(denrip->f);
  fscanf(denrip->f, "%ld%ld%*[^\n]", &LINK->wantminRi, &LINK->wantmaxRi);
  getc(denrip->f);

  if (LINK->wantfrom < rangefrom) {
    printf("wantfrom cannot be less than %ld\n", (long)rangefrom);
    halt();
  }

  if (LINK->wantto > rangeto) {
    printf("wantto cannot be more than %ld\n", (long)rangeto);
    halt();
  }

  if (LINK->wantfrom > LINK->wantto) {
    printf("wantfrom cannot be more than wantto\n");
    halt();
  }

  if (LINK->wantminRi < minRi) {
    printf("wantminRi cannot be less than %ld\n", (long)minRi);
    halt();
  }

  if (LINK->wantmaxRi > maxRi) {
    printf("maxRi cannot be more than %ld\n", (long)maxRi);
    halt();
  }

  if (LINK->wantminRi > LINK->wantmaxRi) {
    printf("wantminRi cannot be more than wantmaxRi\n");
    halt();
  }

}

Local Void writeparameters(out, LINK)
_TEXT *out;
struct LOC_themain *LINK;
{
  fprintf(out->f, "* denri %4.2f\n", version);
  fprintf(out->f, "* %9.2f parameterversion for denri\n",
	  LINK->parameterversion);
  fprintf(out->f, "* %9ld numbertoprocess\n", LINK->numbertoprocess);
  fprintf(out->f, "* %4ld %4ld wantfrom, wantto\n",
	  LINK->wantfrom, LINK->wantto);
  fprintf(out->f, "* %4ld %4ld wantminRi, wantmaxRi\n",
	  LINK->wantminRi, LINK->wantmaxRi);
}

Local Void clearbins(LINK)
struct LOC_themain *LINK;
{
  /* clear the bins out */
  long FORLIM, FORLIM1;

  printf("clearing bins ...\n");
  FORLIM = LINK->wantto;
  for (LINK->position = LINK->wantfrom; LINK->position <= FORLIM; LINK->
	 position++) {
    FORLIM1 = LINK->wantmaxRi;
    for (LINK->d = LINK->wantminRi; LINK->d <= FORLIM1; LINK->d++)
      LINK->bins[LINK->position - rangefrom][LINK->d - minRi] = 0;
  }
  printf("cleared\n");
}

Local Void fillonedatum(LINK)
struct LOC_themain *LINK;
{
  /* read and store one data item */
  if (P_peek(LINK->data->f) == '*') {
    fscanf(LINK->data->f, "%*[^\n]");
    getc(LINK->data->f);
    return;
  }
  skipcolumn(LINK->data);
  skipcolumn(LINK->data);
  skipcolumn(LINK->data);
  fscanf(LINK->data->f, "%ld", &LINK->position);
  skipcolumn(LINK->data);
  skipcolumn(LINK->data);
  skipcolumn(LINK->data);
  skipcolumn(LINK->data);
  fscanf(LINK->data->f, "%lg", &LINK->deltaRi);
  fscanf(LINK->data->f, "%*[^\n]");
  getc(LINK->data->f);
  LINK->d = (long)floor(LINK->deltaRi + 0.5);

  /*
        if (d < minRi) or (d > maxRi)
           or
           (position < rangefrom) or (position > rangeto)
        then begin
           writeln(output,'position = ',position:1,' d=',d:1);
           writeln(output,'is out of the storage area');
halt
        end
        else begin
writeln(output);
writeln(output,'position = ',position:1,' d=',d:1);
writeln(output,' bins[position, d] = ', bins[position, d]:1);
writeln(output,'succ is: ',succ(bins[position, d]):1);
writeln(output,'filled ',position:1,' ',d:1);
  */

  if (LINK->d >= LINK->wantminRi && LINK->d <= LINK->wantmaxRi &&
      LINK->position >= LINK->wantfrom && LINK->position <= LINK->wantto)
    LINK->bins[LINK->position - rangefrom][LINK->d - minRi]++;
}

Local Void fillbins(LINK)
struct LOC_themain *LINK;
{
  /* fill the bins with data */
  long FORLIM;

  printf("filling bins ...\n");
  if (*LINK->data->name != '\0') {
    if (LINK->data->f != NULL)
      LINK->data->f = freopen(LINK->data->name, "r", LINK->data->f);
    else
      LINK->data->f = fopen(LINK->data->name, "r");
  } else
    rewind(LINK->data->f);
  if (LINK->data->f == NULL)
    _EscIO2(FileNotFound, LINK->data->name);
  RESETBUF(LINK->data->f, Char);
  if (LINK->numbertoprocess < 0) {
    while (!BUFEOF(LINK->data->f))
      fillonedatum(LINK);
  } else {
    FORLIM = LINK->numbertoprocess;
    for (LINK->c = 1; LINK->c <= FORLIM; LINK->c++)
      fillonedatum(LINK);
  }
  printf("filled\n");
}

Local Void process(LINK)
struct LOC_themain *LINK;
{
  /* process and output results */
  long b;   /* bin value */
  long sum;   /* for determining the sum of data at a position */
  double r;   /* bin ratio = b/sum */
  double hue;   /* conversion of r to a color */
  double saturation;   /* conversion of r to a color */
  double brightness;   /* conversion of r to a color */
  long FORLIM, FORLIM1;

  /*zzz*/
  printf("processing...\n");
  fprintf(LINK->den->f, "* columns:\n");
  fprintf(LINK->den->f, "* 1: position\n");
  fprintf(LINK->den->f, "* 2: bits\n");
  fprintf(LINK->den->f, "* 3: total counts at this position\n");
  fprintf(LINK->den->f,
    "* 4: (counts at this bit-position) / (total counts at this position)\n");
  fprintf(LINK->den->f, "* 5: 1\n");
  fprintf(LINK->den->f, "* 6: 1\n");
  fprintf(LINK->den->f, "* 7: hue\n");
  fprintf(LINK->den->f, "* 8: saturation\n");
  fprintf(LINK->den->f, "* 9: brightness\n");
  FORLIM = LINK->wantto;
  for (LINK->position = LINK->wantfrom; LINK->position <= FORLIM; LINK->
	 position++) {
    sum = 0;
    FORLIM1 = LINK->wantmaxRi;
    for (LINK->d = LINK->wantminRi; LINK->d <= FORLIM1; LINK->d++)
      sum += LINK->bins[LINK->position - rangefrom][LINK->d - minRi];
    if (sum > 0) {
      FORLIM1 = LINK->wantmaxRi;
      for (LINK->d = LINK->wantminRi; LINK->d <= FORLIM1; LINK->d++) {
	b = LINK->bins[LINK->position - rangefrom][LINK->d - minRi];
	/*zzz*/
	if (b == 0) {
	  brightness = 0.9;
	  saturation = 0.0;
	} else {
	  r = (double)b / sum;
	  hue = 0.84 * r + 0.16;
	  saturation = 1.0;
	  brightness = 1.0;
	}
	/* 1 */
	/* 2 */
	/* 3 */
	/* 4 */
	/* 5 */
	/* 6 */
	/* 7 */
	/* 8 */
	fprintf(LINK->den->f,
		"%4ld %3ld %*ld %*ld %*d %*d %*.*f %*.*f %*.*f\n",
		LINK->position, LINK->d, wid, sum, wid, b, wid, 1, wid, 1,
		wid, dec, hue, wid, dec, saturation, wid, dec, brightness);
	    /* 9 */
      }
    }
  }
  printf("processed\n");
}

#undef wid
#undef dec

Local Void mkxyplop(xyplop, LINK)
_TEXT *xyplop;
struct LOC_themain *LINK;
{
  /* make the xyplop file */
  long ix, iy;   /* intervals for x and y */

  ix = LINK->wantto - LINK->wantfrom;
  while (ix > 20)
    ix = (long)(ix / 10.0);
  iy = (long)((LINK->wantmaxRi - LINK->wantminRi) / 10.0);
  if (*xyplop->name != '\0') {
    if (xyplop->f != NULL)
      xyplop->f = freopen(xyplop->name, "w", xyplop->f);
    else
      xyplop->f = fopen(xyplop->name, "w");
  } else {
    if (xyplop->f != NULL)
      rewind(xyplop->f);
    else
      xyplop->f = tmpfile();
  }
  if (xyplop->f == NULL)
    _EscIO2(FileNotFound, xyplop->name);
  SETUPBUF(xyplop->f, Char);
  fprintf(xyplop->f, "3 6       zerox zeroy         \n");
  fprintf(xyplop->f, "x %ld %ld zx min max          \n",
	  LINK->wantfrom, LINK->wantto);
  fprintf(xyplop->f, "y %ld %ld  zy min max          \n",
	  LINK->wantminRi, LINK->wantmaxRi);
  fprintf(xyplop->f, "%ld %ld 2 2 xinterval yinterval \n", ix, iy);
  fprintf(xyplop->f, "4 4       xwidth    ywidth    \n");
  fprintf(xyplop->f, "0 0       xdecimal  ydecimal  \n");
  fprintf(xyplop->f, "15 15     xsize     ysize     \n");

  fprintf(xyplop->f, "position, bases\n");
  fprintf(xyplop->f, "deltaRi, bits.");
  fprintf(xyplop->f, "  number of sequences at position 0: %ld\n",
	  LINK->bins[-rangefrom][-minRi]);

  fprintf(xyplop->f, "a         zc                  \n");
  fprintf(xyplop->f, "n 2       zxl base            \n");
  fprintf(xyplop->f, "n 2       zyl base            \n");
  fprintf(xyplop->f, "          ********************\n");
  fprintf(xyplop->f, "1 2       xcolumn   ycolumn   \n");
  fprintf(xyplop->f, "0         symbol column       \n");
  fprintf(xyplop->f, "5  6      xscolumn  yscolumn  \n");
  fprintf(xyplop->f, "7 8 9     hue saturation brightnes\n");
  fprintf(xyplop->f, "          ********************\n");
  fprintf(xyplop->f, "R         symbol-to-plot      \n");
  fprintf(xyplop->f, "i         symbol-flag         \n");
  fprintf(xyplop->f, "-1.0      symbol sizex        \n");
  fprintf(xyplop->f, "-1.0      symbol sizey        \n");
  fprintf(xyplop->f, "n         connection size     \n");
  fprintf(xyplop->f, "n  0.125  linetype  size      \n");
  fprintf(xyplop->f, "          ********************\n");
  fprintf(xyplop->f, ".\n");
  fprintf(xyplop->f, "          **** version 8.50 of\n");
  fprintf(xyplop->f, "l 0 0 0.125  User defined line\n");
  /*
writeln(xyplop,'l -1 0 0.125  User defined line');
*/
}


/* end module skipblanks version = 4.20; (@ of prgmod.p 1997 March 15 */

/* begin module denri.themain */
Static Void themain(data_, denrip, den_, xyplop)
_TEXT *data_, *denrip, *den_, *xyplop;
{
  /* the main procedure of the program */
  struct LOC_themain V;
  _TEXT TEMP;

  V.data = data_;
  V.den = den_;
  printf("denri %4.2f\n", version);
  readparameters(denrip, &V);
  if (*V.den->name != '\0') {
    if (V.den->f != NULL)
      V.den->f = freopen(V.den->name, "w", V.den->f);
    else
      V.den->f = fopen(V.den->name, "w");
  } else {
    if (V.den->f != NULL)
      rewind(V.den->f);
    else
      V.den->f = tmpfile();
  }
  if (V.den->f == NULL)
    _EscIO2(FileNotFound, V.den->name);
  SETUPBUF(V.den->f, Char);
  writeparameters(V.den, &V);
  TEMP.f = stdout;
  *TEMP.name = '\0';
  writeparameters(&TEMP, &V);
  clearbins(&V);
  fillbins(&V);
  process(&V);
  mkxyplop(xyplop, &V);
  printf("done\n");
}


/* end module denri.themain */

main(argc, argv)
int argc;
Char *argv[];
{
  PASCAL_MAIN(argc, argv);
  if (setjmp(_JL1))
    goto _L1;
  xyplop.f = NULL;
  strcpy(xyplop.name, "xyplop");
  den.f = NULL;
  strcpy(den.name, "den");
  denrip.f = NULL;
  strcpy(denrip.name, "denrip");
  data.f = NULL;
  strcpy(data.name, "data");
  themain(&data, &denrip, &den, &xyplop);
_L1:
  if (data.f != NULL)
    fclose(data.f);
  if (denrip.f != NULL)
    fclose(denrip.f);
  if (den.f != NULL)
    fclose(den.f);
  if (xyplop.f != NULL)
    fclose(xyplop.f);
  exit(EXIT_SUCCESS);
}






















/* End. */