/* Output from p2c 1.21alpha-07.Dec.93, the Pascal-to-C translator */ /* From input file "evd.p" */ #include /* evolution display Dr. Thomas D. Schneider National Cancer Institute Laboratory of Experimental and Computational Biology Frederick, Maryland 21702-1201 toms@ncifcrf.gov permanent email: toms@alum.mit.edu http://www.lecb.ncifcrf.gov/~toms/ module libraries required: delman, prgmods, delmods, ev */ /* end of program */ /* begin module version */ #define version 2.48 /* of evd.p 2005 Jul 20 2005 Mar 23, 2.45: upgrade for GPC 2.43, 2002 Apr 4: compute thermodynamic information of matrix 2.42, 2002 Apr 4: make gpc compatable: variable x local to dononsites 2.41, 2002 Mar 28: remove namelength const and alpha type which were not used 2.38, 2000 Jul 18: invent ev.run script 2.37, 2000 Jul 17: add Schneider.oxyr ref. 2.35, 2000 Jul 16: upgrade documentation 2.34, 2000 Jul 16: upgrade documentation 2.31, 2000 Jul 9: upgrade documentation for publication release 2.30, 1999 Jul 26: need to evaluate if selecting was off for proper display 2.27, 1999 Jul 23: (+++) now have parens between bases 2.23, 1999 Jul 22: tidy up; fix location of threshold in evfeatures 2.21, 1999 Feb 18: upgrade to gaussian noise 2.20, 1999 Feb 17: upgrade to more modern time modules, solve y2k problem 1998 Mar 7 upgrade to allow random placement of sites 1988 oct 27 previous changes 1984 apr 19 origin */ #define updateversion 2.26 /* defines lowest acceptable current parameter file */ /* end module version */ /* begin module describe.evd */ /* name evd: evolution display synopsis evd(all: in, evdp: in, display: out, sites: out, genomes: out, evfeatures: out, evdp: in, output: out); files all: the all file from program ev. It contains all the genomes of the evolving creatures, and other parameters and data. It is created by a binary dump by the ev program for the sake of speed and so is not readable by people and so probably cannot be successfully run on another kind of computer. The evd program therefore is needed to interpret the all file. evdp: parameters. one per line: firstcreature: the number of the first creature to display secondcreature: the number of the last creature to display non-site features: if the first character is 'n' then non-sites that are recognized by the weight matrix are shown as features in evfeatures If evdp is empty, the defaults are: 1/1/- (first creature only, no non-site features) The creature of rank 1 makes the fewest mistakes; number 2 makes more, etc. display: a marked display of the genomes and other data. sites: raw sequences of the sites (and 5 bases around each), the sequences are separated by periods, and different creatures are separated by blank lines. The current method for using this to create a sequence logo is described in the 'see also' section. genomes: raw sequences of the genomes of the creatures, separated by periods. evfeatures: The features in this genome in the form used by the lister program. output: messages to the user. description The purpose of the ev program is to evolve sites; the purpose of evd is to display an intermediate or final result of an evolutionary run. The genomes are displayed with the locations of the recognizer gene and its sites: a-- c-- g-- t-- is the region encoding the weights for one of the recognizer fingers. The numbers underneath are the weights for each base. TTTT marks the threshold for the matrix. The number underneath is the threshold. (-------------) is a site. The number underneath is the evaluation of the site by the weight matrix. If the site is recognized (the evaluation is greater or equal to the threshold) then + signs are used. The sequences around each site are written to the file 'sites', and the entire genomes are written to the file 'genomes'. This allows analysis by other programs. EXPERIMENTAL: Thermodynamic Probabilities and Information of Weight Matrix Thermodynamic probabilities are computed by assuming that the weight matrix values represent energies, which may not be true. However, given this assumption, we can compute the corresponding probabilities in a Boltzmann distribution by first computing the partition function Q = sum_b exp(weight_b) where b is one of the four bases and weight_b is the weight at some position l for base b. Then the probability for base b is exp(weight_b)/Q. The uncertainty and information are then computed in the normal way. As a technical note, weight values are often high integers, such as 400. This will exceed the capability of the exponentiation. To avoid this, the absolute value of the weights at all positions are taken and the largest weight is used to normalize the entire matrix to the range -1 to +1. 2002 April 4: The thermodynamic computation was implemented. Notably, for the standard evp.selection at 1000 generations, where Rf = 4.0 and Rs = 4.71035+/-0.29733, the computed value is 0.98 bits. I do not yet know what this discrepancy means. However, if the normalization is set to be between -0.5 and +0.5, then the computed value increases to 3.08 bits. This may indicate that the computation is meaningless or that something else has to be done ... see also {The "Evolution of Biological Information" paper (with active hypertext links in references):} http://www-lecb.ncifcrf.gov/~toms/paper/ev/ {Example parameter file:} evdp {Program for the evolution of binding sites (this creates the all file):} ev.p {Program to display the genomes marked by the sites:} {* individual information version:} lister.p {* public version:} listerx.p {These programs require the Delila book format rather than the simple sequence in the genomes file. To convert, use the} makebk.p {program: cp genomes sequ # copy the sites file to the sequ file makebk <} makebkp {# make the delila book with} makebkp {and} makebk.p http://www.lecb.ncifcrf.gov/~toms/icons/donor.pure.gif {TO MAKE LOGOS: Briefly, after running ev to evolve binding sites, use this} evd.p {program to get the binding site sequences. Then copy the sites file to the sequ file and use the} makebk.p {program (preferably in automatic mode) to create a Delila book. Then use standard delila system programs to create the logo, as described at} http://www.lecb.ncifcrf.gov/~toms/delila.html{.} {Assuming that you have all the necessary input files, in detail the procedure is: ev # evolve the creatures with} ev.p { evd # make the display files with this program,} evd.p { cp sites sequ # copy the sites file to the sequ file (a unix command) makebk <} makebkp {# make the delila book with} makebkp {and} makebk.p { encode # be sure to use the f mode in} encodep {with} encode.p { rseq # compute information using the} rseq.p {program dalvec # make the symvec using the} dalvec.p {program makelogo # make the logo from the symvec with} makelogo.p{. Descriptions of the required input files are given in the documentation of each program and a general review is given in the paper} http://www-lecb.ncifcrf.gov/~toms/paper/oxyr/{. An example script that performs the steps above and also creates all the necessary input files is} run.ev{. {A movie of binding site evolution created using these steps is at} http://www.lecb.ncifcrf.gov/~toms/paper/ev/movie/ author Thomas Dana Schneider bugs none known (well, actually it's full of them... ;-) Evd cannot handle non-random site placement in the display file because the drawing mechanism was not designed to do so. The lister or listerx program has to be used. A warning is provided to output. technical notes */ /* end module describe.evd */ /* a handy line for indicating that the program is under construction: http://www.lecb.ncifcrf.gov/~toms/icons/construction.gif {Under Construction} */ /* begin module ev.const */ #define maxgenomesize 5000 /* was 1040 maximum genome size */ /* 1048576 gives segmentation fault */ #define maxgamma 128 /* maximum number of sites. this is used as the array bound on the locations of the sites */ #define maxbugs 128 /* maximum number of creatures */ #define maxwidth 20 /* maximum width of a recognizer */ #define maxbpi 10 /* one more than the maximum bases per integer. this value determines how fine the fingers of the recognizer are able to feel */ #define extra 5 /* the number of extra bases around each site printed into file sequ. */ #define col 14 /* column width for the data display */ #define dec 5 /* decimal places for the data display */ #define displayingtypes 8 /* number of ways to display things in list */ #define mutatehalf false /* if true, only newly replicated bugs are mutated, thereby preserving good strains. If false, all bugs are mutated every generation. */ /* end module ev.const version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module evd.const */ #define deflinewidth 30 /* bases per line, maximum */ /* end module evd.const */ /* begin module datetime.const */ #define datetimearraylength 19 /* length of dataarray for dates, It is just long enough to include the 4 digit year - solving the year 2000 problem: 1980/06/09 18:49:11 123456789 123456789 1 2 */ /* end module datetime.const version = 'cdatemod.p 1.19 1999Dec13'; */ /* begin module datetime.type */ /* array for dates */ typedef Char datetimearray[datetimearraylength]; /* end module datetime.type version = 'cdatemod.p 1.19 1999Dec13'; */ /* begin module ev.type */ typedef enum { a, c, g, t } base; /* the symbols in the genomes of the bugs */ typedef struct misc { /* various precalculated data */ long power[maxbpi][(long)t - (long)a + 1]; /* powers of 4. power[x,a] = 0 power[x,c] = 4 ^ x. power[x,g] = 2*power[x,c] power[x,t] = 3*power[x,c] this forms a lookup-table for quick translation */ long negative; /* the smallest value in the range 0..4^bpi that represents a negative value in two s complement notation */ long twonegative; /* twice negative. this is the actual amount to be subtracted */ long endscang; /* g-1; the end of a scan of the genome */ double ehnb; /* the sampling bias factor for the number of sites, gamma */ double varhnb; /* the variance of ehnb */ double rstdev; /* the total standard deviation of the sampling bias, sqrt(varhnb*width) */ double hg; /* genomic uncertainty */ double eofnforhg; /* e(n) for Hg */ double ln2; /* the natural log of 2 for calculating things in bits */ double rfrequency; /* the information needed to locate sites */ double evversion; /* capture the version number of ev */ } misc; /* the genome of a creature is set up as a record to allow indexing to one bug. this should be faster. */ /* p2c: evd.p, line 286: Note: * Field width for chromosome assumes enum base has 4 elements [105] */ typedef uchar chromosome[(maxgenomesize + 3) / 4]; typedef struct genometype { chromosome genome; /* current composition of the genome: */ long composition[(long)t - (long)a + 1]; } genometype; /* the mistakes each creature has made. The first dimension is the rank number. The second dimension has: 1. bug identification number 2. mistakes made. This allows sorting using quicksort. */ typedef uchar position; /* the position type for quicksort */ typedef enum { bugid, mistakes } rankthings; /* ranktype = array[1..maxbugs,bugid..mistakes] of integer; uuu */ typedef long ranktype[maxbugs][2]; typedef struct population { /* information about the entire population */ long bugs; /* number of creatures */ /* genome: integer; (* number of potential binding sites *) */ long Genome; /* number of potential binding sites. 2005 Jun 15. This was formerly called 'genome' but that is not a very good name. 'G' is worse though because it conflicts with 'g' in the genome. So use 'Genome'. */ long gamma; /* number of sites */ Char randomsites; /* r: random site placement with overlap n: random site placement, no overlap a: array placement */ long width; /* width of sites in bases */ long bpi; /* bases per integer in recognizer gene wmatrix */ long bpthreshold; /* bases per threshold integer */ long genomesize; /* the total length of the genome, calculated as g+width-1+extra */ long site[maxgamma]; /* locations of the sites in the genomes. The first base of a site that starts at base x is at base x+1. See function recognize. */ genometype creature[maxbugs]; /* all the bugs */ ranktype rank; /* the rankings of the bugs by their mistakes */ /* this array records the locations of the sites to make testing faster */ uchar sitelocations[(maxgenomesize + 7) / 8]; long halfofbugs; /* half of the bugs, (bugs div 2) */ long firstmutated; /* the first bug mutated. This is either 1 or halfofbugs+1 */ } population; typedef struct everything { misc precalc; /* lots of precalculated stuff */ datetimearray datetimestamp; /* the date and time used to identify the moment of creation on all output */ long mu; /* mutation rate in hits per creature per generation */ double seed; /* random number generator seed */ double initialseed; /* the seed read from the evp */ population p; long generation; /* cycle number */ long cycles; /* cycles to run */ long displayinterval; /* every display-th generation is shown */ /* NB: if more types are added to this list, increase the constant displayingtypes in module ev.const */ boolean displayaverage; /* display average mistakes */ boolean displaychange; /* display change of mistakes */ boolean displaygenomich; /* display genomic uncertainty */ boolean displayindividual; /* display individual mistakes */ boolean displayorbits; /* display Rindividual orbits */ boolean displayrsequence; /* display Rsequence */ boolean displaystatus; /* display to output */ boolean displayonemistake; /* display the best individuals mistakes */ long previousmistakes; /* number of mistakes made in previous display (not in previous generation, it's only updated when the display is activated according to displayinterval */ boolean selecting; /* If true, then the organisms are sorted by their mistakes. If false, then the organisms are randomly sorted. Normally this should be 'true'. */ long storagefrequency; /* frequency (generations) with which to store the everything record. */ double sigma; /* standard deviation of noise */ boolean noise; /* true if sigma > 0; used for speeding the code */ boolean haveone; /* have a gaussian value; used by the gaussian procedure to switch between gaussianX and gaussianY */ double gaussianX, gaussianY; /* gaussian values for the gaussian procedure */ Char specialrule; /* SPECIAL RULE for ties: r - random survival of one organism b - both survive (original method) s - one survives depending on the sorting algorithm. This is an arbitrary function of how the quicksort routine works and is therefore generally not a good way to make the decision. */ } everything; /* the allfile allows the program to be stopped and restarted again: */ typedef struct allfile { FILE *f; FILEBUFNC(f,everything); Char name[_FNSIZE]; } allfile; typedef struct wmatrix { long matrix[(long)t - (long)a + 1][maxwidth]; /* the w matrix */ /* zzz consider switching to zero based method: matrix: array[a..t, 0..maxwidth] of integer; (* the w matrix *) */ long threshold; /* the lower threshold for function */ } wmatrix; /* an array type for reals, such as f(B,L) and Ri(B,L) */ typedef double blarray[(long)t - (long)a + 1][maxwidth]; /* end module ev.type version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module evd.var */ Static allfile all; /* the data read in */ Static _TEXT evdp; /* parameters */ Static _TEXT display; /* the data redisplayed */ Static _TEXT sites; /* the raw sequences of the sites and surround */ Static _TEXT genomes; /* the raw sequences of the genomes */ Static _TEXT evfeatures; /* features for the lister program */ Static everything e; /* all the data */ Static jmp_buf _JL1; /* end module evd.var */ /* 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 = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module writedatetime */ Static Void writedatetime(thefile, adatetime) _TEXT *thefile; Char *adatetime; { /* expand the date and time out and print in the file */ long index; /* index of datetime */ for (index = 0; index < datetimearraylength; index++) putc(adatetime[index], thefile->f); } /* end module writedatetime version = 'cdatemod.p 1.19 1999Dec13'; */ /* begin module basetochar */ 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; } /* end module basetochar version = 3.97; (@ of ev.p 2005 Jul 20 */ 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; } /* ************************************************************************ */ /* ************************************************************************ */ /* begin module ev.decode */ Static long decode(e, bug, z, basesperinteger, precalcnegative, precalctwonegative) everything *e; long bug, z, basesperinteger, precalcnegative, precalctwonegative; { /* decode the genetic sequence into an integer. Starting at position z, the ; basesperinteger bases are translated into a two's complement integer. precalcnegative and precalctwonegative are precalculated values for forming the two's complement. see procedure makeprecalc. */ long val = 0; /* intermediate values of the decoding */ long place; /* place value of a value */ long x; /* location on the genome */ genometype *WITH; /* convert one stretch of the genome into an integer. the width of the stretch is basesperinteger. */ x = z + basesperinteger; WITH = &e->p.creature[bug-1]; for (place = 1; place <= basesperinteger; place++) val += e->precalc.power[place-1] [P_getbits_UB(WITH->genome, x - place - 1, 1, 3) - (long)a]; /* convert to two's complement notation */ if (val >= precalcnegative) val -= precalctwonegative; return val; } /* end module ev.decode version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module ev.translate */ Static Void translate(e, bug, z, w) everything *e; long bug, z; wmatrix *w; { /* convert the gene of the bug at location z into its recognizer w matrix */ long x = z; /* location on the genome */ long l; /* location on the recognizer */ base b; /* base touched by a finger of recognizer at l */ genometype *WITH; long FORLIM; /* writeln(list,'* bug = ',bug:1,' '); */ /* base 1 of the gene is at position z */ WITH = &e->p.creature[bug-1]; FORLIM = e->p.width; for (l = 0; l < FORLIM; l++) { for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) { w->matrix[(long)b - (long)a] [l] = decode(e, bug, x, e->p.bpi, e->precalc.negative, e->precalc.twonegative); x += e->p.bpi; /*;writeln(list,'* x=',x:2, ' w(',ord(b):1,',',l:2,') = ',val:4); */ } } /* decode the threshold from sequence just following the rest of the gene */ /* WARNING: if bpi <> bpthreshold then negative and twonegative are not correct below. See makeprecalc. */ w->threshold = decode(e, bug, x, e->p.bpthreshold, e->precalc.negative, e->precalc.twonegative); /* writeln(output,'threshold for bug ',bug:2,'=',w.threshold:3); */ } /* end module ev.translate version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module ev.score */ Static double score(w, width, creature, x, e) wmatrix *w; long width; uchar *creature; long x; everything *e; { /* Evaluate the w matrix placed on the genome at position x. Note: x+1 is the first base touched by the w matrix. That is, x is the zero of the w, and the first base of the w is at x+1. Note that the score function does not imply how the recognition occurs. That is done by the recognize function. */ double realval = 0.0; /* the value of the site evaluated by w */ long integerval; /* the value of the site evaluated by w */ long l; /* index to w matrix */ if (!e->noise) { integerval = 0; for (l = 0; l < width; l++) integerval += w->matrix[P_getbits_UB(creature, l + x, 1, 3) - (long)a] [l]; return integerval; } e->noise = false; /* */ } /* end module ev.score version = 3.97; (@ of ev.p 2005 Jul 20 */ /* the following segment of code, when placed into the score function, will print out the values used, and their sum ;for l:=1 to width do write(output,' ',w[creature[l+x],l]:3); ;writeln(output,' = ',val:5); */ /* begin module ev.recognize */ Static boolean recognize(w, width, creature, x, e) wmatrix *w; long width; uchar *creature; long x; everything *e; { /* Evaluate the w matrix placed on the genome at position x. Note: x+1 is the first base touched by the w matrix. That is, x is the zero of the w, and the first base of the w is at x+1. The test for recognition is that score value is greater than the threshold. */ /* */ return (score(w, width, creature, x, e) >= w->threshold); /* */ } /* end module ev.recognize version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module ev.evaluate */ Static Void evaluate(e, bug) everything *e; long bug; { /* evaluate the particular bug and put the number of mistakes into the rank array. */ wmatrix w; /* the recognizer; translation of the gene */ long m = 0; /* the current number of mistakes counted */ long x; /* position on the genome */ boolean recognized; /* a site was recognized at x */ genometype *WITH; long FORLIM; /*writeln(output,'in evaluate, bug: ',bug:1);*/ /* translate the recognizer gene into the w matrix */ translate(e, bug, 1L, &w); /* scan the genome */ WITH = &e->p.creature[bug-1]; FORLIM = e->precalc.endscang; for (x = 0; x <= FORLIM; x++) { recognized = recognize(&w, e->p.width, WITH->genome, x, e); if (P_getbits_UB(e->p.sitelocations, x - 1, 0, 3) != recognized) m++; } /* record the results */ e->p.rank[bug-1][(long)bugid] = bug; /* identify the bug */ e->p.rank[bug-1][(long)mistakes] = m; /*;writeln(output, 'bug ',bug:1, ' makes ',m:1,' mistakes');*/ } /* Local variables for listparams: */ struct LOC_listparams { _TEXT *list; } ; Local Void truthof(x, LINK) boolean x; struct LOC_listparams *LINK; { /* print the truth value of x to file list*/ if (x) fprintf(LINK->list->f, "* true "); else fprintf(LINK->list->f, "* false"); } /* end module ev.evaluate version = 3.97; (@ of ev.p 2005 Jul 20 */ /* ************************************************************************ */ /* begin module ev.listparams */ Static Void listparams(list_, e) _TEXT *list_; everything *e; { /* create the header of the file and list the parameters */ struct LOC_listparams V; population *WITH; V.list = list_; WITH = &e->p; fprintf(V.list->f, "*\n"); fprintf(V.list->f, "* %*ld number of creatures\n", col, WITH->bugs); fprintf(V.list->f, "* %*ld number of potential binding sites, Genome\n", col, WITH->Genome); fprintf(V.list->f, "* %*ld number of sites per creature, gamma\n", col, WITH->gamma); fprintf(V.list->f, "* %*ld width of the recognizer in bases\n", col, WITH->width); fprintf(V.list->f, "* %*ld bases per integer of the recognizer\n", col, WITH->bpi); fprintf(V.list->f, "* %*ld mutation rate in hits per creature per generation\n", col, e->mu); fprintf(V.list->f, "* %*.12f initial seed for the random number generator (given in evp)\n", col, e->initialseed); fprintf(V.list->f, "* %*ld cycles\n", col, e->cycles); fprintf(V.list->f, "* %*ld display interval\n", col, e->displayinterval); truthof(e->displayaverage, &V); fprintf(V.list->f, " : listing average mistakes\n"); truthof(e->displaychange, &V); fprintf(V.list->f, " : listing changes of mistakes\n"); truthof(e->displayindividual, &V); fprintf(V.list->f, " : listing individuals\n"); truthof(e->displaystatus, &V); fprintf(V.list->f, " : listing status to output\n"); truthof(e->displayorbits, &V); fprintf(V.list->f, " : listing Rindividual orbits\n"); truthof(e->displayrsequence, &V); fprintf(V.list->f, " : listing Rsequence\n"); truthof(e->displaygenomich, &V); fprintf(V.list->f, " : listing genomic uncertainty\n"); truthof(e->displayonemistake, &V); fprintf(V.list->f, " : listing mistakes of best individual\n"); fprintf(V.list->f, "* %*.12f current seed for the random number generator\n", col, e->seed); truthof(e->selecting, &V); fprintf(V.list->f, " : selecting by mistakes\n"); fprintf(V.list->f, "* %*ld storage frequency\n", col, e->storagefrequency); fprintf(V.list->f, "* %*.*f sigma, standard deviation of noise\n", col, dec, e->sigma); if (e->noise) fprintf(V.list->f, "* %*c Noise is modeled\n", col, ' '); else fprintf(V.list->f, "* %*c No noise is modeled\n", col, ' '); fprintf(V.list->f, "* %c SPECIAL RULE for ties: r: random, b: both, s: sort\n", e->specialrule); } /* end module ev.listparams version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module ev.geteverything */ Static Void geteverything(e, all) everything *e; allfile *all; { /* get everything in the all file. Note that this does not do a reset */ fread(e, sizeof(everything), 1, all->f); } /* end module ev.geteverything version = 3.97; (@ of ev.p 2005 Jul 20 */ /* begin module evd.getevdparams */ Static Void getevdparams(evdp, first, last, nonsites, e) _TEXT *evdp; long *first, *last; boolean *nonsites; everything *e; { /* get the parameters: first and last creatures to display and whether to display nonsites to evfeatures */ double parameterversion; /* parameter version number */ if (*evdp->name != '\0') { if (evdp->f != NULL) evdp->f = freopen(evdp->name, "r", evdp->f); else evdp->f = fopen(evdp->name, "r"); } else rewind(evdp->f); if (evdp->f == NULL) _EscIO2(FileNotFound, evdp->name); RESETBUF(evdp->f, Char); if (BUFEOF(evdp->f)) { *first = 1; *last = 1; *nonsites = false; return; } fscanf(evdp->f, "%lg%*[^\n]", ¶meterversion); getc(evdp->f); if ((long)floor(100 * parameterversion + 0.5) < (long)floor(100 * updateversion + 0.5)) { printf("You have an old parameter file!\n"); printf("parameterversion is %4.2f\n", parameterversion); printf("updateversion is %4.2f\n", updateversion); halt(); } fscanf(evdp->f, "%ld%*[^\n]", first); getc(evdp->f); fscanf(evdp->f, "%ld%*[^\n]", last); getc(evdp->f); if (BUFEOF(evdp->f)) *nonsites = false; else { if (P_peek(evdp->f) == 'n') *nonsites = true; else *nonsites = false; fscanf(evdp->f, "%*[^\n]"); getc(evdp->f); } if (*first > *last) { printf(" first creature (=%ld) must be less than or equal to last (=%ld)\n", *first, *last); halt(); } if (*first < 1) { printf(" first parameter must be positive\n"); halt(); } if (*last > e->p.bugs) { printf(" warning: there are only %ld creatures\n", e->p.bugs); *last = e->p.bugs; } } /* end module evd.getevdparams */ /* begin module evd.showandtell */ Static Void showandtell(display, e) _TEXT *display; everything e; { /* show and tell about a bunch of interesting facts */ long area; /* area in bits devoted to the recognizer gene */ long accuracy; /* maximum number of bits of accuracy that the recognizer gene could handle */ /* tell some facts about general things */ fprintf(display->f, " from version %4.2f of ev\n\n", e.precalc.evversion); fprintf(display->f, " date-time stamp: "); writedatetime(display, e.datetimestamp); fprintf(display->f, "\n\n %*ld generation\n\n", col, e.generation); /* tell some facts about the population */ fprintf(display->f, " rfrequency = %10.5f bits per site\n\n", e.precalc.rfrequency); area = e.p.width * e.p.bpi * 8; accuracy = (long) floor(log(e.p.width * exp(e.p.bpi * log(4.0))) / e.precalc.ln2 + 0.5); fprintf(display->f, " bpi = bases per integer\n"); fprintf(display->f, " recognizer gene area = (4 x width) x (2 x bpi) = %ld bits / recognizer gene\n", area); fprintf(display->f, " recognizer accuracy = log2 (width x 4^bpi) = %ld bits\n", accuracy); fprintf(display->f, " recognizer gene effectiveness = accuracy/area = %5.1f%%\n\n", accuracy * 100.0 / area); fprintf(display->f, " maximum site information = 2 x width = %5ld bits/site\n", e.p.width * 2); fprintf(display->f, " probable site density = rfrequency/(2*width) = %5.3f bits/site\n", e.precalc.rfrequency / (e.p.width * 2)); if (e.precalc.rfrequency > e.p.width * 2) { fprintf(display->f, "\n WARNING: requested site density exceeds 1.0\n"); printf(" WARNING: requested site density exceeds 1.0\n"); putc('\n', display->f); } putc('\n', display->f); } /* end module evd.showandtell */ /* begin module evd.showsites */ Static Void showsites(sites, e, first, last) _TEXT *sites; everything *e; long first, last; { /* show the sequences of the sites for bugs of rank first to last */ long bug; /* index to bugs */ long s; /* index to sites of a bug */ long x; /* index to genome of a bug */ population *WITH; genometype *WITH1; long FORLIM1, FORLIM2; if (*sites->name != '\0') { if (sites->f != NULL) sites->f = freopen(sites->name, "w", sites->f); else sites->f = fopen(sites->name, "w"); } else { if (sites->f != NULL) rewind(sites->f); else sites->f = tmpfile(); } if (sites->f == NULL) _EscIO2(FileNotFound, sites->name); SETUPBUF(sites->f, Char); WITH = &e->p; for (bug = first - 1; bug < last; bug++) { WITH1 = &WITH->creature[WITH->rank[bug][(long)bugid] - 1]; FORLIM1 = WITH->gamma; for (s = 0; s < FORLIM1; s++) { FORLIM2 = WITH->site[s] + WITH->width + extra; /* writeln(output,'zero of site[',s:1,'] is at coordinate ',site[s]:1); */ /*zzz the first position is site[s]+1!!!*/ for (x = WITH->site[s] - extra; x < FORLIM2; x++) fputc(basetochar((base)P_getbits_UB(WITH1->genome, x, 1, 3)), sites->f); fprintf(sites->f, ".\n"); } } if (first != last) /* separate sets */ putc('\n', sites->f); } /* end module evd.showsites */ /* begin module evd.showgenomes */ Static Void showgenomes(genomes, e, first, last) _TEXT *genomes; everything *e; long first, last; { /* show the sequences for bugs of rank first to last of the genomes */ long bug; /* index to bugs */ long x; /* index to genome of a bug */ population *WITH; genometype *WITH1; long FORLIM1; if (*genomes->name != '\0') { if (genomes->f != NULL) genomes->f = freopen(genomes->name, "w", genomes->f); else genomes->f = fopen(genomes->name, "w"); } else { if (genomes->f != NULL) rewind(genomes->f); else genomes->f = tmpfile(); } if (genomes->f == NULL) _EscIO2(FileNotFound, genomes->name); SETUPBUF(genomes->f, Char); WITH = &e->p; for (bug = first - 1; bug < last; bug++) { WITH1 = &WITH->creature[WITH->rank[bug][(long)bugid] - 1]; FORLIM1 = WITH->genomesize; for (x = 1; x <= FORLIM1; x++) { if (x % (deflinewidth * 2) == 0) putc('\n', genomes->f); /* p2c: evd.p, line 853: * Note: Using % for possibly-negative arguments [317] */ fputc(basetochar((base)P_getbits_UB(WITH1->genome, x - 1, 1, 3)), genomes->f); } fprintf(genomes->f, ".\n\n"); } } /* end module evd.showgenomes */ /* begin module leftjustify */ Static Void leftjustify(afile, number, width) _TEXT *afile; long number, width; { /* print the number into afile in a field of some width, but left justified */ long numberlength; /* number of characters taken up by a number */ /* put spaces on the right side */ /* add 1 to abs(number) to assure that ln can't bomb */ numberlength = (long)floor(1.5 + log(labs(number) + 1.0) / log(10.0) + 0.5); if (width < numberlength) { /* won't fit */ printf( "the number (%*ld) cannot fit in the display properly ... sorry if there is a mess\n", (int)width, number); fprintf(afile->f, "%*ld", (int)width, number); /* try it anyway... */ return; } fprintf(afile->f, "%*ld", (int)numberlength, number); if (width > numberlength) fprintf(afile->f, "%*c", (int)(width - numberlength), ' '); /* put blank padding */ } /* nrange = 0.5; (* normalization range is to between -nrange to +nrange *) */ #define nrange 1.0 /* normalization range is to between -nrange to +nrange */ #define tolerance 1e-10 /* tolerance for sum of probabilities */ /* end module leftjustify */ /* begin module evd.thermo */ Static Void thermo(display, bug, e) _TEXT *display; long bug; everything *e; { /* compute thermo and information values for the weight matrix of the bug */ base b; /* an index (with l) for the weight matrix */ double boundary = 0.0; /* the largest absolute value of the weight matrix at l */ double h; /* uncertainty computed from p */ double p; /* computed probablity - approximately a frequency */ double psum; /* sum of p to check program */ long l; /* an index (with b) for the weight matrix */ double Q; /* the partition function */ double Rsl; /* information content at position l in bits */ double Rs = 0.0; /* total information content in bits */ wmatrix w; /* the wmatrix from a bug */ population *WITH; genometype *WITH1; long FORLIM; WITH = &e->p; /*zzzqqq*/ WITH1 = &WITH->creature[WITH->rank[bug-1][(long)bugid] - 1]; /* translate this bug's matrix: */ translate(e, WITH->rank[bug-1][(long)bugid], 1L, &w); fprintf(display->f, " Thermodynamic probabilities and information of weight matrix\n"); FORLIM = WITH->width; /* scale into a range from -1 to +1 so that exp does not croak */ for (l = 0; l < FORLIM; l++) { for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) { /* determine min and max */ if (labs(w.matrix[(long)b - (long)a][l]) > boundary) boundary = labs(w.matrix[(long)b - (long)a][l]); } } fprintf(display->f, " abs(largest weight): %4ld\n", (long)floor(boundary + 0.5)); fprintf(display->f, " normalization range: %4.2f to %4.2f\n", -nrange, nrange); /* adjust the normalization range */ fprintf(display->f, " normalization value: %4ld\n", (long)floor(boundary + 0.5)); fprintf(display->f, "%4c", 'L'); for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) fprintf(display->f, " %4c", basetochar(b)); fprintf(display->f, " | "); for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) fprintf(display->f, " p(%c)", basetochar(b)); fprintf(display->f, " | "); fprintf(display->f, " bits\n"); FORLIM = WITH->width; for (l = 0; l < FORLIM; l++) { fprintf(display->f, "%4ld", l + 1); Q = 0.0; for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) { fprintf(display->f, " %4ld", w.matrix[(long)b - (long)a][l]); Q += exp(w.matrix[(long)b - (long)a][l] / boundary); } fprintf(display->f, " | "); h = 0.0; psum = 0.0; for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) { p = exp(w.matrix[(long)b - (long)a][l] / boundary) / Q; psum += p; h -= p * log(p) / log(2.0); fprintf(display->f, " %4.2f", p); } Rsl = 2 - h; Rs += Rsl; fprintf(display->f, " | %5.2f", Rsl); /* make sure that the probability sum is very close to 1.0 */ if (1.0 - psum > tolerance) fprintf(display->f, " | probabilty sum differs from 1 by % .1E", 1.0 - psum); putc('\n', display->f); } fprintf(display->f, " Thermodynamic Matrix Information = %5.2f bits\n\n", Rs); } #undef nrange #undef tolerance typedef enum { nonsite, begsite, midsite, endsite } sitetype; /* the beginning, middle and end of a site are indicated */ /* end module evd.thermo */ /* begin module evd.showcreatures */ Static Void showcreatures(display, e, first, last) _TEXT *display; everything *e; long first, last; { /* show to the display the creatures of rank first to last */ base b; /* an index (with l) that rotates through the bases for printing the wmatrix values */ long bug; /* index to the bugs */ long linewidth = deflinewidth; /* the number of bases per line */ long l; /* an index (with b) that rotates through the bases for printing the wmatrix values */ sitetype lsitestate; /* like sthisline: the sitestate at the beginning of the line. */ boolean padding; /* true means that the site finger or the site was split between lines. In this case that some padding blanks have to be put at the beginning of the next line, to keep things straight. */ Char sitecharacter; /* the character to print to mark a site: '+' means the site is recognized, '-' means it is not recognized */ sitetype sitestate; /* state of printing a site */ long s; /* which site we are at */ long sthisline; /* the site we are looking at at the start of the current line. this way s can be reset for printing the site evaluations, after being changed for printing the (--) marks */ double squared; /* the sum of squares of the weight matrix values */ double thescore; /* current score for a site by using the wmatrix */ long twobpi; /* twice the bpi, calculated once for speed */ long twowidth; /* twice the width of sites, calculated once for speed */ wmatrix w; /* the wmatrix from a bug */ long wbound; /* the upper boundary of the recognizer gene */ long x, y; /* position in a genome */ long yupper; /* the upper value to plot (never above g) */ population *WITH; genometype *WITH1; long FORLIM1, TEMP; fprintf(display->f, " the creatures that rank from %ld to %ld are displayed\n\n", first, last); WITH = &e->p; if ((linewidth & 3) != 0) { /* adjust linewidth to put sites in evenly */ linewidth -= linewidth & 3; } twobpi = WITH->bpi * 2; twowidth = WITH->width * 2; for (bug = first; bug <= last; bug++) { WITH1 = &WITH->creature[WITH->rank[bug-1][(long)bugid] - 1]; fprintf(display->f, " creature of rank %ld makes %ld mistakes\n", bug, e->p.rank[bug-1][(long)mistakes]); if (WITH->rank[bug-1][(long)bugid] == 1) fprintf(display->f, " this creature makes the fewest mistakes of all\n"); putc('\n', display->f); /* translate this bug's matrix: */ translate(e, WITH->rank[bug-1][(long)bugid], 1L, &w); fprintf(display->f, " the w matrix is:\n"); fprintf(display->f, "%4c", 'L'); for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) fprintf(display->f, " %4c", basetochar(b)); putc('\n', display->f); squared = 0.0; FORLIM1 = WITH->width; for (l = 0; l < FORLIM1; l++) { fprintf(display->f, "%4ld", l + 1); for (b = a; (long)b <= (long)t; b = (base)((long)b + 1)) { fprintf(display->f, " %4ld", w.matrix[(long)b - (long)a][l]); TEMP = w.matrix[(long)b - (long)a][l]; squared = TEMP * TEMP + squared; } putc('\n', display->f); } fprintf(display->f, " threshold = %ld\n\n", w.threshold); fprintf(display->f, " matrix magnitude = %10.5f\n", sqrt(squared)); fprintf(display->f, " sigma = %10.5f\n", e->sigma); fprintf(display->f, " width = %10ld\n", WITH->width); fprintf(display->f, " 4 * width * sigma = %10.5f\n\n", WITH->width * 4 * e->sigma); /*zzzqqq*/ thermo(display, bug, e); x = 1; wbound = WITH->width * WITH->bpi * 4; sitestate = nonsite; s = 1; b = a; l = 1; while (x < WITH->genomesize) { /* the range of the plot is x to yupper */ yupper = x + linewidth - 1; if (yupper > WITH->genomesize) yupper = WITH->genomesize; /* display sequence */ putc(' ', display->f); for (y = x - 1; y < yupper; y++) fprintf(display->f, " %c", basetochar((base)P_getbits_UB(WITH1->genome, y, 1, 3))); fprintf(display->f, " %ld\n", yupper); /* display recognizer gene and the sites */ putc(' ', display->f); sthisline = s; lsitestate = sitestate; for (y = x; y <= yupper; y++) { if (y == WITH->site[s-1] + 1) { sitestate = begsite; if (recognize(&w, WITH->width, WITH1->genome, WITH->site[s-1], e)) sitecharacter = '+'; else sitecharacter = '-'; } else if (y >= WITH->site[s-1] + WITH->width && sitestate != nonsite) { /* the > makes sure we print even if its a mess... */ sitestate = endsite; if (s < WITH->gamma) s++; } if (y <= wbound) { /* display the recognizer gene */ if ((y - 1) % WITH->bpi == 0) { /* if at start of a finger... */ switch (((y - 1) / WITH->bpi) & 3) { case 0: fprintf(display->f, "|a"); break; case 1: fprintf(display->f, " c"); break; case 2: fprintf(display->f, " g"); break; case 3: fprintf(display->f, " t"); break; } } /* p2c: evd.p, line 1104: * Note: Using % for possibly-negative arguments [317] */ else fprintf(display->f, "--"); } else if (y <= wbound + WITH->bpthreshold) { /* display threshold */ if (y == wbound + 1) fprintf(display->f, "|T"); else fprintf(display->f, "TT"); } else { switch (sitestate) { case nonsite: fprintf(display->f, " "); break; case begsite: fprintf(display->f, " ("); sitestate = midsite; break; case midsite: fprintf(display->f, "%c%c", sitecharacter, sitecharacter); break; case endsite: fprintf(display->f, "%c)", sitecharacter); sitestate = nonsite; /* display sites */ break; } } } fprintf(display->f, "\n "); /* display w matrix components and site evaluations */ s = sthisline; sitestate = lsitestate; for (y = x; y <= yupper; y++) { if (y == WITH->site[s-1] + 1) sitestate = begsite; else if (y >= WITH->site[s-1] + WITH->width && sitestate != nonsite) { /* the > makes sure we print even if its a mess... */ sitestate = endsite; if (s < WITH->gamma) s++; /* point to next site */ } if (y <= wbound) { /* display w matrix components */ if ((y - 1) % WITH->bpi == 0) { /* if at start of a finger... */ if (b == a) putc('|', display->f); else putc(' ', display->f); /* see explanation of padding below */ leftjustify(display, w.matrix[(long)b - (long)a][l-1], twobpi - 1); if (y > yupper - WITH->bpi) padding = true; else padding = false; if (b != t) b = (base)((long)b + 1); /* rotate to next base */ else { l++; /* rotate to next position */ b = a; } } /* p2c: evd.p, line 1152: * Note: Using % for possibly-negative arguments [317] */ else if (padding) fprintf(display->f, " "); if (y == wbound) padding = false; /* reset for sites */ } else if (y == wbound + 1) { /* display threshold */ putc('|', display->f); leftjustify(display, w.threshold, WITH->bpthreshold * 2 - 1); } else if (y > wbound + WITH->bpthreshold) { /* display site scores, once beyond the threshold... */ /* The number is always put at the beginning of the site. This requres the spaces to be put on the right. When a site is split between lines the missing spaces are padded into the beginning of the next line, using variable padding */ switch (sitestate) { case nonsite: fprintf(display->f, " "); break; case begsite: thescore = score(&w, WITH->width, WITH1->genome, WITH->site[s-1], e); putc(' ', display->f); leftjustify(display, (long)floor(thescore + 0.5), twowidth - 1); if (y > yupper - WITH->width) padding = true; else padding = false; sitestate = midsite; break; case midsite: /* filler stuff at the begin of a line replaces material lost at the end of previous line */ if (padding) fprintf(display->f, " "); break; case endsite: if (padding) { fprintf(display->f, " "); /* finish off */ padding = false; /* only pad the first site */ } sitestate = nonsite; break; } } } fprintf(display->f, "\n\n"); x += linewidth; } } } /* Local variables for showfeatures: */ struct LOC_showfeatures { _TEXT *features; everything *e; long ascore; /* a score of a position in the genome */ boolean atsite; /* we are at a site */ long bug; /* index to bugs */ boolean done; /* done searching for whether x is at a site */ wmatrix w; /* the wmatrix from a bug */ } ; Local Void dononsites(LINK) struct LOC_showfeatures *LINK; { /* do non-sites */ long x; /* index to genome of a bug (gpc demands that index be local) */ long s; /* index to genome of a bug (gpc demands that index be local) */ population *WITH; genometype *WITH1; long FORLIM; WITH = &LINK->e->p; WITH1 = &WITH->creature[WITH->rank[LINK->bug-1][(long)bugid] - 1]; /* ok, now add sites located that are NOT on the list */ fprintf(LINK->features->f, "* Sites not at correct places (nonsite features)\n"); FORLIM = LINK->e->precalc.endscang; for (x = 0; x <= FORLIM; x++) { /* search through the sites to see if this is one of them */ s = 1; LINK->done = false; LINK->atsite = false; while (!LINK->done) { if (WITH->site[s-1] == x) { LINK->atsite = true; LINK->done = true; } if (s == WITH->gamma) LINK->done = true; else s++; } if (!LINK->atsite) { LINK->ascore = (long)floor(score(&LINK->w, WITH->width, WITH1->genome, x, LINK->e) + 0.5); if (LINK->ascore > LINK->w.threshold) { fprintf(LINK->features->f, "@ p%ld %ld +1 \"EvSite", LINK->bug, x + 1); putc('!', LINK->features->f); fprintf(LINK->features->f, "\" \"%ld evaluation: %ld\"\n", x, LINK->ascore); } } } putc('\n', LINK->features->f); } /* dononsites */ /* end module evd.showcreatures */ /* begin module evd.showfeatures */ Static Void showfeatures(features_, e_, first, last, nonsites) _TEXT *features_; everything *e_; long first, last; boolean nonsites; { /* show the sequences of the sites for bugs of rank first to last */ struct LOC_showfeatures V; long s; /* index to sites of a bug */ long x; /* index to genome of a bug */ long wbound; /* the upper boundary of the recognizer gene on the genome */ population *WITH; genometype *WITH1; long FORLIM1; V.features = features_; V.e = e_; if (*V.features->name != '\0') { if (V.features->f != NULL) V.features->f = freopen(V.features->name, "w", V.features->f); else V.features->f = fopen(V.features->name, "w"); } else { if (V.features->f != NULL) rewind(V.features->f); else V.features->f = tmpfile(); } if (V.features->f == NULL) _EscIO2(FileNotFound, V.features->name); SETUPBUF(V.features->f, Char); WITH = &V.e->p; fprintf(V.features->f, "* evd %4.2f evdfeatures: features for the lister program\n", version); fprintf(V.features->f, "define \"EvSite+\" \"+\" \"()\" \"()\" -0.5 %1.1f\n", WITH->width - 1 + 0.5); fprintf(V.features->f, "define \"EvSite-\" \"-\" \"()\" \"()\" -0.5 %1.1f\n", WITH->width - 1 + 0.5); fprintf(V.features->f, "define \"EvSite!\" \"!\" \"()\" \"()\" -0.5 %1.1f\n", WITH->width - 1 + 0.5); fprintf(V.features->f, "define \"Threshold\" \"-\" \"--\" \"--\" -0 %ld\n", WITH->bpi - 1); fprintf(V.features->f, "define \"A\" \"-\" \"--\" \"--\" 0 %ld\n", WITH->bpi - 1); fprintf(V.features->f, "define \"C\" \"-\" \"--\" \"--\" 0 %ld\n", WITH->bpi - 1); fprintf(V.features->f, "define \"G\" \"-\" \"--\" \"--\" 0 %ld\n", WITH->bpi - 1); fprintf(V.features->f, "define \"T\" \"-\" \"--\" \"--\" 0 %ld\n", WITH->bpi - 1); for (V.bug = first; V.bug <= last; V.bug++) { WITH1 = &WITH->creature[WITH->rank[V.bug-1][(long)bugid] - 1]; /* translate this bug's matrix: */ translate(V.e, WITH->rank[V.bug-1][(long)bugid], 1L, &V.w); FORLIM1 = WITH->width; /* write weight matrix */ for (x = 0; x < FORLIM1; x++) { fprintf(V.features->f, "@ p%ld %ld +1 \"A\" \"%ld\"\n", V.bug, WITH->bpi * x * 4 + 1, V.w.matrix[0][x]); fprintf(V.features->f, "@ p%ld %ld +1 \"C\" \"%ld\"\n", V.bug, WITH->bpi * x * 4 + WITH->bpi + 1, V.w.matrix[(long)c - (long)a][x]); fprintf(V.features->f, "@ p%ld %ld +1 \"G\" \"%ld\"\n", V.bug, WITH->bpi * x * 4 + WITH->bpi * 2 + 1, V.w.matrix[(long)g - (long)a][x]); fprintf(V.features->f, "@ p%ld %ld +1 \"T\" \"%ld\"\n", V.bug, WITH->bpi * x * 4 + WITH->bpi * 3 + 1, V.w.matrix[(long)t - (long)a][x]); } /* write threshold */ wbound = WITH->width * WITH->bpi * 4; fprintf(V.features->f, "@ p%ld %ld +1 \"Threshold\" \"%ld\"\n", V.bug, wbound + 1, V.w.threshold); FORLIM1 = WITH->gamma; /* write the sites */ for (s = 1; s <= FORLIM1; s++) { V.ascore = (long)floor(score(&V.w, WITH->width, WITH1->genome, WITH->site[s-1], V.e) + 0.5); /*zzz NOTE LOCATION OF SITE IS AT +1 OF WEIGHT MATRIX!! */ fprintf(V.features->f, "@ p%ld %ld +1 \"EvSite", V.bug, WITH->site[s-1] + 1); if (V.ascore < V.w.threshold) putc('-', V.features->f); else putc('+', V.features->f); fprintf(V.features->f, "\" \"%ld evaluation: %ld\"\n", s, V.ascore); } putc('\n', V.features->f); if (nonsites) dononsites(&V); } } /* end module evd.showfeatures */ /* begin module evd.themain */ Static Void themain(all, evdp, display, evfeatures) allfile *all; _TEXT *evdp, *display, *evfeatures; { /* the main procedure of evd */ long first, last; /* what creatures to display */ boolean nonsites; /* whether to display non-sites in evfeatures */ printf(" evd %4.2f\n", version); if (*display->name != '\0') { if (display->f != NULL) display->f = freopen(display->name, "w", display->f); else display->f = fopen(display->name, "w"); } else { if (display->f != NULL) rewind(display->f); else display->f = tmpfile(); } if (display->f == NULL) _EscIO2(FileNotFound, display->name); SETUPBUF(display->f, Char); fprintf(display->f, " evd %4.2f evolution of sites display\n", version); /* obtain the all data */ if (*all->name != '\0') { if (all->f != NULL) all->f = freopen(all->name, "rb", all->f); else all->f = fopen(all->name, "rb"); } else rewind(all->f); if (all->f == NULL) _EscIO2(FileNotFound, all->name); RESETBUF(all->f, everything); if (BUFEOF(all->f)) { printf(" all file is empty\n"); halt(); } geteverything(&e, all); /* grab the data into e */ if (e.p.randomsites != 'a') { printf("***************************************************\n"); printf("* WARNING: evd does not know how to display sites *\n"); printf("* when they are not in a regular array! *\n"); printf("* The display is not reliable. Use lister. *\n"); printf("***************************************************\n"); } listparams(display, &e); putc('\n', display->f); getevdparams(evdp, &first, &last, &nonsites, &e); showandtell(display, e); showcreatures(display, &e, first, last); showsites(&sites, &e, first, last); showgenomes(&genomes, &e, first, last); showfeatures(evfeatures, &e, first, last, nonsites); } /* end module evd.themain */ main(argc, argv) int argc; Char *argv[]; { PASCAL_MAIN(argc, argv); if (setjmp(_JL1)) goto _L1; evfeatures.f = NULL; strcpy(evfeatures.name, "evfeatures"); genomes.f = NULL; strcpy(genomes.name, "genomes"); sites.f = NULL; strcpy(sites.name, "sites"); display.f = NULL; strcpy(display.name, "display"); evdp.f = NULL; strcpy(evdp.name, "evdp"); all.f = NULL; strcpy(all.name, "all"); themain(&all, &evdp, &display, &evfeatures); _L1: if (all.f != NULL) fclose(all.f); if (evdp.f != NULL) fclose(evdp.f); if (display.f != NULL) fclose(display.f); if (sites.f != NULL) fclose(sites.f); if (genomes.f != NULL) fclose(genomes.f); if (evfeatures.f != NULL) fclose(evfeatures.f); exit(EXIT_SUCCESS); } /* evd */ /* End. */