1 #define VERSION "1.50b"
2 /* version 1.50b (c) Copyright 1993 by Joseph Felsenstein.
3    Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont,
4    Andrew Keeffe.
5    Permission is granted to copy and use this program provided no fee is
6    charged for it and provided that this copyright notice is not removed.
7 */
8 
9 /* Modified by Peter Beerli and Jon Yamato for Lamarc package  */
10 
11 /* machine-specific stuff:
12    based on a number of factors in the library stdlib.h, we will try
13    to determine what kind of machine/compiler this program is being
14    built on.  However, it doesn't always succeed.  However, if you have
15    ANSI conforming C, it will probably work.
16 
17  we will try to figure out machine type
18  based on defines in stdio, and compiler-defined things as well.: */
19 
20 #ifndef FLUC_INCLUDE
21 #define FLUC_INCLUDE
22 #endif
23 
24 #ifdef  GNUDOS
25 #define DJGPP
26 #define DOS
27 #endif
28 
29 #ifdef THINK_C
30 #define MAC
31 #endif
32 
33 /*      for Metrowerks compiler
34         GENERATINGPOWERPC       - Compiler is generating PowerPC
35                                   instructions
36         GENERATING68K           - Compiler is generating 68k
37                                   family instructions
38         GENERATING68881         - Compiler is generating mc68881
39                                   floating point instructions
40 
41         SystemSevenFiveOrLater  - Compiled code will only be run on a
42                                   System 7.5 or later Macintosh
43         SystemSevenOrLater      - Compiled code will only be run
44                                   on a System 7.0 or later Macintosh
45         SystemSixOrLater        - Compiled code will only be run
46                                   on a System 6.0 or later Macintosh */
47 #ifdef GENERATINGPOWERPC
48 #define MAC
49 #endif
50 
51 #ifdef GENERATING68K
52 #define MAC
53 #endif
54 
55 #ifdef GENERATING68881
56 #define MAC
57 #endif
58 
59 #ifdef SystemSevenFiveOrLater
60 #define MAC
61 #endif
62 
63 #ifdef SystemSevenOrLater
64 #define MAC
65 #endif
66 
67 #ifdef SystemSixOrLater
68 #define MAC
69 #endif
70 
71 #include <stdio.h>
72 #include <stdlib.h>
73 
74 #ifdef __CMS_OPEN
75 #define CMS
76 #define EBCDIC true
77 #define INFILE "infile data"
78 #define OUTFILE "outfile data"
79 #define TREEFILE "treefile data"
80 #define FONTFILE "fontfile data"
81 #define PLOTFILE "plotfile data"
82 #define INTREE "intree data"
83 #define OUTTREE "outtree data"
84 #else
85 #define EBCDIC false
86 #define INFILE "infile"
87 #define OUTFILE "outfile"
88 #define TREEFILE "treefile"
89 #define FONTFILE "fontfile" /* on unix this might be /usr/local/lib/fontfile */
90 #define PLOTFILE "plotfile"
91 #define INTREE "intree"
92 #define OUTTREE "outtree"
93 #endif
94 
95 #ifdef L_ctermid            /* try and detect for sysV or V7. */
96 #define SYSTEM_FIVE
97 #endif
98 
99 #ifdef sequent
100 #define SYSTEM_FIVE
101 #endif
102 
103 #ifndef SYSTEM_FIVE
104 # include<stdlib.h>
105 # if defined(_STDLIB_H_) || defined(_H_STDLIB) || defined(H_SCCSID) || defined(unix)
106 # define UNIX
107 # define MACHINE_TYPE "BSD Unix C"
108 # endif
109 #endif
110 
111 
112 #ifdef __STDIO_LOADED
113 # ifdef __DECC /* DEC C version provided by David Mathog */
114 # define MACHINE_TYPE "OpenVMS DECC"
115 # else
116 # define VMS /* probably obsolete VAX C version */
117 # define MACHINE_TYPE "VAX/VMS C"
118 # define printf vax_printf_is_broken
119 # define fprintf vax_fprintf_is_broken
120 void vax_printf_is_broken(const char *fmt,...);
121 void vax_fprintf_is_broken(FILE *fp,const char *fmt,...);
122 void vax_tweak_fmt(char *);
123 # endif
124 #endif
125 
126 #ifdef __WATCOMC__
127 #define QUICKC
128 #define WATCOM
129 #define DOS
130 #include "graph.h"
131 #endif
132 /* watcom-c has graphics library calls that are almost identical to    *
133  * quick-c, so the "QUICKC" symbol name stays.                         */
134 
135 
136 #ifdef _QC
137 #define MACHINE_TYPE "MS-DOS / Quick C"
138 #define QUICKC
139 #include "graph.h"
140 #define DOS
141 #endif
142 
143 #ifdef _DOS_MODE
144 #define MACHINE_TYPE "MS-DOS /Microsoft C "
145 #define DOS           /* DOS is  always defined if  on a dos machine */
146 #define MSC           /* MSC is defined for microsoft C              */
147 #endif
148 
149 #ifdef __MSDOS__      /* TURBO c compiler, ONLY (no other DOS C compilers) */
150 #define DOS
151 #define TURBOC
152 #include<stdlib.h>
153 #include<graphics.h>
154 #endif
155 
156 #ifdef DJGPP          /* DJ's gnu  C/C++ port */
157 #include<graphics.h>
158 #endif
159 
160 #ifndef MACHINE_TYPE
161 #define MACHINE_TYPE "ANSI C"
162 #endif
163 
164 #ifdef DOS
165 #define MALLOCRETURN void
166 #else
167 #define MALLOCRETURN void
168 #endif
169 #ifdef VMS
170 #define signed /* signed doesn't exist in VMS */
171 #endif
172 
173 /* default screen types */
174 #ifdef DOS
175 #define IBMCRT true
176 #define ANSICRT false
177 #else
178 #ifdef MAC
179 #define IBMCRT false
180 #define ANSICRT false
181 #else
182 #define IBMCRT false
183 #define ANSICRT true
184 #endif
185 #endif
186 
187 #ifdef DJGPP
188 #undef MALLOCRETURN
189 #define MALLOCRETURN void
190 #endif
191 
192 
193 /* includes: */
194 #ifdef UNIX
195 #include<strings.h>
196 #else
197 #include<string.h>
198 #endif
199 
200 #include <math.h>
201 #include <ctype.h>
202 
203 #ifdef MAC
204 #include "mac_interface.h"
205 #endif
206 
207 #define FClose(file) if (file) fclose(file) ; file=NULL
208 #define Malloc(x) mymalloc((long)x)
209 
210 typedef void *Anyptr;
211 #define Signed     signed
212 #define Const     const
213 #define Volatile  volatile
214 #define Char        char      /* Characters (not bytes) */
215 #define Static     static     /* Private global funcs and vars */
216 #define Local      static     /* Nested functions */
217 
218 typedef unsigned char boolean;
219 
220 #define true    1
221 #define false   0
222 #define SETBITS 32
223 
224 #ifdef MAC
225 MALLOCRETURN    *mymalloc(long);
226 #else
227 MALLOCRETURN    *mymalloc();
228 #endif
229 
230 #include "nu_constants.h"
231 #include <float.h>
232 
233 typedef double          sitelike[4];
234 typedef sitelike        *ratelike;
235 typedef ratelike        *phenotype;
236 typedef char            **sequence;
237 typedef double          *contribarr;
238 typedef char            naym[NMLNGTH];
239 typedef long            longer[5];
240 
241 typedef struct _node {
242   struct _node *next, *back;
243   boolean tip;
244   long number;
245   phenotype x;
246   naym nayme;
247   boolean top;
248   double v, tyme, length;
249   long xcoord, ycoord, ymin, ymax;
250 } node;
251 
252   /* this is a list describing the tree as a series of horizontal
253      slices of the tree.  Each slice starts at a node, and extends
254      down until the next cut is encountered. */
255 typedef struct tlist {
256   struct tlist *prev, *succ;
257   node *eventnode;   /* points to a random tip for the tips */
258   double age;        /* tyme from tree top                  */
259   long numbranch;
260   node **branchlist;
261 } tlist;
262 
263 typedef struct tree {
264   node   **nodep;
265   double likelihood;
266   node *root;
267   tlist *tymelist;
268 } tree;
269 
270 typedef struct rlrec {
271   double *val;
272 } rlrec;
273 
274 typedef struct valrec {
275   double rat_xi, rat_xv, zz, z1, y1, ww1, zz1, ww2, zz2, z1zz, z1yy,
276      xiz1, xiy1_xv, ww1zz1, vv1zz1, ww2zz2, vv2zz2;
277 } valrec;
278 
279 typedef struct dnadata {
280   long numseq;          /* number of individuals */
281   long sites;           /* number of base pairs */
282   sequence seqs;        /* DNA sequences (**char) */
283   double freqa, freqc, freqg, freqt, freqr, freqy, freqar, freqcy,
284     freqgr, freqty;    /* base frequencies */
285   double xi;           /* transition pool prob */
286   double xv;           /* transversion pool prob */
287   double fracchange;   /* probability of change */
288   double ttratio;      /* transition/transversion ratio */
289 } dnadata;
290 
291 typedef struct option_struct {
292   boolean ctgry, watt, printdata, usertree, progress,
293     treeprint, interleaved, ibmpc, ansi, growthused, autocorr, same_ne,
294     same_mu, freqsfrom, interact;
295   long steps[NUM_TYPE_CHAINS], numchains[NUM_TYPE_CHAINS],
296     increm[NUM_TYPE_CHAINS], numout[NUM_TYPE_CHAINS];
297 } option_struct;
298 
299 #define FILEP(A,B) (menu ? (A) : (B))
300 #define ERRFILE FILEP(stderr,simlog)
301 #define REF_CHAIN(A) (((A) < op->numchains[0]) ? 0 : (A)+1-op->numchains[0])
302 #define TYPE_CHAIN(A) (((A) < op->numchains[0]) ? 0 : 1)
303 #define MAX(A,B) (((A) > (B)) ? (A) : (B))
304 #define EXP(A) (((A) > EXPMAX) ? printf("\n\nTOO BIG\n") : exp(A))
305 
306 void openfile(FILE **fp, char *filename, char *mode, char *application,
307    char *perm);
308 double randum(void);
309 boolean readseedfile(long *inseedptr);
310 void printtymelist(void);
311 double to_real(double tyme);
312 double to_magic(double tyme);
313 double lengthof(node *p);
314 double howlong(tlist *t);
315 node *findtop(node *p);
316 
317 void VarMalloc(node *p, boolean allokate);
318 void newnode(node **p);
319 void freenode(node *p);
320 void newtymenode(tlist **t);
321 void freetymenode(tlist *t);
322 void freetymelist(tlist *t);
323 
324 void hookup(node *p, node *q);
325 void atr(node *p);
326 
327 double zerocollis(long numl, long numother, double length);
328 double onecollis(long numl, long numother, double length);
329 double twocollis(long numother, double length);
330 
331 tlist *gettymenode(long target);
332 tlist *gettyme(node *p, node *daughter1, node *daughter2,
333    node *ans);
334 void inserttymelist(node *prime);
335 void subtymelist(node *ndonor, node *nrecip);
336 void ltov(node *p);
337 
338 void getnums(void);
339 long boolcheck(char ch);
340 boolean booleancheck(char *var, char *value);
341 boolean numbercheck(char *var, char *value);
342 void readparmfile(void);
343 void getoptions(void);
344 
345 void firstinit(void);
346 void locusinit(void);
347 void inputoptions(void);
348 void setuptree(void);
349 void freetree(void);
350 
351 boolean sitecompare(long site1, long site2);
352 void makesiteptr(void);
353 void getinput(void);
354 
355 double watterson(void);
356 
357 void orient(node *p);
358 void finishsetup(node *p);
359 void initbranchlist(void);
360 void inittable(void);
361 void initweightrat(void);
362 
363 void treeout(node *p, long s, FILE **usefile);
364 
365 void evaluate(tree *tr, boolean first);
366 void nuview(node *p);
367 void update(node *p);
368 void smooth(node *p);
369 void localsmooth(node *p);
370 
371 boolean testratio(void);
372 
373 void seekch(char c);
374 void getch(char *c);
375 void processlength(node *p);
376 void addelement(node *p, long *nextnode);
377 void treeread(void);
378 
379 void treevaluate(void);
380 void localevaluate(node *p, node *pansdaught);
381 
382 void copynode(node *source, node *target);
383 
384 void joinnode(double length, node *p, node *q);
385 void constructtree(long numtips, double branchlength);
386 
387 void updateslide(node *p, boolean wanted);
388 void rebuildbranch(void);
389 
390 boolean setlength(long numl, long numother, double tstart, double tlength,
391    node *p);
392 boolean setlength2(long numother, double tstart, double tlength, node *p,
393    node *q);
394 void updatebranch(node *oldans, node *oldp, node *prime);
395 
396 long counttymelist(tlist *first, tlist *last);
397 
398 boolean accept_tree(node *prime, node *primans, node *ans,
399    node *pansdaught, node *newbr[], tlist *tplus);
400 
401 boolean slide(void);
402 void maketree(void);
403 void finalfree(void);
404 
405 void memerror(void);
406 
407 int main(int argc, char *argv[]);
408 int eof(FILE *f);
409 int eoln(FILE *f);
410 
411 /* functions donated by getdata.c */
412 extern void setupdata(dnadata **dna, long sites, long numseq);
413 extern void getdata(tree *curtree, dnadata *dna, option_struct *op,
414   FILE *infile, FILE *outfile);
415 extern void freedata(dnadata *dna);
416 extern void makevalues(dnadata *dna, long categs, tree *curtree);
417 extern void empiricalfreqs(tree *curtree, dnadata *dna, long *weight);
418 extern void getbasefreqs(dnadata *dna, option_struct *op, double
419   locus_ttratio, FILE *outfile);
420 
421 /* functions donated by wrap.c */
422 /* Rich Miller, 4/13/97 */
423 /*lint -esym(652,calloc,realloc) */
424 #if 0
425 #define calloc(nelem,size)  E_calloc((long)(nelem), size, __FILE__, __LINE__)
426 #define realloc(ptr,size)   E_realloc(ptr, (long)(size), __FILE__, __LINE__)
427 #define E_fopen(fn,mode)    Err_fopen(fn, mode, __FILE__, __LINE__)
428 #define strtok(str,delim)   E_strtok(str,delim,__FILE__, __LINE__)
429 
430 void *E_calloc(long, size_t, char *, int);
431 void *E_realloc(void *, long, char *, int);
432 FILE *Err_fopen(char *, char *, char *, int);
433 char *E_strtok(char *newstring, const char *delimiters, char
434   *CallerFile, int CallerLine);
435 #endif
436 
437 /* function donated by fluc_modellike.c */
438 extern void fluc_estimate(long chain, boolean locusend);
439 
440 /* function donated by coal_modellike.c */
441 extern double coal_singlechain(long chain, boolean chend, boolean rend);
442 extern void coal_curveplot(void);
443