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


#include <p2c/p2c.h>


/* aran: aligned random sequences
   by Thomas Schneider
   copyright 1990
   module libraries: delman, delmods, prgmods */

/* end of program */

/* begin module version */

#define version         1.16
/* of aran.p 1994 Sep 5
origin 1990 Oct 2 from alist */
/* end module version */

/* begin module describe.aran */
/*
name
   aran: aligned random sequences

synopsis
   aran(book: in, aranp: in, list: out, sequ: out, output: out)

files
   book: the book generated by Delila

   aranp:  Parameters to control the program.
      The FIRST LINE must contain one real number which is the degree of
      conservation.  For example, if this is 0.85, then each base will have 85%
      chance of being the same, while the other bases will be 5% each.

      The SECOND LINE must contain the number of sequences to generate.

   list: details of the run.

   sequ: the aligned sequences, for input to makebk

   output: messages to the user

description
   Aran takes a sequence as a starting point and generates random
   sequences from it.  The program simulates a very simple dirty synthesis of
   the sequence.  The synthesis is to be mostly the bases given in the
   sequence.  The probability of conserving each base (f) is defined in the
   parameter file.  If a particular base is not conserved, then the other
   three bases are assigned probabilities of (1-f)/3.

documentation
   delman.use.aligned.books

see also
   alist.p

author
   Thomas D. Schneider

bugs
   See alist

technical notes
   The program constant seqmax defines the length of the longest sequence
   that can be created.

*/
/* end module describe.aran */

/* begin module aran.const */

#define seqmax          100   /* maximum size of the sequence array */
/* end module aran.const */

/* begin module book.const */
/* constants needed for book manipulations */

#define dnamax          3000   /* length of dna arrays */
#define namelength      20   /* maximum key name length */
#define linelength      80   /* maximum line readable in book */


/* end module book.const version = 'delmod 6.58 90 Jan 2 tds/gds' */


/* begin module book.type */
/* types needed for book manipulations */

typedef long chset[5];


/* types defined in book definition */

typedef Char alpha[namelength];   /* this is not alfa */

/* name is a left justified string with blanks following the
   characters */

typedef struct name {
  alpha letters;
  /* zero means an unspecified structure */
  char length;
} name;

typedef struct line {
  /* a line of characters */
  Char letters[linelength];
  char length;
  struct line *next;
} line;

typedef enum {
  plus, minus, dircomplement, dirhomologous
} direction;
typedef enum {
  linear, circular
} configuration;
typedef enum {
  on, off
} state;


typedef struct header {
  /* header of key */
  name keynam;   /* key name of structure */
  line *fulnam;   /* full name of structure */
  /* note key */
  line *note;
} header;


/* base types */
typedef enum {
  a, c, g, t
} base;
typedef short dnarange;

/* p2c: aran.p, line 120:
 * Note: Field width for seq assumes enum base has 4 elements [105] */
typedef uchar seq[(dnamax + 3) / 4];

typedef struct dnastring {
  seq part;
  dnarange length;
  struct dnastring *next;
} dnastring;


typedef struct orgkey {
  /* organism key */
  header hea;
  /* genetic map units */
  line *mapunit;
} orgkey;

typedef struct chrkey {
  /* chromosome key */
  header hea;
  double mapbeg;   /* number of genetic map beginning */
  /* number of genetic map ending */
  double mapend;
} chrkey;

typedef struct piekey {
  /* piece key */
  header hea;
  double mapbeg;   /* genetic map beginning */
  configuration coocon;   /* configruation (circular/linear) */
  direction coodir;   /* direction (+/-) relative to genetic map */
  long coobeg;   /* beginning nucleotide */
  long cooend;   /* ending nucleotide */
  configuration piecon;   /* configruation (circular/linear) */
  direction piedir;   /* direction (+/-) relative to coordinates */
  long piebeg;   /* beginning nucleotide */
  long pieend;   /* ending nucleotide */
} piekey;


typedef struct piece {
  piekey key;
  dnastring *dna;
} piece;

typedef struct reference {
  name pienam;   /* name of piece referred to */
  double mapbeg;   /* genetic map beginning */
  direction refdir;   /* direction relative to coordinates */
  long refbeg;   /* beginning nucleotide */
  long refend;   /* ending nucleotide */
} reference;

typedef struct genkey {
  /* gene key */
  header hea;
  reference ref;
} genkey;

typedef struct trakey {
  /* transcript key */
  header hea;
  reference ref;
} trakey;

typedef struct markey {
  /* marker key */
  header hea;
  reference ref;
  state sta;
  line *phenotype;
  struct marker *next;
} markey;

typedef struct marker {
  markey key;
  dnastring *dna;
} marker;

/* end module book.type version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module aran.type */

typedef struct seqmatrix {
  double data[(long)t - (long)a + 1][seqmax];   /* creation matrix */
  long length;   /* portion of the matrix in use */
} seqmatrix;

/* end module aran.type */


Static _TEXT book;   /* the book to be randomized */
Static _TEXT list;   /* details of the run */
Static _TEXT sequ;   /* sequences generated */
Static _TEXT aranp;   /* parameters */

/* begin module book.var */
/* ************************************************************************ */
/* global variables needed for book manipulations */

/* free storage: */
Static line *freeline;   /* unused lines */
Static dnastring *freedna;   /* unused dnas */

Static boolean readnumber;
/* whether to read a number from the notes, or
                        to read in the notes */
Static long number;   /* the number of the item just read */
Static boolean numbered;   /* true when the item just read is numbered */
Static boolean skipunnum;


Static jmp_buf _JL1;


/* a control variable to allow skipping of
                       un-numbered items in the book */

/* ************************************************************************ */
/* end module book.var version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* 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 = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module copyaline */
Static Void copyaline(fin, fout)
_TEXT *fin, *fout;
{
  /* copy a line from file fin to file fout */
  while (!P_eoln(fin->f)) {
    putc(P_peek(fin->f), fout->f);
    getc(fin->f);
  }
  fscanf(fin->f, "%*[^\n]");
  getc(fin->f);
  putc('\n', fout->f);
}  /* copyaline */


/* end module copyaline version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* ************************************************************************ */

/* begin module package.nextbase */
/* ************************************************************************ */
/* begin module package.getpiece */
/* ************************************************************************ */
/* begin module package.brpiece */
/* ************************************************************************ */
/* begin module book.basis */
/* procedures needed for book manipulations */

/* get procedures should be used for all linked lists of records */

Static Void getline(l)
line **l;
{
  /* obtain a line from the free line list or by making a new one */
  if (freeline != NULL) {
    *l = freeline;
    freeline = freeline->next;
  } else
    *l = (line *)Malloc(sizeof(line));
  (*l)->length = 0;
  (*l)->next = NULL;
}


Static Void getdna(l)
dnastring **l;
{
  if (freedna != NULL) {
    *l = freedna;
    freedna = freedna->next;
  } else
    *l = (dnastring *)Malloc(sizeof(dnastring));
  (*l)->length = 0;
  (*l)->next = NULL;
}


/* clear procedures should be called each time the records are no longer needed
         failure to do this may result in a stack overflow. */

Static Void clearline(l)
line **l;
{
  /* return a line to the free line list */
  line *lptr;

  if (*l == NULL)
    return;
  lptr = *l;
  *l = (*l)->next;
  lptr->next = freeline;
  freeline = lptr;
}


Static Void cleardna(l)
dnastring **l;
{
  dnastring *lptr;

  if (*l == NULL)
    return;
  lptr = *l;
  *l = (*l)->next;
  lptr->next = freedna;
  freedna = lptr;
}


Static Void clearheader(h)
header *h;
{
  /* clear the header h (remove lines to free storage) */
  clearline(&h->fulnam);
  while (h->note != NULL)
    clearline(&h->note);
}


Static Void clearpiece(p)
piece **p;
{
  /* clear the dna of the piece */
  while ((*p)->dna != NULL)
    cleardna(&(*p)->dna);
  clearheader(&(*p)->key.hea);
}


Static base chartobase(ch)
Char ch;
{
  /* convert a character into a base */
  base Result;

  switch (ch) {

  case 'a':
    Result = a;
    break;

  case 'c':
    Result = c;
    break;

  case 'g':
    Result = g;
    break;

  case 't':
    Result = t;
    break;
  }
  return Result;
}


Static Char basetochar(ba)
base ba;
{
  /* convert a base into a character */
  Char Result;

  switch (ba) {

  case a:
    Result = 'a';
    break;

  case c:
    Result = 'c';
    break;

  case g:
    Result = 'g';
    break;

  case t:
    Result = 't';
    break;
  }
  return Result;
}


Static base complement(ba)
base ba;
{
  /* take the complement of ba */
  base Result;

  switch (ba) {

  case a:
    Result = t;
    break;

  case c:
    Result = g;
    break;

  case g:
    Result = c;
    break;

  case t:
    Result = a;
    break;
  }
  return Result;
}


Static long pietoint(p, pie)
long p;
piece *pie;
{
  /* p is a coordinate on the piece.
     we want to transform p into a number
     from 1 to n: an internal coordinate system for
     easy manipulation of piece coordinates */
  long i;   /* an intermediate value */
  piekey *WITH;

  WITH = &pie->key;
  switch (WITH->piedir) {

  case plus:
    if (p >= WITH->piebeg)
      i = p - WITH->piebeg + 1;
    else
      i = p - WITH->coobeg + WITH->cooend - WITH->piebeg + 2;
    break;

  case minus:
    if (p <= WITH->piebeg)
      i = WITH->piebeg - p + 1;
    else
      i = WITH->cooend - p + WITH->piebeg - WITH->coobeg + 2;
    break;
  }
  return i;
}


Static long inttopie(i, pie)
long i;
piece *pie;
{
  /* i is in the range 1 to some maximum.  it is an internal coordinate
     system for the program.  we want to do a
     coordinate transformation to obtain
     a value in the range of the piece called pie:
     i=1 corresponds to piebeg and
     i=its maximum corresponds to pieend */
  long p;   /* an intermediate value */
  piekey *WITH;

  WITH = &pie->key;
  switch (WITH->piedir) {

  case plus:
    p = WITH->piebeg + i - 1;
    if (p > WITH->cooend) {
      if (WITH->coocon == circular)
	p += WITH->coobeg - WITH->cooend - 1;
    }
    break;

  case minus:
    p = WITH->piebeg - i + 1;
    if (p < WITH->coobeg) {
      if (WITH->coocon == circular)
	p += WITH->cooend - WITH->coobeg + 1;
    }
    break;
  }
  return p;
}


Static long piecelength(pie)
piece *pie;
{
  /* return the length of the dna in pie */
  return (pietoint(pie->key.pieend, pie));
}


/* end module book.basis version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.getto */
Static Char getto(thefile, ch)
_TEXT *thefile;
long *ch;
{
  /* search the file for a character in the first line which is a
     member of the set ch. */
  Char achar = ' ';

  while ((!P_inset(achar, ch)) & (!BUFEOF(thefile->f))) {
    fscanf(thefile->f, "%c%*[^\n]", &achar);
    getc(thefile->f);
    if (achar == '\n')
      achar = ' ';
  }
  if (P_inset(achar, ch))
    return achar;
  else
    return ' ';
}


/* end module book.getto version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.skipstar */
Static Void skipstar(thefile)
_TEXT *thefile;
{
  /* skip start of line (or star = '*'). */
  if (P_peek(thefile->f) != '*') {
    printf(" procedure skipstar: bad book\n");
    printf(" \"*\" expected as first character on the line, but \"%c\" was found\n",
	   P_peek(thefile->f));
    halt();
  }
  getc(thefile->f);   /* skip the star */

  if (P_peek(thefile->f) != ' ') {   /* skip the blank */
    printf(" procedure skipstar: bad book\n");
    printf(" \"* \" expected on a line but \"*%c\" was found\n",
	   P_peek(thefile->f));
    halt();
  }
  getc(thefile->f);
}  /* skipstar */


/* end module book.skipstar version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brreanum */
Static Void brreanum(thefile, reanum)
_TEXT *thefile;
double *reanum;
{
  /* read a real number from the file */
  skipstar(thefile);
  fscanf(thefile->f, "%lg%*[^\n]", reanum);
  getc(thefile->f);
}


/* end module book.brreanum version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brnumber */
Static Void brnumber(thefile, num)
_TEXT *thefile;
long *num;
{
  /* read a number from the file */
  skipstar(thefile);
  fscanf(thefile->f, "%ld%*[^\n]", num);
  getc(thefile->f);
}


/* end module book.brnumber version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brname */
Static Void brname(thefile, nam)
_TEXT *thefile;
name *nam;
{
  /* read a name from the file */
  long i;   /* an index to the name */
  Char c_;   /* a character read */

  skipstar(thefile);
  nam->length = 0;
  do {
    nam->length++;
    c_ = getc(thefile->f);
    if (c_ == '\n')
      c_ = ' ';
    nam->letters[nam->length - 1] = c_;
  } while (!(P_eoln(thefile->f) || nam->length >= namelength ||
	     nam->letters[nam->length - 1] == ' '));
  if (nam->letters[nam->length - 1] == ' ')
    nam->length--;
  if (nam->length < namelength) {
    for (i = nam->length; i < namelength; i++)
      nam->letters[i] = ' ';
  }
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
}  /* brname */


/* end module book.brname version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brline */
Static Void brline(thefile, l)
_TEXT *thefile;
line **l;
{
  /* read a line from the file */
  long i = 0;
  long j;
  Char acharacter;
  long FORLIM;

  skipstar(thefile);
  while (!P_eoln(thefile->f)) {
    i++;
    acharacter = getc(thefile->f);
    if (acharacter == '\n')
      acharacter = ' ';
    (*l)->letters[i-1] = acharacter;
  }
  if (i < (*l)->length) {
    FORLIM = (*l)->length;
    for (j = i; j < FORLIM; j++)
      (*l)->letters[j] = ' ';
  }
  (*l)->length = i;
  (*l)->next = NULL;
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
}


/* end module book.brline version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brdirect */
Static Void brdirect(thefile, direct)
_TEXT *thefile;
direction *direct;
{
  /* read a direction */
  Char ch;

  skipstar(thefile);
  fscanf(thefile->f, "%c%*[^\n]", &ch);
  getc(thefile->f);
  if (ch == '\n')
    ch = ' ';
  if (ch == '+')
    *direct = plus;
  else
    *direct = minus;
}


/* end module book.brdirect version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brconfig */
Static Void brconfig(thefile, config)
_TEXT *thefile;
configuration *config;
{
  /* read a configuration */
  Char ch;

  skipstar(thefile);
  fscanf(thefile->f, "%c%*[^\n]", &ch);
  getc(thefile->f);
  if (ch == '\n')
    ch = ' ';
  if (ch == 'l')
    *config = linear;
  else
    *config = circular;
}


/* end module book.brconfig version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brnotenumber */
Static Void brnotenumber(thefile, note)
_TEXT *thefile;
line **note;
{
  /* book note reading to obtain the number of the object.
the procedure returns the value of the number as a global.
(this is not such a good practice, but we are stuck with it for now.) */
  *note = NULL;
  numbered = false;
  number = 0;
  /* force number to zero if there
     is no number at all */
  /* the next character is n or * depending on whether there are notes */
  if (P_peek(thefile->f) != 'n')
    return;
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
  if (P_peek(thefile->f) == 'n') {
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);
    return;
  }
  skipstar(thefile);
  if (!P_eoln(thefile->f)) {
    if (P_peek(thefile->f) == '#') {
      numbered = true;
      getc(thefile->f);   /* move past the number symbol */
      fscanf(thefile->f, "%ld", &number);
    }
  }
  do {
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);
  } while (P_peek(thefile->f) != 'n');
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
}  /* brnotenumber */


/* end module book.brnotenumber version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brnote */
Static Void brnote(thefile, note)
_TEXT *thefile;
line **note;
{
  /* read note key */
  line *newnote;   /* the new note */
  line *previousnote;   /* the last line of the notes */

  *note = NULL;
  if (P_peek(thefile->f) != 'n')  /* enter note */
    return;
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
  if (P_peek(thefile->f) == 'n') {  /* abort null note (n/n) */
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);
    return;
  }
  getline(note);
  newnote = *note;
  while (P_peek(thefile->f) != 'n') {  /* wait until end of note */
    brline(thefile, &newnote);
    previousnote = newnote;
    /* get next note */
    getline(&newnote->next);
    newnote = newnote->next;
  }
  /* last note was not used, so: */
  clearline(&newnote);
  previousnote->next = NULL;
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
}  /* brnote */


/* end module book.brnote version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brheader */
Static Void brheader(thefile, hea)
_TEXT *thefile;
header *hea;
{
  /* read the header of a key. */
  /* read key name */
  brname(thefile, &hea->keynam);

  /* read full name */
  getline(&hea->fulnam);
  brline(thefile, &hea->fulnam);

  /* read note key */
  if (readnumber)
    brnotenumber(thefile, &hea->note);
  else
    brnote(thefile, &hea->note);
}


/* end module book.brheader version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brpiekey */
Static Void brpiekey(thefile, pie)
_TEXT *thefile;
piekey *pie;
{
  /* read piece key */
  brheader(thefile, &pie->hea);
  brreanum(thefile, &pie->mapbeg);
  brconfig(thefile, &pie->coocon);
  brdirect(thefile, &pie->coodir);
  brnumber(thefile, &pie->coobeg);
  brnumber(thefile, &pie->cooend);
  brconfig(thefile, &pie->piecon);
  brdirect(thefile, &pie->piedir);
  brnumber(thefile, &pie->piebeg);
  brnumber(thefile, &pie->pieend);
}


/* end module book.brpiekey version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brdna */
Static Void brdna(thefile, dna)
_TEXT *thefile;
dnastring **dna;
{
  /* read in dna from thefile */
  /* note: if the dna were circularized, by linking the last dnastring
to the first, then the cleardna routine could not clear properly,
and would loop forever... there is no reason to do that, since a simple
mod function will allow one to access the circle. */
  Char ch;
  dnastring *workdna;
  long SET[5];
  long TEMP;

  getdna(dna);
  workdna = *dna;
  ch = getto(thefile, P_addset(P_expset(SET, 0L), 'd'));
  ch = getc(thefile->f);   /* skipstar */
  if (ch == '\n')
    ch = ' ';
  while (ch == '*') {
    ch = getc(thefile->f);   /* skip blank */
    if (ch == '\n')
      ch = ' ';
    do {
      ch = getc(thefile->f);
      if (ch == '\n')
	ch = ' ';
      if (ch == 't' || ch == 'g' || ch == 'c' || ch == 'a') {
	if (workdna->length == dnamax) {
	  getdna(&workdna->next);
	  workdna = workdna->next;
	}
	workdna->length++;
	TEMP = workdna->length - 1;
	P_clrbits_B(workdna->part, TEMP, 1, 3);
	P_putbits_UB(workdna->part, TEMP, (int)chartobase(ch), 1, 3);
      }
    } while (!P_eoln(thefile->f));
    fscanf(thefile->f, "%*[^\n]");
    getc(thefile->f);   /* go to next line */
    ch = getc(thefile->f);   /* ch is either '*' or 'd' */
    if (ch == '\n')
      ch = ' ';
  }
  fscanf(thefile->f, "%*[^\n]");
  getc(thefile->f);
}


/* end module book.brdna version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brpiece */
Static Void brpiece(thefile, pie)
_TEXT *thefile;
piece **pie;
{
  /* read in a piece */
  brpiekey(thefile, &(*pie)->key);
  if (numbered || !skipunnum)
    brdna(thefile, &(*pie)->dna);
}


/* end module book.brpiece version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.brinit */
Static Void brinit(book)
_TEXT *book;
{
  /* check that the book is ok to read, and
     set up the global variables for br routines */
  /* halt if the book is bad (first word is 'halt') or the first
     character is not * */
  if (*book->name != '\0') {
    if (book->f != NULL)
      book->f = freopen(book->name, "r", book->f);
    else
      book->f = fopen(book->name, "r");
  } else
    rewind(book->f);
  if (book->f == NULL)
    _EscIO2(FileNotFound, book->name);
  RESETBUF(book->f, Char);
  if (!BUFEOF(book->f)) {
    /* check for the date line */
    if (P_peek(book->f) != '*') {
      if (P_peek(book->f) != 'h')
	printf(" this is not the first line of a book:\n");
      else
	printf(" bad book:\n");
      putchar(' ');

      while (!(P_eoln(book->f) | BUFEOF(book->f))) {
	putchar(P_peek(book->f));
	getc(book->f);
      }
      putchar('\n');
      halt();
    }
  } else {
    printf(" book is empty\n");
    halt();
  }

  /* initialize free storage */
  freeline = NULL;
  freedna = NULL;

  readnumber = true;   /* usually we read in numbers for items */
  number = 0;   /* arbitrary value */
  numbered = false;   /* the piece has no number (none yet read in) */
  skipunnum = false;
}  /* brinit */


/* end module book.brinit version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* ************************************************************************ */
/* end module package.brpiece version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.getpiece */
Static Void getpiece(thefile, pie)
_TEXT *thefile;
piece **pie;
{
  /* move to and read in the next piece in the book */
  Char ch;
  long SET[5];

  ch = getto(thefile, P_addset(P_expset(SET, 0L), 'p'));
      /* get to the next p(iece) in the book */
  if (ch != ' ') {
    brpiece(thefile, pie);
    ch = getto(thefile, P_addset(P_expset(SET, 0L), 'p'));
	/* read past closing p */
  }
}


/* end module book.getpiece version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* ************************************************************************ */
/* end module package.getpiece version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module book.stepbase */
Static base stepbase(startdna, dna, d)
dnastring *startdna, **dna;
dnarange *d;
{
  /* advance d by one base in dna and then return the base at the new d.
     (this means that one should initialize d to zero)
     if we go past the last base, we restart at startdna.
     note: d is not the number of the base... it is used as a
     record for stepbase.  do not mess with it, and do not use it to find
     out what base you are on.  use a separate counter.  */
  long TEMP;

  if (*d != dnamax && *d != (*dna)->length) {
    (*d)++;
    TEMP = *d - 1;
    return ((base)P_getbits_UB((*dna)->part, TEMP, 1, 3));
  }
  *d = 1;
  *dna = (*dna)->next;
  if (*dna == NULL)
    *dna = startdna;
  TEMP = *d - 1;
  return ((base)P_getbits_UB((*dna)->part, TEMP, 1, 3));
}


/* end module book.stepbase version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module nextbase */
Static Void initnextbase(book, pie, lastbase, endofbook)
_TEXT *book;
piece **pie;
boolean *lastbase, *endofbook;
{
  /* initialize variables for function nextbase.
     book is the book to be read,
     pie is the piece,
     lastbase is the flag that is true when we are at the
     last base of a piece,
     and endofbook is true if we are at the end of the book (see nextbase). */
  piece *WITH;
  header *WITH1;

  brinit(book);

  *pie = (piece *)Malloc(sizeof(piece));
  WITH = *pie;

  WITH1 = &WITH->key.hea;
  WITH1->fulnam = NULL;
  WITH1->note = NULL;
  WITH->dna = NULL;
  *lastbase = true;   /* this will trigger reading of the next piece */
  if (BUFEOF(book->f))
    *endofbook = true;
  else
    *endofbook = false;
}


Static base nextbase(book, pie, dnalink, dnalinkspot, dnaspot, length,
		     lastbase, endofbook)
_TEXT *book;
piece **pie;
dnastring **dnalink;
dnarange *dnalinkspot;
long *dnaspot, *length;
boolean *lastbase, *endofbook;
{
  /* the book being read */
  /* the next variables can be ignored */
  /* the piece */
  /* the current link we are on */
  /* the spot in the dnalink */
  /* these are useful to the general user: */
  /* integer in 1 to length, which base this is */
  /* length of this piece */
  /* true if the base was the last one on
                   the piece. */
  /* true when we have reached
                         the end of the book */
  /* the user simply declares variables for
     book, pie, dnalink, dnalinkspot (of the appropriate type)
note: you can convert the valuespot to published coordinates by
        pubcoords:= inttopie(dnaspot,pie);

warning: if the end of the book has been reached, then endofbook
is true, but the value returned by the function has no meaning.
  */
  base Result;

  if (!*lastbase) {
    Result = stepbase((*pie)->dna, dnalink, dnalinkspot);
    (*dnaspot)++;
    if (*dnaspot == *length)
      *lastbase = true;
    return Result;
  } else {
    clearpiece(pie);
    getpiece(book, pie);
    if (!BUFEOF(book->f)) {
      *dnalink = (*pie)->dna;
      *dnalinkspot = 0;
      *dnaspot = 0;
      *length = piecelength(*pie);
      *lastbase = false;
      *endofbook = false;
      return (nextbase(book, pie, dnalink, dnalinkspot, dnaspot, length,
		       lastbase, endofbook));
    } else {   /* we are no longer at the last base */
      *endofbook = true;   /* we have reached the end of the book */
      *lastbase = false;
      return a;   /* a fake value */
    }
  }
  /* we are at the last base of the previous piece */
  return Result;
}


#define pow14           16384

#define pow15           32768L
#define pow22           4194304L
#define pow23           8388608L


/* end module nextbase version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* ************************************************************************ */
/* end module package.nextbase version = 'delmod 6.58 90 Jan 2 tds/gds' */

/* begin module random */
Static Void random(seed)
double *seed;
{
  /* random generator 2.
     version = 1.01 of random.2 1990 Oct 2
     origin 1986 December 31
     Test this routine with the program tstrnd.
     written by David Masternarde */

  /* This random number generator is based on a shift register with a single bit
of feedback, as described in Electronics for Neurobiologists, by Brown, Maxfield
and Moraff, MIT press 1973, referencing Random Process Simulation and
Measurement by Korn, McGraw-Hill 1966.  The random seed rand, a number between 0
and 1 exclusive, is converted to an integer between 1 and 2**23-1, inclusive.
This 23-bit number is shifted right one bit and the output of the last (23rd)
bit and the 9th bit are added modulo 2 (exclusive orred) and fed back into the
new first bit.  This is done between 4 and 11 times, depending on the last 3
bits of the original number.  The result is converted back to a real number
between 0 and 1 from which the 23 bit integer can be recovered on the next call.
The 23-bit shift register goes through all 2**23-1 values before repeating;
the repetition frequency of this algorithm could be less or greater depending
on the seed, because of the random number of multiple shifts per call. */
  /* powers of 2 */
  long iseed;   /* integer shift register */
  long i, nrep;   /* index, number of times to do shift */

  iseed = (long)floor(*seed * pow23 + 0.5);   /* convert to 23 bit number */
  if (iseed < 1 || iseed >= pow23)
    iseed = 1;
  nrep = (iseed & 7) + 4;   /* do it 4 to 11 times based on last 3 bits */
  for (i = 1; i <= nrep; i++) {
    /* if last bit and 9th bit are equal, feed back a 0, otherwise a 1 */
    if ((iseed & 1) == ((iseed & (pow15 - 1)) >= pow14))
      iseed /= 2;
    else
      iseed = iseed / 2 + pow22;
  }
  *seed = (double)iseed / pow23;
}  /* random */

#undef pow14
#undef pow15
#undef pow22
#undef pow23


/* end module random */

/* begin module aran.readparameters */
Static Void readparameters(aranp, fraction, number)
_TEXT *aranp;
double *fraction;
long *number;
{
  /* read the parameters */
  if (*aranp->name != '\0') {
    if (aranp->f != NULL)
      aranp->f = freopen(aranp->name, "r", aranp->f);
    else
      aranp->f = fopen(aranp->name, "r");
  } else
    rewind(aranp->f);
  if (aranp->f == NULL)
    _EscIO2(FileNotFound, aranp->name);
  RESETBUF(aranp->f, Char);
  if (BUFEOF(aranp->f)) {
    printf("aranp must contain two parmeters\n");
    halt();
  }
  fscanf(aranp->f, "%lg%*[^\n]", fraction);
  getc(aranp->f);

  if ((unsigned)(*fraction) > 1.0) {
    printf("The fraction given in aranp must be between 0 and 1\n");
    halt();
  }

  if (BUFEOF(aranp->f)) {
    printf("aranp must contain the number of sequences to generate\n");
    halt();
  }

  fscanf(aranp->f, "%ld%*[^\n]", number);
  getc(aranp->f);

  if (*number < 1) {
    printf("at least one sequence must be generated\n");
    halt();
  }
}


/* end module aran.readparameters */

/* begin module aran.getseqmatrix */
Static Void getseqmatrix(book, fout, conserved, unconserved, sm)
_TEXT *book, *fout;
double conserved, unconserved;
seqmatrix *sm;
{
  /* Get the sequence matrix.  The procedure is a modification of the
demonextbase routine in delmod.p */
  base b;   /* a base of the sequence and index to the sm matrix */
  long l;   /* index to the sm matrix */

  /* these variables are all defined in nextbase */
  piece *pie;
  dnastring *dnalink;
  dnarange dnalinkspot;
  long dnaspot, length;
  boolean lastbase, endofbook;

  /* this variable is needed to catch the value of nextbase,
     so that we can check that we have not reached the end of the book. */
  Char character;
  long FORLIM;

  /* clear the matrix completely */
  for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) {
    for (l = 0; l < seqmax; l++)
      sm->data[(long)b - (long)a][l] = unconserved;
  }

  l = 0;
  /* obtain the sequence, and store it in the matrix */
  initnextbase(book, &pie, &lastbase, &endofbook);
  do {
    b = nextbase(book, &pie, &dnalink, &dnalinkspot, &dnaspot, &length,
		 &lastbase, &endofbook);
    if (length > seqmax) {
      printf("sequence length (%ld) exceeds seqmax constant (%ld)\n",
	     length, (long)seqmax);
      halt();
    }
    l++;
    sm->data[(long)b - (long)a][l-1] = conserved;
	/* insert the base into the matrix */

    character = basetochar(b);
    putc(character, fout->f);
    /*
          if not endofbook
             then writeln(fout,' ',character,' ',dnaspot:1,
                                 ' ',inttopie(dnaspot,pie):1);
          if lastbase
             then writeln(fout,' end of a piece ',length:1,'bp')
    */
  } while (!lastbase);

  sm->length = l;
  fprintf(fout->f, ".\n");
  fprintf(fout->f, "%ld bases\n", sm->length);

  fprintf(fout->f, "----\n");

  fprintf(fout->f, "The frequency matrix:\n");
  FORLIM = sm->length;
  for (l = 0; l < FORLIM; l++) {
    for (b = a; (long)b <= (long)t; b = (base)((long)b + 1))
      fprintf(fout->f, " %4.2f", sm->data[(long)b - (long)a][l]);
    putc('\n', fout->f);
  }
  fprintf(fout->f, "----\n");

  /* now sum across the matrix so things go faster during generation */
  fprintf(fout->f, "Sum across the matrix for faster sequence generation:\n");
  FORLIM = sm->length;
  for (l = 0; l < FORLIM; l++) {
    for (b = c; (long)b <= (long)t; b = (base)((long)b + 1))
      sm->data[(long)b - (long)a][l] += sm->data[(long)b - (long)a - 1][l];
  }

  FORLIM = sm->length;
  for (l = 0; l < FORLIM; l++) {
    for (b = a; (long)b <= (long)t; b = (base)((long)b + 1))
      fprintf(fout->f, " %4.2f", sm->data[(long)b - (long)a][l]);
    putc('\n', fout->f);
  }
  fprintf(fout->f, "----\n");
}  /* getseqmatrix */


#define pagewidth       60   /* width of a page for the sequence output */


/* end module aran.getseqmatrix */

/* begin module aran.generate */
Static Void generate(seed, sm, number, sequ)
double *seed;
seqmatrix sm;
long number;
_TEXT *sequ;
{
  /* generate a number of sequences to the sequ file, given the
matrix sm.  seed is a random number between 0 and 1 */
  long l;   /* index to the sm matrix */
  long n;   /* index to sequences */

  if (*sequ->name != '\0') {
    if (sequ->f != NULL)
      sequ->f = freopen(sequ->name, "w", sequ->f);
    else
      sequ->f = fopen(sequ->name, "w");
  } else {
    if (sequ->f != NULL)
      rewind(sequ->f);
    else
      sequ->f = tmpfile();
  }
  if (sequ->f == NULL)
    _EscIO2(FileNotFound, sequ->name);
  SETUPBUF(sequ->f, Char);
  for (n = 1; n <= number; n++) {
    for (l = 1; l <= sm.length; l++) {
      if (l % pagewidth == 0)
	putc('\n', sequ->f);
/* p2c: aran.p, line 1010:
 * Note: Using % for possibly-negative arguments [317] */
      random(seed);   /* get a random number */
      if (sm.data[0][l-1] > *seed)
	putc('a', sequ->f);
      else if (sm.data[(long)c - (long)a][l-1] > *seed)
	putc('c', sequ->f);
      else if (sm.data[(long)g - (long)a][l-1] > *seed)
	putc('g', sequ->f);
      else
	putc('t', sequ->f);   /* must be */
    }
    fprintf(sequ->f, ".\n");
  }
}

#undef pagewidth


/* end module aran.generate */

/* begin module aran.themain */
Static Void themain(book, list, sequ, aranp)
_TEXT *book, *list, *sequ, *aranp;
{
  /* the main procedure of the program */
  double conserved;   /* fraction of conserved bases */
  long number;   /* number of sequences to generate */
  double unconserved;   /* fraction of each base not conserved */
  double seed = 0.5;   /* a random number to start the random series */
  seqmatrix sm;   /* the matrix for generating the sequences */

  printf("aran %4.2f\n", version);
  if (*sequ->name != '\0') {
    if (sequ->f != NULL)
      sequ->f = freopen(sequ->name, "w", sequ->f);
    else
      sequ->f = fopen(sequ->name, "w");
  } else {
    if (sequ->f != NULL)
      rewind(sequ->f);
    else
      sequ->f = tmpfile();
  }
  if (sequ->f == NULL)
    _EscIO2(FileNotFound, sequ->name);
  SETUPBUF(sequ->f, Char);
  readparameters(aranp, &conserved, &number);
  unconserved = (1.0 - conserved) / 3.0;

  if (*list->name != '\0') {
    if (list->f != NULL)
      list->f = freopen(list->name, "w", list->f);
    else
      list->f = fopen(list->name, "w");
  } else {
    if (list->f != NULL)
      rewind(list->f);
    else
      list->f = tmpfile();
  }
  if (list->f == NULL)
    _EscIO2(FileNotFound, list->name);
  SETUPBUF(list->f, Char);
  fprintf(list->f, "aran %4.2f\n", version);
  fprintf(list->f, "  Conserved fraction: %4.2f\n", conserved);
  fprintf(list->f, "unConserved fraction: %4.2f\n", unconserved);
  getseqmatrix(book, list, conserved, unconserved, &sm);
  fprintf(list->f, "%ld sequences generated\n", number);

  generate(&seed, sm, number, sequ);
}


/* end module aran.themain */

main(argc, argv)
int argc;
Char *argv[];
{
  PASCAL_MAIN(argc, argv);
  if (setjmp(_JL1))
    goto _L1;
  aranp.f = NULL;
  strcpy(aranp.name, "aranp");
  sequ.f = NULL;
  strcpy(sequ.name, "sequ");
  list.f = NULL;
  strcpy(list.name, "list");
  book.f = NULL;
  strcpy(book.name, "book");
  themain(&book, &list, &sequ, &aranp);
_L1:
  if (book.f != NULL)
    fclose(book.f);
  if (list.f != NULL)
    fclose(list.f);
  if (sequ.f != NULL)
    fclose(sequ.f);
  if (aranp.f != NULL)
    fclose(aranp.f);
  exit(EXIT_SUCCESS);
}  /* aran */



/* End. */