source: branches/items/GDE/TREEPUZZLE/src/ml.h

Last change on this file was 655, checked in by westram, 22 years ago

Mac OSX patches from Ben Hines
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 KB
Line 
1/*
2 * ml.h
3 *
4 *
5 * Part of TREE-PUZZLE 5.0 (June 2000)
6 *
7 * (c) 1999-2000 by Heiko A. Schmidt, Korbinian Strimmer,
8 *                  M. Vingron, and Arndt von Haeseler
9 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
10 *
11 * All parts of the source except where indicated are distributed under
12 * the GNU public licence.  See http://www.opensource.org for details.
13 */
14
15
16#ifndef _ML_
17#define _ML_
18
19/* definitions */
20
21#define MINTS 0.20  /* Ts/Tv parameter */
22#define MAXTS 30.0
23#define MINYR 0.10  /* Y/R Ts parameter */
24#define MAXYR 6.00
25#define MINFI 0.00  /* fraction invariable sites */
26#define MAXFI 0.99  /* only for input */
27#define MINGE 0.01  /* rate heterogeneity parameter */
28#define MAXGE 0.99
29#define MINCAT 4    /* discrete Gamma categories */
30#define MAXCAT 16
31
32#define RMHROOT    5.0     /* upper relative bound for height of root         */
33#define MAXARC     900.0   /* upper limit on branch length (PAM) = 6.0        */
34#define MINARC     0.001   /* lower limit on branch length (PAM) = 0.00001    */
35#define EPSILON    0.0001  /* error in branch length (PAM) = 0.000001         */
36#define HEPSILON   0.0001  /* error in node and root heights                  */
37#define MAXIT      100     /* maximum number of iterates of smoothing         */
38#define MINFDIFF   0.00002 /* lower limit on base frequency differences       */
39#define MINFREQ    0.0001  /* lower limit on base frequencies = 0.01%         */
40#define NUMQBRNCH  5       /* number of branches in a quartet                 */
41#define NUMQIBRNCH 1       /* number of internal branches in a quartet        */
42#define NUMQSPC    4       /* number of sequences in a quartet                */
43
44/* 2D minimisation */
45#define PEPS1  0.01    /* epsilon substitution process estimation */
46#define PEPS2  0.01    /* epsilon rate heterogeneity estimation  */
47
48/* quartet series */
49#define MINPERTAXUM 2
50#define MAXPERTAXUM 6
51#define TSDIFF 0.20
52#define YRDIFF 0.10
53
54/* type definitions */
55
56typedef struct node
57{
58        struct node *isop;
59        struct node *kinp;
60        int descen;
61        int number;
62        double length;
63        double lengthc;
64        double varlen;
65        double height;
66        double varheight;
67        ivector paths;
68        cvector eprob;
69        dcube partials; /* partial likelihoods */
70        char *label;  /* internal labels */
71} Node;
72
73typedef struct tree
74{
75        Node *rootp;
76        Node **ebrnchp;        /* list of pointers to external branches           */
77        Node **ibrnchp;        /* list of pointers to internal branches           */
78        double lklhd;          /* total log-likelihood                            */
79        double lklhdc;         /* total log-likelihood clock                      */
80        dmatrix condlkl;       /* likelihoods for each pattern and non-zero rate  */
81        double rssleast;
82} Tree;
83
84
85/* global variables */
86
87EXTERN Node *chep;         /* pointer to current height node                  */
88EXTERN Node *rootbr;       /* pointer to root branch                          */
89EXTERN Node **heights;     /* pointer to height nodes in unrooted tree        */
90EXTERN int Numhts;         /* number of height nodes in unrooted tree         */
91EXTERN double hroot;       /* height of root                                  */
92EXTERN double varhroot;    /* variance of height of root                      */
93EXTERN double maxhroot;    /* maximal height of root                          */
94EXTERN int locroot;        /* location of root                                */
95EXTERN int numbestroot;    /* number of best locations for root               */
96EXTERN int clockmode;      /* clocklike vs. nonclocklike computation          */
97EXTERN cmatrix Identif;    /* sequence names                                  */
98EXTERN cmatrix Seqchar;    /* ML sequence data                                */
99EXTERN cmatrix Seqpat;     /* ordered site patterns                           */
100EXTERN ivector constpat;   /* indicates constant site patterns                */
101EXTERN cvector seqchi;
102EXTERN cvector seqchj;
103EXTERN dcube partiali;
104EXTERN dcube partialj;
105EXTERN dcube ltprobr;      /* transition probabilites (for all non-zero rates */
106EXTERN dmatrix Distanmat;  /* matrix with maximum likelihood distances        */
107EXTERN dmatrix Evec;       /* Eigenvectors                                    */
108EXTERN dmatrix Ievc;       /* Inverse eigenvectors                            */
109EXTERN double TSparam;     /* Ts/Tv parameter                                 */
110EXTERN double tsmean, yrmean;
111EXTERN double YRparam;     /* Y/R Ts parameter                                */
112EXTERN double geerr;       /* estimated error of rate heterogeneity           */
113EXTERN double Geta;        /* rate heterogeneity parameter                    */
114EXTERN double fracconst;   /* fraction of constant sites                      */
115EXTERN double fracconstpat;/* fraction of constant patterns                   */
116EXTERN double Proportion;  /* for tree drawing                                */
117EXTERN double tserr;       /* estimated error of TSparam                      */
118EXTERN double yrerr;       /* estimated error of YRparam                      */
119EXTERN double fracinv;     /* fraction of invariable sites                    */
120EXTERN double fierr;       /* estimated error of fracinv                      */
121EXTERN dvector Brnlength;
122EXTERN dvector Distanvec;
123EXTERN dvector Eval;       /* Eigenvalues of 1 PAM rate matrix                */
124EXTERN dvector Freqtpm;    /* base frequencies                                */
125EXTERN dvector Rates;      /* rate of each of the categories                  */
126EXTERN dmatrix iexp;
127EXTERN imatrix Basecomp;   /* base composition of each taxon                  */
128EXTERN ivector usedtaxa;   /* list needed in the input treefile procedure     */
129EXTERN int numtc;          /* auxiliary variable for printing rooted tree     */
130EXTERN int qcalg_optn;     /* use quartet subsampling algorithm               */
131EXTERN int approxp_optn;   /* approximate parameter estimation                */
132EXTERN int chi2fail;       /* flag for chi2 test                              */
133EXTERN int Converg;        /* flag for ML convergence (no clock)              */
134EXTERN int Convergc;       /* flag for ML convergence (clock)                 */
135EXTERN int data_optn;      /* type of sequence input data                     */
136EXTERN int Dayhf_optn;     /* Dayhoff model                                   */
137EXTERN int HKY_optn;       /* use HKY model                                   */
138EXTERN int Jtt_optn;       /* JTT model                                       */
139EXTERN int blosum62_optn;  /* BLOSUM 62 model                                 */
140EXTERN int mtrev_optn;     /* mtREV model                                     */
141EXTERN int cprev_optn;     /* cpREV model                                     */
142EXTERN int vtmv_optn;      /* VT model                                        */
143EXTERN int wag_optn;       /* WAG model                                       */
144EXTERN int Maxsite;        /* number of ML characters per taxum               */
145EXTERN int Maxspc;         /* number of sequences                             */
146EXTERN int mlmode;         /* quartet ML or user defined tree ML              */
147EXTERN int nuc_optn;       /* nucleotide (4x4) models                         */
148EXTERN int Numbrnch;       /* number of branches of current tree              */
149EXTERN int numcats;        /* number of rate categories                       */
150EXTERN int Numconst;       /* number of constant sites                        */
151EXTERN int Numconstpat;    /* number of constant patterns                     */
152EXTERN int Numibrnch;      /* number of internal branches of current tree     */
153EXTERN int Numitc;         /* number of ML iterations assumning clock         */
154EXTERN int Numit;          /* number of ML iterations if there is convergence */
155EXTERN int Numptrn;        /* number of site patterns                         */
156EXTERN int Numspc;         /* number of sequences of current tree             */
157EXTERN int optim_optn;     /* optimize model parameters                       */
158EXTERN int grate_optim;    /* optimize Gamma rate heterogeneity parameter     */
159EXTERN int SH_optn;        /* SH nucleotide (16x16) model                     */
160EXTERN int TN_optn;        /* use TN model                                    */
161EXTERN int tpmradix;       /* number of different states                      */
162EXTERN int fracinv_optim;  /* optimize fraction of invariable sites           */
163EXTERN int typ_optn;       /* type of PUZZLE analysis                         */
164EXTERN ivector Weight;     /* weight of each site pattern                     */
165EXTERN Tree *Ctree;        /* pointer to current tree                         */
166EXTERN ulivector badtaxon; /* involment of each taxon in a bad quartet        */
167EXTERN int qca, qcb, qcc, qcd; /* quartet currently optimized                 */
168EXTERN ivector Alias;      /* link site -> corresponding site pattern         */
169EXTERN ivector bestrate;   /* optimal assignment of rates to sequence sites   */
170
171EXTERN int bestratefound;
172
173/* function prototypes of all ml function */
174
175void convfreq(dvector);
176void a_radixsort(cmatrix, ivector, int, int, int *);
177void condenceseq(cmatrix, ivector, cmatrix, ivector, int, int, int);
178void countconstantsites(cmatrix, ivector, int, int, int *, int*);
179void evaluateseqs(void);
180void elmhes(dmatrix, ivector, int);
181void eltran(dmatrix, dmatrix, ivector, int);
182void mcdiv(double, double, double, double, double *, double *);
183void hqr2(int, int, int, dmatrix, dmatrix, dvector, dvector);
184void onepamratematrix(dmatrix);
185void eigensystem(dvector, dmatrix);
186void luinverse(dmatrix, dmatrix, int);
187void checkevector(dmatrix, dmatrix, int);
188void tranprobmat(void);
189void tprobmtrx(double, dmatrix);
190double comptotloglkl(dmatrix);
191void allsitelkl(dmatrix, dvector);
192double pairlkl(double);
193double mldistance(int, int);
194void initdistan(void);
195void computedistan(void);
196void productpartials(Node *);
197void partialsinternal(Node *);
198void partialsexternal(Node *);
199void initpartials(Tree *);
200double intlkl(double);
201void optinternalbranch(Node *);
202double extlkl(double);
203void optexternalbranch(Node *);
204void finishlkl(Node *);
205double optlkl(Tree *);
206double treelkl(Tree *);
207void luequation(dmatrix, dvector, int);
208void lslength(Tree *, dvector, int, int, dvector);
209
210void getusertree(FILE *, cvector, int);
211Node *internalnode(Tree *, char **, int *);
212void constructtree(Tree *, cvector);
213void removebasalbif(cvector);
214void makeusertree(FILE *);
215Tree *new_tree(int, int, cmatrix);
216Tree *new_quartet(int, cmatrix);
217void free_tree(Tree *, int);
218void make_quartet(int, int, int, int);
219void changedistan(dmatrix, dvector, int);
220double quartet_lklhd(int, int, int, int);
221double quartet_alklhd(int, int, int, int);
222void readusertree(FILE *);
223double usertree_lklhd(void);
224double usertree_alklhd(void);
225void mlstart(void);
226void distupdate(int, int, int, int);
227void mlfinish(void);
228void prbranch(Node *, int, int, int, ivector, ivector, FILE *);
229void getproportion(double *, dvector, int);
230void prtopology(FILE *);
231void fputphylogeny(FILE *);
232void resulttree(FILE *);
233void njtree(FILE *);
234void njdistantree(Tree *);
235void findbestratecombination(void);
236void printbestratecombination(FILE *);
237int checkedge(int);
238void fputsubstree(FILE *, Node *);
239void fputrooted(FILE *, int);
240void findheights(Node *);
241void initclock(int);
242double clock_alklhd(int);
243double heightlkl(double);
244void optheight(void);
245double rheightlkl(double);
246void optrheight(void);
247double clock_lklhd(int);
248int findrootedge(void);
249void resultheights(FILE *);
250
251double homogentest(int);
252void YangDiscreteGamma(double, int, double *);
253void updaterates(void);
254void computestat(double *, int, double *, double *);
255double quartetml(int, int, int, int);
256double opttsq(double);
257double optyrq(double);
258void optimseqevolparamsq(void);
259double opttst(double);
260double optyrt(double);
261void optimseqevolparamst(void);
262double optfi(double);
263double optge(double);
264void optimrateparams(void);
265
266int gettpmradix(void);
267void rtfdata(dmatrix, double *);
268int code2int(cvector);
269char *int2code(int);
270
271void jttdata(dmatrix, double *);
272void dyhfdata(dmatrix, double *);
273void mtrevdata(dmatrix, double *);
274void cprev45data(dmatrix, double *);
275void blosum62data(dmatrix, double *);
276void vtmvdata(dmatrix, double *);
277void wagdata(dmatrix, double *);
278
279#endif
Note: See TracBrowser for help on using the repository browser.