ml.h 15 KB
Newer Older
Andreas Tille's avatar
Andreas Tille committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373
/*
 * ml.h
 *
 *
 * Part of TREE-PUZZLE 5.2 (July 2004)
 *
 * (c) 2003-2004 by Heiko A. Schmidt, Korbinian Strimmer, and Arndt von Haeseler
 * (c) 1999-2003 by Heiko A. Schmidt, Korbinian Strimmer,
 *                  M. Vingron, and Arndt von Haeseler
 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
 *
 * All parts of the source except where indicated are distributed under
 * the GNU public licence.  See http://www.opensource.org for details.
 *
 * ($Id$)
 *
 */


#ifndef _ML_
#define _ML_

/* values for mlmode */
#define ML_QUART 1
#define ML_TREE  2
#define NJ_TREE  3

#if 0
/* values for rhettype */
#define MLUNIFORMRATE 0
#define MLGAMMARATE   1
#define MLTWORATE     2
#define MLMIXEDRATE   3
#endif

/* values for rhetmode */
#define UNIFORMRATE 0
#define GAMMARATE   1
#define TWORATE     2
#define MIXEDRATE   3


/* definitions */

#define MINTS 0.20  /* Ts/Tv parameter */
#define MAXTS 30.0
#define MINYR 0.10  /* Y/R Ts parameter */
#define MAXYR 6.00
#define MINFI 0.00  /* fraction invariable sites */
#define MAXFI 0.99  /* only for input */
#define MINGE 0.01  /* rate heterogeneity parameter */
#define MAXGE 0.99
#define MINCAT 4    /* discrete Gamma categories */
#define MAXCAT 16

#define RMHROOT         5.0     /* upper relative bound for height of root   */    
#define MAXARC          900.0   /* upper limit on branch length (PAM) = 6.0  */
#define MINARC          0.001   /* lower limit on branch length (PAM) = 0.00001 */

#ifdef USE_ADJUSTABLE_EPS
    /* error in branch length (PAM) = 0.000001   */
#   define EPSILON_BRANCH_DEFAULT  0.0001  
#   define EPSILON_BRANCH epsilon_branch
    EXTERN double epsilon_branch;
#else
#   define EPSILON_BRANCH  0.0001  
#endif

#ifdef USE_ADJUSTABLE_EPS
/* error in node and root heights            */
#   define EPSILON_HEIGHTS_DEFAULT 0.0001
#   define EPSILON_HEIGHTS epsilon_heights
    EXTERN double epsilon_heights;
#else
#   define EPSILON_HEIGHTS 0.0001
#endif

#ifdef USE_ADJUSTABLE_EPS
#   if (EXTERN != extern)
	double epsilon_branch = EPSILON_BRANCH_DEFAULT;
	double epsilon_heights = EPSILON_HEIGHTS_DEFAULT;
#   endif
#endif

#define MAXIT           100     /* maximum number of iterates of smoothing   */
#define MINFDIFF        0.00002 /* lower limit on base frequency differences */
#define MINFREQ         0.0001  /* lower limit on base frequencies = 0.01%   */
#define NUMQBRNCH       5       /* number of branches in a quartet           */
#define NUMQIBRNCH      1       /* number of internal branches in a quartet  */
#define NUMQSPC         4       /* number of sequences in a quartet          */

/* 2D minimisation */

#ifdef USE_ADJUSTABLE_EPS
    /* epsilon substitution process estimation */
#   define EPSILON_SUBSTPARAM_DEFAULT 0.01
#   define EPSILON_SUBSTPARAM epsilon_substparam 
    EXTERN double epsilon_substparam;
#else
#   define EPSILON_SUBSTPARAM 0.01
#endif

#ifdef USE_ADJUSTABLE_EPS
    /* epsilon rate heterogeneity estimation  */
#   define EPSILON_RATEPARAM_DEFAULT  0.01
#   define EPSILON_RATEPARAM  epsilon_rateparam
    EXTERN double epsilon_rateparam;
#else
#   define EPSILON_RATEPARAM 0.01
#endif

#ifdef USE_ADJUSTABLE_EPS
#   if (EXTERN != extern)
	double epsilon_substparam = EPSILON_SUBSTPARAM_DEFAULT;
	double epsilon_rateparam = EPSILON_RATEPARAM_DEFAULT;
#   endif
#endif

/* quartet series */
#define MINPERTAXUM 2
#define MAXPERTAXUM 6
#define TSDIFF 0.20
#define YRDIFF 0.10

/* type definitions */

typedef struct node
{
	struct node *isop;	/* pointers within node (nodes) */
	struct node *kinp;	/* pointers between nodes (branches) */
	int descen;
	int number;
	double length;		/* branch length inferred (no clock) */
	double lengthc;		/* branch length inferred (clock) */
	double lengthext;	/* external length from usertree */
	int    lengthset;	/* length set by usertree */
	double varlen;
	double height;
	double varheight;
	ivector paths;		/* path directions for Korbi's p-step */
	cvector eprob;
	dcube partials;		/* partial likelihoods */
	char *label;		/* internal labels */
} Node;

typedef struct tree
{
	Node    *rootp;
	Node   **ebrnchp;	/* list of pointers to external branches     */
	Node   **ibrnchp;       /* list of pointers to internal branches     */
	double   lklhd;	        /* total log-likelihood                      */
	double   lklhdc;        /* total log-likelihood clock                */
	dmatrix  condlkl;	/* lhs for each pattern and non-zero rate    */
	double   rssleast;
} Tree;


/* global variables */

EXTERN Node    *chep;      /* pointer to current height node                 */
EXTERN Node    *rootbr;    /* pointer to root branch                         */
EXTERN Node   **heights;   /* pointer to height nodes in unrooted tree       */
EXTERN int      Numhts;    /* number of height nodes in unrooted tree        */
EXTERN double   hroot;     /* height of root                                 */
EXTERN double   varhroot;  /* variance of height of root                     */
EXTERN double   maxhroot;  /* maximal height of root                         */
EXTERN int      locroot;   /* location of root                               */
EXTERN int      numbestroot; /* number of best locations for root            */
EXTERN int      clockmode; /* clocklike vs. nonclocklike computation         */
EXTERN cmatrix  Identif;   /* sequence names                                 */
EXTERN cmatrix  Namestr;   /* sequence names (delimited by '\0')             */
EXTERN cmatrix  Seqchar;   /* ML sequence data                               */
EXTERN cmatrix  Seqpat;    /* ordered site patterns                          */
EXTERN ivector  constpat;  /* indicates constant site patterns               */
EXTERN cvector  seqchi;
EXTERN cvector  seqchj;
EXTERN dcube    partiali;
EXTERN dcube    partialj;
EXTERN dcube    ltprobr;   /* transition probabilites (for all non-zero rates*/
EXTERN dmatrix  Distanmat;   /* matrix with maximum likelihood distances     */
EXTERN dmatrix  Evec;        /* Eigenvectors                                 */
EXTERN dmatrix  Ievc;        /* Inverse eigenvectors                         */
EXTERN double   TSparam;     /* Ts/Tv parameter                              */
EXTERN double   GTR_ACrate;  /* A <-> G mutation rate for GTR                */
EXTERN double   GTR_AGrate;  /* A <-> C mutation rate for GTR                */
EXTERN double   GTR_ATrate;  /* A <-> T mutation rate for GTR                */
EXTERN double   GTR_CGrate;  /* C <-> G mutation rate for GTR                */
EXTERN double   GTR_CTrate;  /* C <-> T mutation rate for GTR                */
EXTERN double   GTR_GTrate;  /* G <-> T mutation rate for GTR                */
EXTERN double   tsmean, yrmean;
EXTERN double   YRparam;     /* Y/R Ts parameter                             */
EXTERN double   geerr;       /* estimated error of rate heterogeneity        */
EXTERN double   Geta;        /* rate heterogeneity parameter                 */
EXTERN double   fracconst;   /* fraction of constant sites                   */
EXTERN double   fracconstpat;/* fraction of constant patterns                */
EXTERN double   Proportion;  /* for tree drawing                             */
EXTERN double   tserr;       /* estimated error of TSparam                   */
EXTERN double   yrerr;       /* estimated error of YRparam                   */
EXTERN double   fracinv;     /* fraction of invariable sites                 */
EXTERN double   fierr;       /* estimated error of fracinv                   */
EXTERN dvector  Brnlength;
EXTERN dvector  Distanvec;
EXTERN dvector  Eval;        /* Eigenvalues of 1 PAM rate matrix             */
EXTERN dvector  Freqtpm;     /* base frequencies                             */
EXTERN dvector  Rates;       /* rate of each of the categories               */
EXTERN dmatrix  iexp;
EXTERN imatrix  Basecomp;    /* base composition of each taxon               */
EXTERN ivector  usedtaxa;    /* list needed in the input treefile procedure  */
EXTERN int      numtc;       /* auxiliary variable for printing rooted tree  */
EXTERN int      qcalg_optn;  /* use quartet subsampling algorithm            */
EXTERN int      approxp_optn;/* approximate parameter estimation             */
EXTERN int      chi2fail;    /* flag for chi2 test                           */
EXTERN int      Converg;     /* flag for ML convergence (no clock)           */
EXTERN int      Convergc;    /* flag for ML convergence (clock)              */
EXTERN int      data_optn;   /* type of sequence input data                  */  
EXTERN int      Dayhf_optn;  /* Dayhoff model                                */
EXTERN int      HKY_optn;    /* use HKY model                                */
EXTERN int      Jtt_optn;    /* JTT model                                    */
EXTERN int      print_GTR_optn; /* print (restricted) GTR matrix */
EXTERN int      GTR_optn;    /* GTR model                                    */
EXTERN int      blosum62_optn; /* BLOSUM 62 model                            */
EXTERN int      mtrev_optn;  /* mtREV model                                  */
EXTERN int      cprev_optn;  /* cpREV model                                  */
EXTERN int      vtmv_optn;   /* VT model                                     */
EXTERN int      wag_optn;    /* WAG model                                    */
EXTERN int      Maxsite;     /* number of ML characters per taxum            */
EXTERN int      Maxspc;      /* number of sequences                          */
EXTERN int      mlmode;      /* quartet ML or user defined tree ML           */
EXTERN int      nuc_optn;    /* nucleotide (4x4) models                      */
EXTERN int      Numbrnch;    /* number of branches of current tree           */
EXTERN int      numcats;     /* number of rate categories                    */
EXTERN int      Numconst;    /* number of constant sites                     */
EXTERN int      Numconstpat; /* number of constant patterns                  */
EXTERN int      Numibrnch;   /* number of internal branches of current tree  */
EXTERN int      Numitc;      /* number of ML iterations assumning clock      */
EXTERN int      Numit;       /* number of ML iterations if convergence       */
EXTERN int      Numptrn;     /* number of site patterns                      */
EXTERN int      Numspc;      /* number of sequences of current tree          */
EXTERN int      optim_optn;  /* optimize model parameters                    */
EXTERN int      grate_optim; /* optimize Gamma rate heterogeneity parameter  */
EXTERN int      SH_optn;     /* SH nucleotide (16x16) model                  */
EXTERN int      TN_optn;     /* use TN model                                 */
EXTERN int      tpmradix;    /* number of different states                   */
EXTERN int      fracinv_optim;  /* optimize fraction of invariable sites     */
EXTERN int      typ_optn;    /* type of PUZZLE analysis                      */
EXTERN ivector  Weight;      /* weight of each site pattern                  */
EXTERN Tree    *Ctree;       /* pointer to current tree                      */
EXTERN int      qca, qcb, qcc, qcd; /* quartet currently optimized           */
EXTERN ivector  Alias;       /* link site -> corresponding site pattern      */
EXTERN ivector  bestrate;    /* optimal assignment of rates to sequence sites*/

EXTERN int      bestratefound;  /* best rate per site already reconstructed ? */

/* function prototypes of all ml function */

void convfreq(dvector);
void tp_radixsort(cmatrix, ivector, int, int, int *);
void condenceseq(cmatrix, ivector, cmatrix, ivector, int, int, int);
void countconstantsites(cmatrix, ivector, int, int, int *, int*);
void evaluateseqs(void);
void elmhes(dmatrix, ivector, int);
void eltran(dmatrix, dmatrix, ivector, int);
void mcdiv(double, double, double, double, double *, double *);
void hqr2(int, int, int, dmatrix, dmatrix, dvector, dvector);
void onepamratematrix(dmatrix);
void eigensystem(dvector, dmatrix);
void luinverse(dmatrix, dmatrix, int);
void checkevector(dmatrix, dmatrix, int);
void tranprobmat(void);
void tprobmtrx(double, dmatrix);
double comptotloglkl(dmatrix);
void allsitelkl(dmatrix, dvector);
void writesitelklmatrixphy(FILE *outf, int numut, int numsite, dmatrix allslkl, dmatrix allslklc, int clock);
void writesitelklvectorphy(FILE *outf, char *name, int numsite, dvector aslkl);
/* void writesitelklbin(FILE *, dvector); */ /* TODO */
double pairlkl(double);
double mldistance(int, int, double); /*epe*/
/* double mldistance(int, int); */
void initdistan(void);
void computedistan(void);
void productpartials(Node *);
void partialsinternal(Node *);
void partialsexternal(Node *);
void initpartials(Tree *);
double intlkl(double);
void optinternalbranch(Node *);
double extlkl(double);
void optexternalbranch(Node *);
void finishlkl(Node *);
double optlkl(Tree *);
double treelkl(Tree *);
void luequation(dmatrix, dvector, int);
void lslength(Tree *, dvector, int, int, dvector);
#if PARALLEL
	void lslength_par(Tree *, dvector, int, int, dvector);
#endif
void setextlength(Tree *, int, int, dvector);
void getusertree(FILE *, cvector, int, int);
Node *internalnode(Tree *, char **, int *);
void constructtree(Tree *, cvector, int);
void removebasalbif(cvector);
void makeusertree(FILE *, int);
Tree *new_tree(int, int, cmatrix);
Tree *new_quartet(int, cmatrix);
void free_tree(Tree *, int);
void make_quartet(int, int, int, int);
void changedistan(dmatrix, dvector, int);
double quartet_lklhd(int, int, int, int);
double quartet_alklhd(int, int, int, int);
void readusertree(FILE *, int);
double usertree_lklhd(int, int);
double usertree_alklhd(int, int);
void mlstart(void);
void distupdate(int, int, int, int);
void mlfinish(void);
void prbranch(Node *, int, int, int, ivector, ivector, FILE *);
void getproportion(double *, dvector, int);
void prtopology(FILE *);
void fputphylogeny(FILE *);
void resulttree(FILE *outfp, int usebranch);
void njtree(FILE *);
void njdistantree(Tree *);
void findbestratecombination(void);
void printbestratecombination(FILE *);
void printbestratecombinationtofile(FILE *, int);
int checkedge(int);
void fputsubstree(FILE *, Node *);
void fputrooted(FILE *, int);
void findheights(Node *);
void initclock(int);
double clock_alklhd(int);
double heightlkl(double);
void optheight(void);
double rheightlkl(double);
void optrheight(void);
double clock_lklhd(int);
int findrootedge(void);
void resultheights(FILE *);

/* mlparam.c */
double homogentest(int);
void YangDiscreteGamma(double, int, double *);
void updaterates(void);
void computestat(double *, int, double *, double *);
double quartetml(int, int, int, int);
double opttsq(double);
double optyrq(double);
void optimseqevolparamsquart(void);
double opttst(double);
double optyrt(double);
void optimseqevolparamstree(void);
double optfi(double);
double optge(double);
void optimrateparams(void);

void estimateparametersnotree(void); /* moved from puzzle1.c */
void estimateparameterstree(void);   /* moved from puzzle1.c */
/* mlparam.c end */

int gettpmradix(void);
void rtfdata(dmatrix, double *);
int code2int(cvector, int *, int *);
char *int2code(int);

void jttdata(dmatrix, double *);
void dyhfdata(dmatrix, double *);
void mtrevdata(dmatrix, double *);
void cprev45data(dmatrix, double *);
void blosum62data(dmatrix, double *);
void vtmvdata(dmatrix, double *);
void wagdata(dmatrix, double *);

#endif