src/paml.h

   1 /* paml.h \r
   2 */\r
   3 \r
   4 #if (!defined PAML_H)\r
   5 #define PAML_H\r
   6 \r
   7 \r
   8 #include <stdio.h>\r
   9 #include <stdlib.h>\r
  10 #include <ctype.h>\r
  11 #include <string.h>\r
  12 #include <math.h>\r
  13 #include <search.h>\r
  14 #include <limits.h>\r
  15 #include <float.h>\r
  16 #include <time.h>\r
  17 \r
  18 #define square(a) ((a)*(a))\r
  19 #define FOR(i,n) for(i=0; i<n; i++)\r
  20 #define FPN(file) fputc('\n', file)\r
  21 #define F0 stdout\r
  22 #define min2(a,b) ((a)<(b)?(a):(b))\r
  23 #define max2(a,b) ((a)>(b)?(a):(b))\r
  24 #define swap2(a,b,y) { y=a; a=b; b=y; }\r
  25 #define Pi  3.1415926535897932384626433832795\r
  26 \r
  27 #define beep putchar('\a')\r
  28 #define spaceming2(n) ((n)*((n)*2+9+2)*sizeof(double))\r
  29 \r
  30 int ReadSeq (FILE *fout, FILE *fseq, int cleandata, int locus);\r
  31 int ScanFastaFile(FILE *f, int *ns, int *ls, int *aligned);\r
  32 void EncodeSeqs (void);\r
  33 void SetMapAmbiguity(void);\r
  34 void ReadPatternFreq (FILE* fout, char* fpattf);\r
  35 int Initialize (FILE *fout);\r
  36 int MoveCodonSeq (int ns, int ls, char *z[]);\r
  37 int PatternWeight (void);\r
  38 int PatternWeightJC69like (FILE *fout);\r
  39 int PatternWeightSimple (int CollapsJC);\r
  40 int Site2Pattern (FILE *fout);\r
  41 \r
  42 double getRoot(double (*f)(double), double (*df)(double), double initVal);\r
  43 int f_and_x(double x[], double f[], int n, int fromx, int LastItem);\r
  44 void bigexp(double lnx, double *a, double *b);\r
  45 void SetSeed (int seed, int PrintSeed);\r
  46 double rndu (void);\r
  47 void rndu_vector (double r[], int n);\r
  48 void randorder(int order[], int n, int space[]);\r
  49 #define rnduab(a,b) ((a)+((b)-(a))*rndu())\r
  50 double reflect(double x, double a, double b);\r
  51 #define rndexp(mean) (-(mean)*log(rndu()))\r
  52 double rnduM0V1 (void);\r
  53 double rndNormal(void);\r
  54 double rndBox(void);\r
  55 double rndAirplane(void);\r
  56 double rndParabola(void);\r
  57 double rndBactrian(void);\r
  58 double rndBactrianTriangle(void);\r
  59 double rndBactrianLaplace(void);\r
  60 double rndTriangle(void);\r
  61 double rndLaplace (void);\r
  62 double rndCauchy (void);\r
  63 double rndt2 (void);\r
  64 double rndt4 (void);\r
  65 double rndlogt2 (double loc, double s);\r
  66 double rndlogistic (void);\r
  67 double rndloglogistic (double loc, double s);\r
  68 double rndgamma(double alpha);\r
  69 double rndbeta(double p, double q);\r
  70 int rndpoisson(double m);\r
  71 int rndNegBinomial(double shape, double mean);\r
  72 int SampleCat (double P[], int n, double space[]);\r
  73 int MultiNomialAliasSetTable (int ncat, double prob[], double F[], int L[], double space[]);\r
  74 int MultiNomialAlias (int n, int ncat, double F[], int L[], int nobs[]);\r
  75 int MultiNomial2 (int n, int ncat, double prob[], int nobs[], double space[]);\r
  76 int DiscreteBeta (double freq[], double x[], double p, double q, int K, int UseMedian);\r
  77 int DiscreteGamma (double freqK[], double rK[], double alpha, double beta, int K, int UseMedian);\r
  78 int AutodGamma (double Mmat[], double freqK[], double rK[], double *rho1, double alfa, double rho, int K);\r
  79 \r
  80 double QuantileChi2 (double prob, double v);\r
  81 #define QuantileGamma(prob,alpha,beta) QuantileChi2(prob,2.0*(alpha))/(2.0*(beta))\r
  82 double  PDFGamma(double x, double alpha, double beta);\r
  83 #define CDFGamma(x,alpha,beta) IncompleteGamma((beta)*(x),alpha,LnGamma(alpha))\r
  84 double  PDF_InverseGamma(double x, double alpha, double beta);\r
  85 #define CDF_InverseGamma(x,alpha,beta) (1-CDFGamma(1/(x),alpha,beta))\r
  86 #define CDFChi2(x,v) CDFGamma(x,(v)/2.0,0.5)\r
  87 double PDFBeta(double x, double p, double q);\r
  88 double CDFBeta(double x, double p, double q, double lnbeta);\r
  89 double QuantileBeta(double prob, double p, double q, double lnbeta);\r
  90 double Quantile(double(*cdf)(double x,double par[]), double p, double x, double par[], double xb[2]);\r
  91 double QuantileNormal (double prob);\r
  92 double PDFNormal(double x, double mu, double sigma2);\r
  93 double logPDFNormal(double x, double mu, double sigma2);\r
  94 double CDFNormal(double x);\r
  95 double logCDFNormal(double x);\r
  96 double PDFCauchy(double x, double m, double sigma);\r
  97 double PDFloglogistic(double x, double loc, double s);\r
  98 double PDFlogt2(double x, double loc, double s);\r
  99 double PDFt2(double x, double m, double s);\r
 100 double PDFt4(double x, double m, double s);\r
 101 double PDFt(double x, double loc, double scale, double df, double lnConst);\r
 102 double CDFt(double x, double loc, double scale, double df, double lnbeta);\r
 103 double PDFSkewT(double x, double loc, double scale, double shape, double df);\r
 104 double PDFSkewN(double x, double loc, double scale, double shape);\r
 105 double logPDFSkewN(double x, double loc, double scale, double shape);\r
 106 \r
 107 int StirlingS2(int n, int k);\r
 108 double lnStirlingS2(int n, int k);\r
 109 double LnGamma(double alpha);\r
 110 #define LnBeta(p,q) (LnGamma(p) + LnGamma(q) - LnGamma(p+q))\r
 111 double DFGamma(double x, double alpha, double beta);\r
 112 double IncompleteGamma (double x, double alpha, double ln_gamma_alpha);\r
 113 #define CDFBinormal(h,k,r)  LBinormal(-(h),-(k),r)   /* CDF for bivariate normal */\r
 114 double LBinormal (double h, double k, double r);\r
 115 double logLBinormal (double h, double k, double r);\r
 116 double probBinomial (int n, int k, double p);\r
 117 double probBetaBinomial (int n, int k, double p, double q);\r
 118 long factorial(int n);\r
 119 double Binomial(double n, int k, double *scale);\r
 120 \r
 121 int GaussLegendreRule(double **x, double **w, int order);\r
 122 int GaussLaguerreRule(double **x, double **w, int order);\r
 123 double NIntegrateGaussLegendre(double(*fun)(double x), double a, double b, int order);\r
 124 \r
 125 int ScatterPlot (int n, int nseries, int yLorR[], double x[], double y[],\r
 126     int nrow, int ncol, int ForE);\r
 127 void rainbowRGB (double temperature, int *R, int *G, int *B);\r
 128 void GetIndexTernary(int *ix, int *iy, double *x, double *y, int itriangle, int K);\r
 129 \r
 130 int CodeChara (char b, int seqtype);\r
 131 int dnamaker (char z[], int ls, double pi[]);\r
 132 int picksite (char z[], int ls, int begin, int gap, char buffer[]);\r
 133 int transform (char z[], int ls, int direction, int seqtype);\r
 134 int RemoveIndel(void);\r
 135 int f_mono_di (FILE *fout, char z[], int ls, int iring, double fb1[], double fb2[], double CondP[]);\r
 136 int PickExtreme (FILE *fout, char z[], int ls, int iring, int lfrag, int ffrag[]);\r
 137 \r
 138 int print1seq (FILE*fout, char *z, int ls, int pose[]);\r
 139 void printSeqs(FILE *fout, int *pose, char keep[], int format);\r
 140 int printPatterns(FILE *fout);\r
 141 void printSeqsMgenes (void);\r
 142 int printsma (FILE*fout, char*spname[], unsigned char*z[], int ns, int l, int lline, int gap, int seqtype, \r
 143     int transformed, int simple, int pose[]);\r
 144 int printsmaCodon (FILE *fout, unsigned char * z[], int ns, int ls, int lline, int simple);\r
 145 int zztox ( int n31, int l, char z1[], char z2[], double *x );\r
 146 int testXMat (double x[]);\r
 147 double SeqDivergence (double x[], int model, double alpha, double *kapa);\r
 148 int symtest (double freq[], int n, int nobs, double space[], double *chisym, \r
 149     double* chihom);\r
 150 int dSdNNG1986(char *z1, char *z2, int lc, int icode, int transfed, double *dS, double *dN, double *Ssites, double *Nsites);\r
 151 int difcodonNG(char codon1[], char codon2[], double *SynSite,double *AsynSite, \r
 152     double *SynDif, double *AsynDif, int transfed, int icode);\r
 153 int difcodonLWL85 (char codon1[], char codon2[], double sites[3], double sdiff[3], \r
 154                    double vdiff[3], int transfed, int icode);\r
 155 int testTransP (double P[], int n);\r
 156 double testDetailedBalance (double P[], double pi[], int n);\r
 157 void pijJC69 (double pij[2], double t);\r
 158 int PMatK80 (double P[], double t, double kapa);\r
 159 int PMatT92 (double P[], double t, double kappa, double pGC);\r
 160 int PMatTN93 (double P[], double a1t, double a2t, double bt, double pi[]);\r
 161 int PMatUVRoot (double P[],double t,int n,double U[],double V[],double Root[]);\r
 162 int PMatCijk (double PMat[], double t);\r
 163 int PMatQRev(double P[], double pi[], double t, int n, double space[]);\r
 164 int EvolveHKY85 (char source[], char target[], int ls, double t, \r
 165     double rates[], double pi[], double kapa, int isHKY85);\r
 166 double DistanceIJ (int is, int js, int model, double alpha, double *kappa);\r
 167 int DistanceMatNuc (FILE *fout, FILE*f2base, int model, double alpha);\r
 168 int EigenQREVbase (FILE* fout, double kappa[], double pi[], \r
 169                    int *nR, double Root[], double Cijk[]);\r
 170 int  DistanceMatNG86 (FILE *fout, FILE*fds, FILE*fdn, FILE*dt, double alpha);\r
 171 int  setmark_61_64 (void);\r
 172 \r
 173 int BootstrapSeq (char* seqfilename);\r
 174 int rell(FILE*flnf, FILE*fout, int ntree);\r
 175 int print1site (FILE*fout, int h);\r
 176 int MultipleGenes (FILE* fout, FILE*fpair[], double space[]);\r
 177 int lfunRates (FILE* fout, double x[], int np);\r
 178 int AncestralSeqs (FILE *fout, double x[]);\r
 179 void ListAncestSeq(FILE *fout, char *zanc);\r
 180 int ChangesSites(FILE*fout, int coding, char *zanc);\r
 181 \r
 182 int NucListall(char b, int *nb, int ib[4]);\r
 183 char *getcodon(char codon[], int icodon);\r
 184 char *getAAstr(char *AAstr, int iaa);\r
 185 int Codon2AA(char codon[3], char aa[3], int icode, int *iaa);\r
 186 int DNA2protein(char dna[], char protein[], int lc, int icode);\r
 187 int printcu (FILE *f1, double fcodon[], int icode);\r
 188 int printcums (FILE *fout, int ns, double fcodons[], int code);\r
 189 int QtoPi (double Q[], double pi[], int n, double *space);\r
 190 int PtoPi (double P[], double pi[], int n, double *space);\r
 191 int PtoX (double P1[], double P2[], double pi[], double X[]);\r
 192 \r
 193 void starttimer(void);\r
 194 char* printtime(char timestr[]);\r
 195 void sleep2(int wait);\r
 196 char *strc (int n, int c);\r
 197 int printdouble(FILE*fout, double a);\r
 198 void strcase (char *str, int direction);\r
 199 void error2(char * message);\r
 200 int  indexing(double x[], int n, int index[], int descending, int space[]);\r
 201 int binarysearch (const void *key, const void *base, size_t n, size_t size, int(*compare)(const void *, const void *), int *found);\r
 202 FILE *gfopen(char *filename, char *mode);\r
 203 int  appendfile(FILE*fout, char*filename);\r
 204 \r
 205 int zero (double x[], int n);\r
 206 double sum (double x[], int n);\r
 207 int fillxc (double x[], double c, int n);\r
 208 int xtoy (double x[], double y[], int n);\r
 209 int abyx (double a, double x[], int n);\r
 210 int axtoy(double a, double x[], double y[], int n);\r
 211 int axbytoz(double a, double x[], double b, double y[], double z[], int n);\r
 212 int identity (double x[], int n);\r
 213 double distance (double x[], double y[], int n);\r
 214 double innerp(double x[], double y[], int n);\r
 215 double norm(double x[], int n);\r
 216 \r
 217 double Det3x3 (double x[3*3]);\r
 218 int matby (double a[], double b[], double c[], int n,int m,int k);\r
 219 int matbytransposed (double a[], double b_transposed[], double c[], int n, int m, int k);\r
 220 int matIout (FILE *fout, int x[], int n, int m);\r
 221 int matout (FILE *file, double x[], int n, int m);\r
 222 int matout2 (FILE *fout, double x[], int n, int m, int wid, int deci);\r
 223 int mattransp1 (double x[], int n);\r
 224 int mattransp2 (double x[], double y[], int n, int m);\r
 225 int matinv (double x[], int n, int m, double space[]);\r
 226 int matexp (double A[], int n, int nTaylorTerms, int nSquares, double space[]);\r
 227 int matsqrt (double A[], int n, double work[]);\r
 228 int CholeskyDecomp (double A[], int n, double L[]);\r
 229 int eigenQREV (double Q[], double pi[], int n, \r
 230                double Root[], double U[], double V[], double spacesqrtpi[]);\r
 231 int eigenRealSym(double A[], int n, double Root[], double work[]);\r
 232 \r
 233 int MeanVar (double x[], int n, double *mean, double *var);\r
 234 int variance (double x[], int n, int nx, double mx[], double vx[]);\r
 235 int correl (double x[], double y[], int n, double *mx, double *my,\r
 236             double *vxx, double *vxy, double *vyy, double *r);\r
 237 int comparefloat  (const void *a, const void *b);\r
 238 int comparedouble (const void *a, const void *b);\r
 239 double Eff_IntegratedCorrelationTime(double x[], int n, double *mx, double *varx);\r
 240 int HPDinterval(double x[], int n, double HPD[2], double alpha);\r
 241 int DescriptiveStatistics(FILE *fout, char infile[], int nbin, int propternary, int SkipColumns);\r
 242 int DescriptiveStatisticsSimple (FILE *fout, char infile[], int SkipColumns);\r
 243 int splitline (char line[], int fields[]);\r
 244 int scanfile (FILE*fin, int *nrecords, int *nx, int *HasHeader, char line[], int ifields[]);\r
 245 \r
 246 double bound (int nx, double x0[], double p[], double x[],\r
 247     int (*testx) (double x[], int nx));\r
 248 int gradient (int n, double x[], double f0, double g[], \r
 249     double (* fun)(double x[],int n), double space[], int Central);\r
 250 int Hessian (int nx, double x[], double f, double g[], double H[],\r
 251     double (*fun)(double x[], int n), double space[]);\r
 252 int HessianSKT2004 (double xmle[], double lnLm, double g[], double H[]);\r
 253 \r
 254 int H_end (double x0[], double x1[], double f0, double f1, double e1, double e2, int n);\r
 255 double LineSearch(double(*fun)(double x),double *f,double *x0,double xb[2],double step,double e);\r
 256 double LineSearch2 (double(*fun)(double x[],int n), double *f, double x0[], \r
 257     double p[], double h, double limit, double e, double space[], int n);\r
 258 \r
 259 void xtoFreq(double x[], double freq[], int n);\r
 260 \r
 261 \r
 262 int SetxBound (int np, double xb[][2]);\r
 263 int ming1 (FILE *fout, double *f, double (* fun)(double x[], int n),\r
 264     int (*dfun) (double x[], double *f, double dx[], int n),\r
 265     int (*testx) (double x[], int n),\r
 266     double x0[], double space[], double e, int n);\r
 267 int ming2 (FILE *fout, double *f, double (*fun)(double x[], int n),\r
 268     int (*dfun)(double x[], double *f, double dx[], int n),\r
 269     double x[], double xb[][2], double space[], double e, int n);\r
 270 int minB (FILE*fout, double *lnL,double x[],double xb[][2],double e, double space[]);\r
 271 int minB2 (FILE*fout, double *lnL,double x[],double xb[][2],double e, double space[]);\r
 272 \r
 273 \r
 274 int Newton (FILE *fout, double *f, double (* fun)(double x[], int n),\r
 275     int (* ddfun) (double x[], double *fx, double dx[], double ddx[], int n),\r
 276     int (*testx) (double x[], int n),\r
 277     double x0[], double space[], double e, int n);\r
 278 \r
 279 int nls2 (FILE *fout, double *sx, double * x0, int nx,\r
 280     int (* fun)(double x[], double y[], int nx, int ny),\r
 281     int (* jacobi)(double x[], double J[], int nx, int ny),\r
 282     int (*testx) (double x[], int nx),\r
 283     int ny, double e);\r
 284 \r
 285 int ResetFinetuneSteps(FILE *fout, double Pjump[], double finetune[], int nsteps);\r
 286 \r
 287 /* tree structure functions in treesub.c */\r
 288 void NodeToBranch (void);\r
 289 void BranchToNode (void);\r
 290 void ClearNode (int inode);\r
 291 int ReadTreeN (FILE *ftree, int *haslength, int *haslabel, int copyname, int popline);\r
 292 int ReadTreeB (FILE *ftree, int popline);\r
 293 int OutTreeN (FILE *fout, int spnames, int printopt);\r
 294 int OutTreeB (FILE *fout);\r
 295 int DeRoot (void);\r
 296 int GetSubTreeN (int keep[], int space[]);\r
 297 void printtree (int timebranches);\r
 298 void PointconPnodes (void);\r
 299 int SetBranch (double x[]);\r
 300 int DistanceMat (FILE *fout, int ischeme, double alfa, double *kapa);\r
 301 int LSDistance (double * ss, double x[], int (*testx)(double x[],int np));\r
 302 \r
 303 int StepwiseAdditionMP (double space[]);\r
 304 double MPScoreStepwiseAddition (int is, double space[], int save);\r
 305 int AddSpecies (int species, int ib);\r
 306 int StepwiseAddition (FILE* fout, double space[]);\r
 307 int readx(double x[], int *fromfile);\r
 308 double TreeScore(double x[], double space[]);\r
 309 \r
 310 int PopEmptyLines (FILE* fseq, int lline, char line[]);\r
 311 int hasbase (char *str);\r
 312 int blankline (char *str);\r
 313 \r
 314 void BranchLengthBD(int rooted, double birth, double death, double sample, \r
 315      double mut);\r
 316 int RandomLHistory (int rooted, double space[]);\r
 317 \r
 318 void Tree2Partition (char partition[]);\r
 319 int Partition2Tree (char splits[], int lsplit, int ns, int nsplit, double label[]);\r
 320 void CladeSupport (FILE *fout, char treef[], int getSnames, char mastertreef[], int pick1tree);\r
 321 int GetNSfromTreeFile(FILE *ftree, int *ns, int *ntree);\r
 322 int NSameBranch (char partition1[],char partition2[], int nib1,int nib2, int IBsame[]);\r
 323 \r
 324 int QTN93 (int model, double Q[], double kappa1, double kappa2, double pi[]);\r
 325 int RootTN93 (int ischeme, double kapa1, double kapa2, double pi[], double *scalefactor, double Root[]);\r
 326 int EigenTN93 (int ischeme, double kapa1, double kapa2, double pi[], int *nR, double Root[], double Cijk[]);\r
 327 \r
 328 int DownStatesOneNode (int ison, int father);\r
 329 int DownStates (int inode);\r
 330 int PathwayMP (FILE *fout, double space[]);\r
 331 double MPScore (double space[]);\r
 332 double RemoveMPNinfSites (double *nsiteNinf);\r
 333 int MarkStopCodons(void);\r
 334 \r
 335 int MPInformSites (void);\r
 336 int CollapsNode (int inode, double x[]);\r
 337 \r
 338 int MakeTreeIb (int ns, int Ib[], int rooted);\r
 339 int GetTreeI (int itree, int ns, int rooted);\r
 340 int NumberTrees (int ns, int rooted);\r
 341 int CountLHistories(void);\r
 342 int ListTrees (FILE* fout, int ns, int rooted);\r
 343 int GetIofTree (int rooted, int keeptree, double space[]);\r
 344 void ReRootTree (int newroot);\r
 345 int NeighborNNI (int i_tree);\r
 346 int GetLHistoryI (int iLH);\r
 347 int GetIofLHistory (void);\r
 348 int CountLHistory(char LHistories[], double space[]);\r
 349 int ReorderNodes (char LHistory[]);\r
 350 \r
 351 int GetSubSeqs(int nsnew);\r
 352 \r
 353 /* functions for evolving sequences */\r
 354 int GenerateSeq (void);\r
 355 int Rates4Sites (double rates[],double alfa,int ncatG,int ls, int cdf,\r
 356     double space[]);\r
 357 void Evolve (int inode);\r
 358 void EvolveJC (int inode);\r
 359 \r
 360 \r
 361 int ReadTreeSeqs(FILE*fout);\r
 362 int UseLocus (int locus, int copyconP, int setmodel, int setSeqName);\r
 363 int GetGtree(int locus);\r
 364 int printGtree(int printBlength);\r
 365 \r
 366 void copySptree(void);\r
 367 void printSptree(void);\r
 368 \r
 369 \r
 370 enum {BASEseq=0, CODONseq, AAseq, CODON2AAseq, BINARYseq, BASE5seq} SeqTypes;\r
 371 \r
 372 enum {PrBranch=1, PrNodeNum=2, PrLabel=4, PrAge=8, PrOmega=16} OutTreeOptions;\r
 373 \r
 374 \r
 375 /* use mean (0; default) for discrete gamma instead of median (1) */\r
 376 #define DGammaUseMedian 0\r
 377 \r
 378 \r
 379 #define FAST_RANDOM_NUMBER\r
 380 \r
 381 #define mBactrian  0.95\r
 382 #define sBactrian  sqrt(1-mBactrian*mBactrian)\r
 383 #define aBox 0.5\r
 384 #define bBox (sqrt(12 - 3*aBox*aBox) - aBox) / 2\r
 385 #define aAirplane 1.0\r
 386 #define aParab 1.0\r
 387 \r
 388 \r
 389 #define MAXNFIELDS 320000\r
 390 \r
 391 #define PAML_RELEASE      0\r
 392 \r
 393 #define FullSeqNames      0   /* 1: numbers at the beginning of sequence name are part of name */\r
 394 \r
 395 #define pamlVerStr "paml version 4.8a, July 2014"\r
 396 \r
 397 #endif\r