4 #if (!defined PAML_H)
\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
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
27 #define beep putchar('\a')
\r
28 #define spaceming2(n) ((n)*((n)*2+9+2)*sizeof(double))
\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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
259 void xtoFreq(double x[], double freq[], int n);
\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
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
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
285 int ResetFinetuneSteps(FILE *fout, double Pjump[], double finetune[], int nsteps);
\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
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
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
310 int PopEmptyLines (FILE* fseq, int lline, char line[]);
\r
311 int hasbase (char *str);
\r
312 int blankline (char *str);
\r
314 void BranchLengthBD(int rooted, double birth, double death, double sample,
\r
316 int RandomLHistory (int rooted, double space[]);
\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
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
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
335 int MPInformSites (void);
\r
336 int CollapsNode (int inode, double x[]);
\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
351 int GetSubSeqs(int nsnew);
\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
357 void Evolve (int inode);
\r
358 void EvolveJC (int inode);
\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
366 void copySptree(void);
\r
367 void printSptree(void);
\r
370 enum {BASEseq=0, CODONseq, AAseq, CODON2AAseq, BINARYseq, BASE5seq} SeqTypes;
\r
372 enum {PrBranch=1, PrNodeNum=2, PrLabel=4, PrAge=8, PrOmega=16} OutTreeOptions;
\r
375 /* use mean (0; default) for discrete gamma instead of median (1) */
\r
376 #define DGammaUseMedian 0
\r
379 #define FAST_RANDOM_NUMBER
\r
381 #define mBactrian 0.95
\r
382 #define sBactrian sqrt(1-mBactrian*mBactrian)
\r
384 #define bBox (sqrt(12 - 3*aBox*aBox) - aBox) / 2
\r
385 #define aAirplane 1.0
\r
389 #define MAXNFIELDS 320000
\r
391 #define PAML_RELEASE 0
\r
393 #define FullSeqNames 0 /* 1: numbers at the beginning of sequence name are part of name */
\r
395 #define pamlVerStr "paml version 4.8a, July 2014"
\r