/* PhyML: a program that computes maximum likelihood phylogenies from DNA or AA homologous sequences. Copyright (C) Stephane Guindon. Oct 2003 onward. All parts of the source except where indicated are distributed under the GNU public licence. See http://www.opensource.org for details. */ #include #ifndef UTILITIES_H #define UTILITIES_H #include #include #include #include #include #include #include #include #include #include #include extern int n_sec1; extern int n_sec2; #define For(i,n) for(i=0; i(b)?(a):(b)) #define MIN(a,b) ((a)<(b)?(a):(b)) #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); #define READ 0 #define WRITE 1 #ifndef isnan # define isnan(x) \ (sizeof (x) == sizeof (long double) ? isnan_ld (x) \ : sizeof (x) == sizeof (double) ? isnan_d (x) \ : isnan_f (x)) static inline int isnan_f (float x) { return x != x; } static inline int isnan_d (double x) { return x != x; } static inline int isnan_ld (long double x) { return x != x; } #endif #ifndef isinf # define isinf(x) \ (sizeof (x) == sizeof (long double) ? isinf_ld (x) \ : sizeof (x) == sizeof (double) ? isinf_d (x) \ : isinf_f (x)) static inline int isinf_f (float x) { return isnan (x - x); } static inline int isinf_d (double x) { return isnan (x - x); } static inline int isinf_ld (long double x) { return isnan (x - x); } #endif #define MCMC_MOVE_RANDWALK_UNIFORM 0 #define MCMC_MOVE_LOG_RANDWALK_UNIFORM 1 #define MCMC_MOVE_RANDWALK_NORMAL 2 #define MCMC_MOVE_LOG_RANDWALK_NORMAL 3 #define MCMC_MOVE_SCALE_THORNE 4 #define MCMC_MOVE_SCALE_GAMMA 5 #define N_MAX_MOVES 50 #define N_MAX_NEX_COM 20 #define T_MAX_NEX_COM 100 #define N_MAX_NEX_PARM 50 #define T_MAX_TOKEN 200 #define N_MAX_MIXT_CLASSES 1000 #define NEXUS_COM 0 #define NEXUS_PARM 1 #define NEXUS_EQUAL 2 #define NEXUS_VALUE 3 #define NEXUS_SPACE 4 #define NNI_MOVE 0 #define SPR_MOVE 1 #define BEST_OF_NNI_AND_SPR 2 #define M_1_SQRT_2_PI 0.398942280401432677939946059934 #define M_SQRT_32 5.656854249492380195206754896838 #define PI 3.14159265358979311600 #define SQRT2PI 2.50662827463100024161 #define LOG2PI 1.83787706640934533908 #define LOG2 0.69314718055994528623 #define LOG_SQRT_2_PI 0.918938533204672741780329736406 #define NORMAL 1 #define EXACT 2 #define PHYLIP 0 #define NEXUS 1 #define YES 1 #define NO 0 #define TRUE 1 #define FALSE 0 #define ON 1 #define OFF 0 #define NT 0 /*! nucleotides */ #define AA 1 /*! amino acids */ #define GENERIC 2 /*! custom alphabet */ #define UNDEFINED -1 #define ACGT 0 /*! A,G,G,T encoding */ #define RY 1 /*! R,Y encoding */ #define INTERFACE_DATA_TYPE 0 #define INTERFACE_MULTIGENE 1 #define INTERFACE_MODEL 2 #define INTERFACE_TOPO_SEARCH 3 #define INTERFACE_BRANCH_SUPPORT 4 #ifndef INFINITY #define INFINITY HUGE #endif #define N_MAX_OPTIONS 100 #define T_MAX_FILE 500 #define T_MAX_LINE 2000000 #define T_MAX_NAME 100 #define T_MAX_ID 20 #define T_MAX_SEQ 2000000 #define T_MAX_OPTION 100 #define T_MAX_LABEL 10 #define T_MAX_STATE 5 #define N_MAX_LABEL 10 #define BLOCK_LABELS 100 #define NODE_DEG_MAX 500 #define BRENT_IT_MAX 500 #define BRENT_CGOLD 0.3819660 #define BRENT_ZEPS 1.e-10 #define MNBRAK_GOLD 1.618034 #define MNBRAK_GLIMIT 100.0 #define MNBRAK_TINY 1.e-20 #define ALPHA_MIN 0.04 #define ALPHA_MAX 100 #define BL_START 1.e-04 #define GOLDEN_R 0.61803399 #define GOLDEN_C (1.0-GOLDEN_R) #define N_MAX_INSERT 20 #define N_MAX_OTU 4000 #define UNLIKELY -1.e10 #define NJ_SEUIL 0.1 #define ROUND_MAX 100 #define DIST_MAX 2.00 #define AROUND_LK 50.0 #define PROP_STEP 1.0 #define T_MAX_ALPHABET 22 #define MDBL_MIN FLT_MIN #define MDBL_MAX FLT_MAX #define POWELL_ITMAX 200 #define LINMIN_TOL 2.0E-04 #define SCALE_POW 10 /*! Scaling factor will be 2^SCALE_POW or 2^(-SCALE_POW) [[ WARNING: SCALE_POW < 31 ]]*/ #define DEFAULT_SIZE_SPR_LIST 20 #define NEWICK 0 #define NEXUS 1 #define OUTPUT_TREE_FORMAT NEWICK #define MAX_PARS 1000000000 #define LIM_SCALE_VAL 1.E-50 /*! Scaling limit (deprecated) */ #define MIN_CLOCK_RATE 1.E-10 #define MIN_VAR_BL 1.E-8 #define MAX_VAR_BL 1.E+3 #define JC69 1 #define K80 2 #define F81 3 #define HKY85 4 #define F84 5 #define TN93 6 #define GTR 7 #define CUSTOM 8 #define WAG 11 #define DAYHOFF 12 #define JTT 13 #define BLOSUM62 14 #define MTREV 15 #define RTREV 16 #define CPREV 17 #define DCMUT 18 #define VT 19 #define MTMAM 20 #define MTART 21 #define HIVW 22 #define HIVB 23 #define CUSTOMAA 24 #define LG 25 // Amino acid ordering: // Ala Arg Asn Asp Cys Gln Glu Gly His Ile Leu Lys Met Phe Pro Ser Thr Trp Tyr Val #define COMPOUND_COR 0 #define COMPOUND_NOCOR 1 #define EXPONENTIAL 2 #define GAMMA 3 #define THORNE 4 #define GUINDON 5 #define STRICTCLOCK 6 #define NONE -1 #define ALRTSTAT 1 #define ALRTCHI2 2 #define MINALRTCHI2SH 3 #define SH 4 #define ABAYES 5 /* /\* Uncomment the lines below to switch to single precision *\/ */ /* typedef float phydbl; */ /* #define LOG logf */ /* #define POW powf */ /* #define EXP expf */ /* #define FABS fabsf */ /* #define SQRT sqrtf */ /* #define CEIL ceilf */ /* #define FLOOR floorf */ /* #define RINT rintf */ /* #define ROUND roundf */ /* #define TRUNC truncf */ /* #define COS cosf */ /* #define SIN sinf */ /* #define TAN tanf */ /* #define SMALL FLT_MIN */ /* #define BIG FLT_MAX */ /* #define SMALL_PIJ 1.E-10 */ /* #define BL_MIN 1.E-5 */ /* #define P_LK_LIM_INF 2.168404e-19 /\* 2^-62 *\/ */ /* #define P_LK_LIM_SUP 4.611686e+18 /\* 2^62 *\/ */ /* Uncomment the line below to switch to double precision */ typedef double phydbl; #define LOG log #define POW pow #define EXP exp #define FABS fabs #define SQRT sqrt #define CEIL ceil #define FLOOR floor #define RINT rint #define ROUND round #define TRUNC trunc #define COS cos #define SIN sin #define TAN tan #define SMALL DBL_MIN #define BIG DBL_MAX #define SMALL_PIJ 1.E-20 #define LOGBIG 690. #define LOGSMALL -690. #if !(defined PHYTIME || defined SERGEII) #define BL_MIN 1.E-8 #define BL_MAX 100. #else #define BL_MIN 1.E-6 #define BL_MAX 1. #endif /* #define P_LK_LIM_INF 7.888609052e-31 */ /* #define P_LK_LIM_MAX 1.267650600e+30 */ /* #define P_LK_LIM_INF 4.909093465e-91 /\* R: format(2^(-300),digits=10) *\/ */ /* #define P_LK_LIM_SUP 2.037035976e+90 /\* R: format(2^(+300),digits=10) *\/ */ #define P_LK_LIM_INF 3.054936e-151 /* 2^-500 */ #define P_LK_LIM_SUP 3.273391e+150 /* 2^500 */ #define T_MAX_XML_TAG 64 /*!********************************************************/ typedef struct __Scalar_Int { int v; struct __Scalar_Int *next; struct __Scalar_Int *prev; }scalar_int; /*!********************************************************/ typedef struct __Scalar_Dbl { phydbl v; bool onoff; struct __Scalar_Dbl *next; struct __Scalar_Dbl *prev; }scalar_dbl; /*!********************************************************/ typedef struct __Vect_Int { int *v; int len; struct __Vect_Int *next; struct __Vect_Int *prev; }vect_int; /*!********************************************************/ typedef struct __Vect_Dbl { phydbl *v; int len; struct __Vect_Dbl *next; struct __Vect_Dbl *prev; }vect_dbl; /*!********************************************************/ typedef struct __String { char *s; int len; struct __String *next; struct __String *prev; }t_string; /*!********************************************************/ typedef struct __Node { struct __Node **v; /*! table of pointers to neighbor nodes. Dimension = 3 */ struct __Node ***bip_node; /*! three lists of pointer to tip nodes. One list for each direction */ struct __Edge **b; /*! table of pointers to neighbor branches */ struct __Node *anc; /*! direct ancestor t_node (for rooted tree only) */ struct __Node *ext_node; struct __Node *match_node; struct __Align *c_seq; /*! corresponding compressed sequence */ struct __Align *c_seq_anc; /*! corresponding compressed ancestral sequence */ struct __Node *next; /*! tree->a_nodes[i]->next <=> tree->next->a_nodes[i] */ struct __Node *prev; /*! See above */ struct __Node *next_mixt; /*! Next mixture tree*/ struct __Node *prev_mixt; /*! Parent mixture tree */ struct __Calibration **calib; short int *calib_applies_to; int n_calib; int *bip_size; /*! Size of each of the three lists from bip_node */ int num; /*! t_node number */ int tax; /*! tax = 1 -> external node, else -> internal t_node */ int check_branch; /*! check_branch=1 is the corresponding branch is labelled with '*' */ char *name; /*! taxon name (if exists) */ char *ori_name; /*! taxon name (if exists) */ phydbl *score; /*! score used in BioNJ to determine the best pair of nodes to agglomerate */ phydbl *l; /*! lengths of the (three or one) branche(s) connected this t_node */ phydbl dist_to_root; /*! distance to the root t_node */ short int common; phydbl y_rank; phydbl y_rank_ori; phydbl y_rank_min; phydbl y_rank_max; int *s_ingrp; /*! does the subtree beneath belong to the ingroup */ int *s_outgrp; /*! does the subtree beneath belong to the outgroup */ int id_rank; /*! order taxa alphabetically and use id_rank to store the ranks */ int rank; int rank_max; }t_node; /*!********************************************************/ typedef struct __Edge { /*! syntax : (node) [edge] (left_1) . .(right_1) \ (left) (right) / \._____________./ / [b_fcus] \ / \ (left_2) . .(right_2) */ struct __Node *left,*rght; /*! t_node on the left/right side of the t_edge */ short int l_r,r_l,l_v1,l_v2,r_v1,r_v2; /*! these are directions (i.e., 0, 1 or 2): */ /*! l_r (left to right) -> left[b_fcus->l_r] = right */ /*! r_l (right to left) -> right[b_fcus->r_l] = left */ /*! l_v1 (left t_node to first t_node != from right) -> left[b_fcus->l_v1] = left_1 */ /*! l_v2 (left t_node to secnd t_node != from right) -> left[b_fcus->l_v2] = left_2 */ /*! r_v1 (right t_node to first t_node != from left) -> right[b_fcus->r_v1] = right_1 */ /*! r_v2 (right t_node to secnd t_node != from left) -> right[b_fcus->r_v2] = right_2 */ struct __NNI *nni; struct __Edge *next; struct __Edge *prev; struct __Edge *next_mixt; struct __Edge *prev_mixt; int num; /*! branch number */ scalar_dbl *l; /*! branch length */ scalar_dbl *best_l; /*! best branch length found so far */ scalar_dbl *l_old; /*! old branch length */ scalar_dbl *l_var; /*! variance of edge length */ scalar_dbl *l_var_old; /*! variance of edge length (previous value) */ int bip_score; /*! score of the bipartition generated by the corresponding edge bip_score = 1 iif the branch is found in both trees to be compared, bip_score = 0 otherwise. */ phydbl *p_lk_left,*p_lk_rght; /*! likelihoods of the subtree on the left and right side (for each site and each relative rate category) */ short int *p_lk_tip_r, *p_lk_tip_l; short int *div_post_pred_left; /*! posterior prediction of nucleotide/aa diversity (left-hand subtree) */ short int *div_post_pred_rght; /*! posterior prediction of nucleotide/aa diversity (rght-hand subtree) */ short int does_exist; int *patt_id_left; int *patt_id_rght; int *p_lk_loc_left; int *p_lk_loc_rght; phydbl *Pij_rr; /*! matrix of change probabilities and its first and secnd derivates */ int *pars_l,*pars_r; /*! parsimony of the subtree on the left and right sides (for each site) */ unsigned int *ui_l, *ui_r; /*! union - intersection vectors used in Fitch's parsimony algorithm */ int *p_pars_l, *p_pars_r; /*! conditional parsimony vectors */ int num_st_left; /*! number of the subtree on the left side */ int num_st_rght; /*! number of the subtree on the right side */ /*! Below are the likelihood scaling factors (used in functions `Get_All_Partial_Lk_Scale' in lk.c */ int *sum_scale_left_cat; int *sum_scale_rght_cat; int *sum_scale_left; int *sum_scale_rght; phydbl bootval; /*! bootstrap value (if exists) */ short int is_alive; /*! is_alive = 1 if this t_edge is used in a tree */ phydbl dist_btw_edges; int topo_dist_btw_edges; int has_zero_br_len; phydbl ratio_test; /*! approximate likelihood ratio test */ phydbl alrt_statistic; /*! aLRT statistic */ char **labels; /*! string of characters that labels the corresponding t_edge */ int n_labels; /*! number of labels */ int n_jumps; /*! number of jumps of substitution rates */ phydbl bin_cod_num; }t_edge; /*!********************************************************/ typedef struct __Tree{ struct __Node *n_root; /*! root t_node */ struct __Edge *e_root; /*! t_edge on which lies the root */ struct __Node **a_nodes; /*! array of nodes that defines the tree topology */ struct __Edge **a_edges; /*! array of edges */ struct __Model *mod; /*! substitution model */ struct __Calign *data; /*! sequences */ struct __Calign *anc_data; /*! ancestral sequences */ struct __Tree *next; /* set to NULL by default. Used for mixture models */ struct __Tree *prev; /* set to NULL by default. Used for mixture models */ struct __Tree *next_mixt; /* set to NULL by default. Used for mixture models */ struct __Tree *prev_mixt; /* set to NULL by default. Used for mixture models */ struct __Tree *mixt_tree; /* set to NULL by default. Used for mixture models */ struct __Option *io; /*! input/output */ struct __Matrix *mat; /*! pairwise distance matrix */ struct __Node **curr_path; /*! list of nodes that form a path in the tree */ struct __SPR **spr_list; struct __SPR *best_spr; struct __Tdraw *ps_tree; /*! structure for drawing trees in postscript format */ struct __T_Rate *rates; /*! structure for handling rates of evolution */ struct __Tmcmc *mcmc; struct __Triplet *triplet_struct; struct __Phylogeo *geo; int is_mixt_tree; int tree_num; /*! tree number. Used for mixture models */ int ps_page_number; /*! when multiple trees are printed, this variable give the current page number */ int depth_curr_path; /*! depth of the t_node path defined by curr_path */ int has_bip; /*!if has_bip=1, then the structure to compare tree topologies is allocated, has_bip=0 otherwise */ int n_otu; /*! number of taxa */ int curr_site; /*! current site of the alignment to be processed */ int curr_catg; /*! current class of the discrete gamma rate distribution */ int n_swap; /*! number of NNIs performed */ int n_pattern; /*! number of distinct site patterns */ int has_branch_lengths; /*! =1 iff input tree displays branch lengths */ int print_boot_val; /*! if print_boot_val=1, the bootstrap values are printed */ int print_alrt_val; /*! if print_boot_val=1, the bootstrap values are printed */ int both_sides; /*! both_sides=1 -> a pre-order and a post-order tree traversals are required to compute the likelihood of every subtree in the phylogeny*/ int num_curr_branch_available; /*!gives the number of the next cell in a_edges that is free to receive a pointer to a branch */ short int *t_dir; int n_improvements; int n_moves; int dp; /*! Data partition */ int s_mod_num; /*! Substitution model number */ int lock_topo; /*! = 1 any subsequent topological modification will be banished */ int write_labels; int write_br_lens; int *mutmap; /*! Mutational map */ phydbl init_lnL; phydbl best_lnL; /*! highest value of the loglikelihood found so far */ int best_pars; /*! highest value of the parsimony found so far */ phydbl c_lnL; /*! loglikelihood */ phydbl old_lnL; /*! old loglikelihood */ phydbl sum_min_sum_scale; /*! common factor of scaling factors */ phydbl *c_lnL_sorted; /*! used to compute c_lnL by adding sorted terms to minimize CPU errors */ phydbl *cur_site_lk; /*! vector of loglikelihoods at individual sites */ phydbl *old_site_lk; /*! vector of likelihoods at individual sites */ phydbl **log_site_lk_cat; /*! loglikelihood at individual sites and for each class of rate*/ phydbl *site_lk_cat; /*! loglikelihood at a single site and for each class of rate*/ phydbl unconstraint_lk; /*! unconstrained (or multinomial) likelihood */ phydbl **log_lks_aLRT; /*! used to compute several branch supports */ phydbl n_root_pos; /*! position of the root on its t_edge */ phydbl size; /*! tree size */ int *site_pars; int c_pars; int *step_mat; int size_spr_list; int perform_spr_right_away; time_t t_beg; time_t t_current; int bl_from_node_stamps; /*! == 1 -> Branch lengths are determined by t_node times */ phydbl sum_y_dist_sq; phydbl sum_y_dist; phydbl tip_order_score; int write_tax_names; int update_alias_subpatt; phydbl geo_mig_sd; /*! standard deviation of the migration step random variable */ phydbl geo_lnL; /*! log likelihood of the phylo-geography model */ int bl_ndigits; phydbl *short_l; /*! Vector of short branch length values */ int n_short_l; /*! Length of short_l */ phydbl norm_scale; short int br_len_recorded; int max_spr_depth; short int apply_lk_scaling; /*! Applying scaling of likelihoods. YES/NO */ phydbl *K;//a vector of the norm.constants for the node times prior. }t_tree; /*!********************************************************/ typedef struct __Super_Tree { struct __Tree *tree; struct __List_Tree *treelist; /*! list of trees. One tree for each data set to be processed */ struct __Calign *curr_cdata; struct __Option **optionlist; /*! list of pointers to input structures (used in supertrees) */ struct __Node ***match_st_node_in_gt; /*! match_st_in_gt_node[subdataset number][supertree t_node number] * gives the t_node in tree estimated from 'subdataset number' that corresponds * to 'supertree t_node number' in the supertree */ struct __Node *****map_st_node_in_gt; /*! mat_st_gt_node[gt_num][st_node_num][direction] gives the * t_node in gt gt_num that maps t_node st_node_num in st. */ struct __Edge ***map_st_edge_in_gt; /*! map_st_gt_br[gt_num][st_branch_num] gives the * branch in gt gt_num that maps branch st_branch_num * in st. */ struct __Edge ****map_gt_edge_in_st; /*! mat_gt_st_br[gt_num][gt_branch_num][] is the list of * branches in st that map branch gt_branch_num * in gt gt_num. */ int **size_map_gt_edge_in_st; /*! size_map_gt_st_br[gt_num][gt_branch_num] gives the * size of the list map_gt_st_br[gt_num][gt_branch_num][] */ struct __Edge ***match_st_edge_in_gt; /*! match_st_edge_in_gt[gt_num][st_branch_num] gives the * branch in gt gt_num that matches branch st_branch_num */ struct __Edge ***match_gt_edge_in_st; /*! match_gt_edge_in_st[gt_num][gt_branch_num] gives the * branch in st that matches branch gt_branch_num */ struct __Node ****closest_match; /*! closest_match[gt_num][st_node_num][dir] gives the * closest t_node in st that matches a t_node in gt gt_num */ int ***closest_dist; /*! closest_dist[gt_num][st_node_num][dir] gives the * number of edges to traverse to get to node * closest_match[gt_num][st_node_num][dir] */ int n_part; /*! number of trees */ phydbl **bl; /*! bl[gt_num][gt_branch_num] gives the length of * branch gt_branch_num */ phydbl **bl_cpy; /*! copy of bl */ phydbl **bl0; /*! bl estimated during NNI (original topo) * See Mg_NNI. */ phydbl **bl1; /*! bl estimated during NNI (topo conf 1) * See Mg_NNI. */ phydbl **bl2; /*! bl estimated during NNI (topo conf 2) * See Mg_NNI. */ int *bl_partition; /*! partition[gt_num] gives the t_edge partition number * gt_num belongs to. */ int n_bl_part; struct __Model **s_mod; /*! substitution model */ int n_s_mod; int lock_br_len; }supert_tree; /*!********************************************************/ typedef struct __List_Tree { /*! a list of trees */ struct __Tree **tree; int list_size; /*! number of trees in the list */ }t_treelist; /*!********************************************************/ typedef struct __Align { char *name; /*! sequence name */ int len; /*! sequence length */ char *state; /*! sequence itself */ int *d_state; /*! sequence itself (digits) */ short int *is_ambigu; /*! is_ambigu[site] = 1 if state[site] is an ambiguous character. 0 otherwise */ }align; /*!********************************************************/ typedef struct __Calign { struct __Align **c_seq; /*! compressed sequences */ phydbl *b_frq; /*! observed state frequencies */ short int *invar; /*! < 0 -> polymorphism observed */ int *wght; /*! # of each site in c_align */ short int *ambigu; /*! ambigu[i]=1 is one or more of the sequences at site i display an ambiguous character */ phydbl obs_pinvar; int n_otu; /*! number of taxa */ int clean_len; /*! uncrunched sequences lenghts without gaps */ int crunch_len; /*! crunched sequences lengths */ int init_len; /*! length of the uncompressed sequences */ int *sitepatt; /*! this array maps the position of the patterns in the compressed alignment to the positions in the uncompressed one */ int format; /*! 0 (default): PHYLIP. 1: NEXUS. */ }calign; /*!********************************************************/ typedef struct __Matrix { /*! mostly used in BIONJ */ phydbl **P,**Q,**dist; /*! observed proportions of transition, transverion and distances between pairs of sequences */ t_tree *tree; /*! tree... */ int *on_off; /*! on_off[i]=1 if column/line i corresponds to a t_node that has not been agglomerated yet */ int n_otu; /*! number of taxa */ char **name; /*! sequence names */ int r; /*! number of nodes that have not been agglomerated yet */ struct __Node **tip_node; /*! array of pointer to the leaves of the tree */ int curr_int; /*! used in the NJ/BIONJ algorithms */ int method; /*! if method=1->NJ method is used, BIONJ otherwise */ }matrix; /*!********************************************************/ typedef struct __RateMatrix { int n_diff_rr; /*! number of different relative substitution rates in the custom model */ vect_dbl *rr; /*! relative rate parameters of the GTR or custom model (given by rr_val[rr_num[i]]) */ vect_dbl *rr_val; /*! relative rate parameters of the GTR or custom model */ vect_int *rr_num; vect_int *n_rr_per_cat; /*! number of rate parameters in each category */ vect_dbl *qmat; vect_dbl *qmat_buff; struct __RateMatrix *next; struct __RateMatrix *prev; }t_rmat; /*!********************************************************/ typedef struct __RAS { /*! Rate across sites */ int n_catg; /*! number of categories in the discrete gamma distribution */ int invar; /*! =1 iff the substitution model takes into account invariable sites */ int gamma_median; /*! 1: use the median of each bin in the discrete gamma distribution. 0: the mean is used */ vect_dbl *gamma_r_proba; /*! probabilities of the substitution rates defined by the discrete gamma distribution */ vect_dbl *gamma_r_proba_unscaled; vect_dbl *gamma_rr; /*! substitution rates defined by the RAS distribution */ vect_dbl *gamma_rr_unscaled; /*! substitution rates defined by the RAS distribution (unscaled) */ scalar_dbl *alpha; /*! gamma shapa parameter */ int free_mixt_rates; scalar_dbl *free_rate_mr; /*! mean relative rate as given by the FreeRate model */ int parent_class_number; scalar_dbl *pinvar; /*! proportion of invariable sites */ short int init_rr; short int init_r_proba; short int normalise_rr; short int *skip_rate_cat; struct __RAS *next; struct __RAS *prev; }t_ras; /*!********************************************************/ typedef struct __EquFreq { /*! Equilibrium frequencies */ vect_dbl *pi; /*! states frequencies */ vect_dbl *pi_unscaled; /*! states frequencies (unscaled) */ struct __EquFreq *next; struct __EquFreq *prev; }t_efrq; /*!********************************************************/ typedef struct __Model { struct __Optimiz *s_opt; /*! pointer to parameters to optimize */ struct __Eigen *eigen; struct __M4 *m4mod; struct __Option *io; struct __Model *next; struct __Model *prev; struct __Model *next_mixt; struct __Model *prev_mixt; struct __RateMatrix *r_mat; struct __EquFreq *e_frq; struct __RAS *ras; t_string *aa_rate_mat_file; FILE *fp_aa_rate_mat; vect_dbl *user_b_freq; /*! user-defined nucleotide frequencies */ t_string *modelname; t_string *custom_mod_string; /*! string of characters used to define custom models of substitution */ int mod_num; /*! model number */ int update_eigen; /*! update_eigen=1-> eigen values/vectors need to be updated */ int whichmodel; int is_mixt_mod; int ns; /*! number of states (4 for ADN, 20 for AA) */ int bootstrap; /*! Number of bootstrap replicates (0 : no bootstrap analysis is launched) */ int use_m4mod; /*! Use a Makrkov modulated Markov model ? */ scalar_dbl *kappa; /*! transition/transversion rate */ scalar_dbl *lambda; /*! parameter used to define the ts/tv ratios in the F84 and TN93 models */ scalar_dbl *br_len_multiplier; vect_dbl *Pij_rr; /*! matrix of change probabilities */ scalar_dbl *mr; /*! mean rate = branch length/time interval mr = -sum(i)(vct_pi[i].mat_Q[ii]) */ short int log_l; /*! Edge lengths are actually log(Edge lengths) if log_l == YES !*/ phydbl l_min; /*! Minimum branch length !*/ phydbl l_max; /*! Maximum branch length !*/ phydbl l_var_sigma; /*! For any edge b we have b->l_var->v = l_var_sigma * (b->l->v)^2 */ phydbl l_var_min; /*! Min of variance of branch lengths (used in conjunction with gamma_mgf_bl == YES) */ phydbl l_var_max; /*! Max of variance of branch lengths (used in conjunction with gamma_mgf_bl == YES) */ int gamma_mgf_bl; /*! P = \int_0^inf exp(QL) p(L) where L=\int_0^t R(s) ds and p(L) is the gamma density. Set to NO by default !*/ int n_mixt_classes; /* Number of classes in the mixture model. */ scalar_dbl *r_mat_weight; scalar_dbl *e_frq_weight; }t_mod; /*!********************************************************/ typedef struct __Eigen{ int size; phydbl *q; /*! matrix which eigen values and vectors are computed */ phydbl *space; int *space_int; phydbl *e_val; /*! eigen values (vector), real part. */ phydbl *e_val_im; /*! eigen values (vector), imaginary part */ phydbl *r_e_vect; /*! right eigen vector (matrix), real part */ phydbl *r_e_vect_im; /*! right eigen vector (matrix), imaginary part */ phydbl *l_e_vect; /*! left eigen vector (matrix), real part */ struct __Eigen *prev; struct __Eigen *next; }eigen; /*!********************************************************/ typedef struct __Option { /*! mostly used in 'help.c' */ struct __Model *mod; /*! pointer to a substitution model */ struct __Tree *tree; /*! pointer to the current tree */ struct __Align **data; /*! pointer to the uncompressed sequences */ struct __Tree *cstr_tree; /*! pointer to a constraint tree (can be a multifurcating one) */ struct __Calign *cdata; /*! pointer to the compressed sequences */ struct __Super_Tree *st; /*! pointer to supertree */ struct __Tnexcom **nex_com_list; struct __List_Tree *treelist; /*! list of trees. */ struct __Option *next; struct __Option *prev; int interleaved; /*! interleaved or sequential sequence file format ? */ int in_tree; /*! =1 iff a user input tree is used as input */ char *in_align_file; /*! alignment file name */ FILE *fp_in_align; /*! pointer to the alignment file */ char *in_tree_file; /*! input tree file name */ FILE *fp_in_tree; /*! pointer to the input tree file */ char *in_constraint_tree_file; /*! input constraint tree file name */ FILE *fp_in_constraint_tree; /*! pointer to the input constraint tree file */ char *out_tree_file; /*! name of the tree file */ FILE *fp_out_tree; char *out_trees_file; /*! name of the tree file */ FILE *fp_out_trees; /*! pointer to the tree file containing all the trees estimated using random starting trees */ char *out_boot_tree_file; /*! name of the tree file */ FILE *fp_out_boot_tree; /*! pointer to the bootstrap tree file */ char *out_boot_stats_file; /*! name of the tree file */ FILE *fp_out_boot_stats; /*! pointer to the statistics file */ char *out_stats_file; /*! name of the statistics file */ FILE *fp_out_stats; char *out_trace_file; /*! name of the file in which the likelihood of the model is written */ FILE *fp_out_trace; char *out_lk_file; /*! name of the file in which the likelihood of the model is written */ FILE *fp_out_lk; char *out_ps_file; /*! name of the file in which tree(s) is(are) written */ FILE *fp_out_ps; char *clade_list_file; int datatype; /*! 0->DNA, 1->AA */ int print_boot_trees; /*! =1 if the bootstrapped trees are printed in output */ int out_stats_file_open_mode; /*! opening file mode for statistics file */ int out_tree_file_open_mode; /*! opening file mode for tree file */ int n_data_sets; /*! number of data sets to be analysed */ int n_trees; /*! number of trees */ int init_len; /*! sequence length */ int n_otu; /*! number of taxa */ int n_data_set_asked; /*! number of bootstrap replicates */ char *nt_or_cd; /*! nucleotide or codon data ? (not used) */ int multigene; /*! if=1 -> analyse several partitions. */ int config_multigene; int n_part; /*! number of data partitions */ int curr_gt; int ratio_test; /*! from 1 to 4 for specific branch supports, 0 of not */ int ready_to_go; int data_file_format; /*! Data format: Phylip or Nexus */ int tree_file_format; /*! Tree format: Phylip or Nexus */ int state_len; int curr_interface; int r_seed; /*! random seed */ int collapse_boot; /*! 0 -> branch length on bootstrap trees are not collapsed if too small */ int random_boot_seq_order; /*! !0 -> sequence order in bootstrapped data set is random */ int print_trace; int print_site_lnl; int m4_model; int rm_ambigu; /*! 0 is the default. 1: columns with ambiguous characters are discarded prior further analysis */ int colalias; int append_run_ID; char *run_id_string; int quiet; /*! 0 is the default. 1: no interactive question (for batch mode) */ int lk_approx; /* EXACT or NORMAL */ char **alphabet; int codpos; int mutmap; char **long_tax_names; char **short_tax_names; int size_tax_names; phydbl *z_scores; phydbl *lat; phydbl *lon; int boot_prog_every; int mem_question; int do_alias_subpatt; struct __Tmcmc *mcmc; struct __T_Rate *rates; }option; /*!********************************************************/ typedef struct __Optimiz { /*! parameters to be optimised (mostly used in 'optimiz.c') */ int print; /*! =1 -> verbose mode */ int opt_alpha; /*! =1 -> the gamma shape parameter is optimised */ int opt_kappa; /*! =1 -> the ts/tv ratio parameter is optimised */ int opt_lambda; /*! =1 -> the F84|TN93 model specific parameter is optimised */ int opt_pinvar; /*! =1 -> the proportion of invariants is optimised */ int opt_state_freq; /*! =1 -> the nucleotide frequencies are optimised */ int opt_rr; /*! =1 -> the relative rate parameters of the GTR or the customn model are optimised */ int opt_subst_param; /*! if opt_topo=0 and opt_subst_param=1 -> the numerical parameters of the model are optimised. if opt_topo=0 and opt_free_param=0 -> no parameter is optimised */ int opt_cov_delta; int opt_cov_alpha; int opt_cov_free_rates; int opt_bl; /*! =1 -> the branch lengths are optimised */ int opt_topo; /*! =1 -> the tree topology is optimised */ int topo_search; phydbl init_lk; /*! initial loglikelihood value */ int n_it_max; /*! maximum bnumber of iteration during an optimisation step */ int last_opt; /*! =1 -> the numerical parameters are optimised further while the tree topology remains fixed */ int random_input_tree; /*! boolean */ int n_rand_starts; /*! number of random starting points */ int brent_it_max; int steph_spr; int user_state_freq; int opt_five_branch; int pars_thresh; int hybrid_thresh; phydbl tree_size_mult; /*! tree size multiplier */ phydbl min_diff_lk_local; phydbl min_diff_lk_global; phydbl min_diff_lk_move; phydbl p_moves_to_examine; int fast_nni; int greedy; int general_pars; int quickdirty; int spr_pars; int spr_lnL; int max_depth_path; int min_depth_path; int deepest_path; phydbl max_delta_lnL_spr; int br_len_in_spr; int opt_free_mixt_rates; int constrained_br_len; int opt_gamma_br_len; int first_opt_free_mixt_rates; int wim_n_rgrft; int wim_n_globl; int wim_max_dist; int wim_n_optim; int wim_n_best; int wim_inside_opt; int opt_rmat_weight; int opt_efrq_weight; int skip_tree_traversal; int serial_free_rates; int curr_opt_free_rates; }t_opt; /*!********************************************************/ typedef struct __NNI{ struct __Node *left; struct __Node *rght; struct __Edge *b; phydbl score; phydbl init_l; phydbl init_lk; phydbl best_l; phydbl lk0,lk1,lk2; phydbl l0,l1,l2; struct __Node *swap_node_v1; struct __Node *swap_node_v2; struct __Node *swap_node_v3; struct __Node *swap_node_v4; int best_conf; /*! best topological configuration : ((left_1,left_2),right_1,right_2) or ((left_1,right_2),right_1,left_2) or ((left_1,right_1),right_1,left_2) */ }nni; /*!********************************************************/ typedef struct __SPR{ struct __Node *n_link; struct __Node *n_opp_to_link; struct __Edge *b_opp_to_link; struct __Edge *b_target; struct __Edge *b_init_target; struct __Node **path; phydbl init_target_l; phydbl init_target_v; phydbl l0,l1,l2; phydbl v0,v1,v2; phydbl lnL; int depth_path; int pars; int dist; struct __SPR *next; struct __SPR *prev; struct __SPR *next_mixt; struct __SPR *prev_mixt; }t_spr; /*!********************************************************/ typedef struct __Triplet{ int size; phydbl *F_bc; phydbl *F_cd; phydbl *F_bd; phydbl ****core; phydbl ***p_one_site; phydbl ***sum_p_one_site; phydbl *pi_bc; phydbl *pi_cd; phydbl *pi_bd; struct __Eigen *eigen_struct; struct __Model *mod; struct __Triplet *next; struct __Triplet *prev; }triplet; /*!********************************************************/ typedef struct __Pnode{ struct __Pnode **next; int weight; int num; }pnode; /*!********************************************************/ typedef struct __M4 { int n_h; /*! number of hidden states */ int n_o; /*! number of observable states */ int use_cov_alpha; int use_cov_free; phydbl **o_mats; /*! set of matrices of substitution rates across observable states */ phydbl *multipl; /*! vector of values that multiply each o_mats matrix */ phydbl *o_rr; /*! relative rates (symmetric) of substitution between observable states */ phydbl *h_rr; /*! relative rates (symmetric) of substitution between hidden states */ phydbl *h_mat; /*! matrix that describes the substitutions between hidden states (aka switches) */ phydbl *o_fq; /*! equilibrium frequencies for the observable states */ phydbl *h_fq; /*! equilibrium frequencies for the hidden states */ phydbl *h_fq_unscaled; /*! unscaled equilibrium frequencies for the hidden states */ phydbl *multipl_unscaled; /*! unscaled vector of values that multiply each o_mats matrix */ phydbl delta; /*! switching rate */ phydbl alpha; /*! gamma shape parameter */ }m4; /*!********************************************************/ typedef struct __Tdraw { phydbl *xcoord; /*! t_node coordinates on the x axis */ phydbl *ycoord; /*! t_node coordinates on the y axis */ phydbl *xcoord_s; /*! t_node coordinates on the x axis (scaled) */ phydbl *ycoord_s; /*! t_node coordinates on the y axis (scaled) */ int page_width; int page_height; int tree_box_width; int *cdf_mat; phydbl *cdf_mat_x; phydbl *cdf_mat_y; phydbl max_dist_to_root; }tdraw; /*!********************************************************/ typedef struct __T_Rate { phydbl lexp; /*! Parameter of the exponential distribution that governs the rate at which substitution between rate classes ocur */ phydbl alpha; phydbl less_likely; phydbl birth_rate; phydbl birth_rate_min; phydbl birth_rate_max; phydbl min_rate; phydbl max_rate; phydbl c_lnL1; phydbl c_lnL2; phydbl c_lnL_rates; /*! Prob(Br len | time stamps, model of rate evolution) */ phydbl c_lnL_times; /*! Prob(time stamps) */ phydbl c_lnL_jps; /*! Prob(# Jumps | time stamps, rates, model of rate evolution) */ phydbl clock_r; /*! Mean substitution rate, i.e., 'molecular clock' rate */ phydbl min_clock; phydbl max_clock; phydbl lbda_nu; phydbl min_dt; phydbl step_rate; phydbl true_tree_size; phydbl p_max; phydbl norm_fact; phydbl nu; /*! Parameter of the Exponential distribution for the corresponding model */ phydbl min_nu; phydbl max_nu; phydbl covdet; phydbl sum_invalid_areas; phydbl *nd_r; /*! Current rates at nodes */ phydbl *br_r; /*! Current rates along edges */ phydbl *nd_t; /*! Current t_node times */ phydbl *triplet; phydbl *true_t; /*! true t_node times (including root node) */ phydbl *true_r; /*! true t_edge rates (on rooted tree) */ phydbl *buff_t; phydbl *buff_r; phydbl *dens; /*! Probability densities of mean substitution rates at the nodes */ phydbl *ml_l; /*! ML t_edge lengths (rooted) */ phydbl *cur_l; /*! Current t_edge lengths (rooted) */ phydbl *u_ml_l; /*! ML t_edge lengths (unrooted) */ phydbl *u_cur_l; /*! Current t_edge lengths (unrooted) */ phydbl *invcov; phydbl *cov_r; phydbl *mean_r; /*! average values of br_r taken across the sampled values during the MCMC */ phydbl *mean_t; /*! average values of nd_t taken across the sampled values during the MCMC */ phydbl *_2n_vect1; phydbl *_2n_vect2; phydbl *_2n_vect3; phydbl *_2n_vect4; short int *_2n_vect5; phydbl *_2n2n_vect1; phydbl *_2n2n_vect2; phydbl *trip_cond_cov; phydbl *trip_reg_coeff; phydbl *cond_var; phydbl *reg_coeff; phydbl *t_prior; phydbl *t_prior_min; phydbl *t_prior_max; phydbl *t_floor; phydbl *t_mean; int *t_ranked; phydbl *mean_l; phydbl *cov_l; phydbl *grad_l; /* gradient */ phydbl inflate_var; phydbl *time_slice_lims; phydbl *survival_rank; phydbl *survival_dur; phydbl *calib_prob; int adjust_rates; /*! if = 1, branch rates are adjusted such that a modification of a given t_node time does not modify any branch lengths */ int use_rates; /*! if = 0, branch lengths are expressed as differences between t_node times */ int bl_from_rt; /*! if =1, branch lengths are obtained as the product of cur_r and t */ int approx; int model; /*! Model number */ int is_allocated; int met_within_gibbs; int update_mean_l; int update_cov_l; int *n_jps; int *t_jps; int n_time_slices; int *n_time_slice_spans; int *curr_slice; int *n_tips_below; short int *t_has_prior; struct __Node **lca; /*! 2-way table of common ancestral nodes for each pari of nodes */ short int *br_do_updt; phydbl *cur_gamma_prior_mean; phydbl *cur_gamma_prior_var; int model_log_rates; short int nd_t_recorded; short int br_r_recorded; int *has_survived; struct __Calibration *calib; int tot_num_cal; int *curr_nd_for_cal; phydbl c_lnL_Hastings_ratio; phydbl *t_prior_min_buff; phydbl *t_prior_max_buff; phydbl *times_partial_proba; }t_rate; /*!********************************************************/ typedef struct __Tmcmc { struct __Option *io; phydbl *tune_move; phydbl *move_weight; phydbl *acc_rate; int *acc_move; int *run_move; int *prev_acc_move; int *prev_run_move; int *num_move; int *move_type; char **move_name; int num_move_nd_r; int num_move_br_r; int num_move_nd_t; int num_move_nu; int num_move_clock_r; int num_move_tree_height; int num_move_subtree_height; int num_move_kappa; int num_move_tree_rates; int num_move_subtree_rates; int num_move_updown_nu_cr; int num_move_updown_t_cr; int num_move_updown_t_br; int num_move_ras; int num_move_cov_rates; int num_move_cov_switch; int num_move_birth_rate; int num_move_jump_calibration; int num_move_geo_lambda; int num_move_geo_sigma; int num_move_geo_tau; int num_move_geo_updown_tau_lbda; int num_move_geo_updown_lbda_sigma; int nd_t_digits; int *monitor; char *out_filename; time_t t_beg; time_t t_cur; time_t t_last_print; FILE *out_fp_stats; FILE *out_fp_trees; FILE *out_fp_means; FILE *out_fp_last; FILE *out_fp_constree; FILE *in_fp_par; int *adjust_tuning; int n_moves; int use_data; int randomize; int norm_freq; int run; int chain_len; int sample_interval; int chain_len_burnin; int print_every; int is_burnin; phydbl max_tune; phydbl min_tune; phydbl *old_param_val; phydbl *new_param_val; phydbl *sum_val; phydbl *sum_valsq; phydbl *sum_curval_nextval; phydbl *ess; phydbl *first_val; int *ess_run; int *start_ess; int is; /* Importance sampling? Yes or NO */ }t_mcmc; /*!********************************************************/ typedef struct __Tpart { int *ns; /*! number of states for each partition (e.g., 2, 4, 3) */ int *cum_ns; /*! cumulative number of states (e.g., 0, 2, 6) */ int ns_max; /*! maximum number of states */ int ns_min; /*! minimum number of states */ int n_partitions; /*! number of partitions */ struct __Eigen *eigen; }part; /*!********************************************************/ typedef struct __Tnexcom { char *name; int nparm; int nxt_token_t; int cur_token_t; struct __Tnexparm **parm; }nexcom; /*!********************************************************/ typedef struct __Tnexparm { char *name; char *value; int nxt_token_t; int cur_token_t; int (*fp)(char *, struct __Tnexparm *, struct __Option *); struct __Tnexcom *com; }nexparm; /*!********************************************************/ typedef struct __ParamInt { int val; }t_param_int; /*!********************************************************/ typedef struct __ParamDbl { phydbl val; }t_param_dbl; /*!********************************************************/ typedef struct __XML_node { struct __XML_attr *attr; // Pointer to the first element of a list of attributes int n_attr; // Number of attributes struct __XML_node *next; // Next sibling struct __XML_node *prev; // Previous sibling struct __XML_node *parent; // Parent of this node struct __XML_node *child; // Child of this node char *id; char *name; char *value; struct __Generic_Data_Structure *ds; // Pointer to a data strucuture. Can be a scalar, a vector, anything. }xml_node; /*!********************************************************/ typedef struct __Generic_Data_Structure { void *obj; struct __Generic_Data_Structure *next; }t_ds; /*!********************************************************/ typedef struct __XML_attr { char *name; char *value; struct __XML_attr *next; // Next attribute struct __XML_attr *prev; // Previous attribute }xml_attr; /*!********************************************************/ typedef struct __Calibration { phydbl *proba; // Probability of this calibration (set by the user and fixed throughout) phydbl lower; // lower bound phydbl upper; // upper bound int cur_applies_to; phydbl calib_proba; struct __Node **all_applies_to; int n_all_applies_to; struct __Calibration *next; struct __Calibration *prev; }t_cal; /*!********************************************************/ typedef struct __Phylogeo{ phydbl *cov; // Covariance of migrations (n_dim x n_dim) phydbl *r_mat; // R matrix. Gives the rates of migrations between locations. See article. phydbl *f_mat; // F matrix. See article. int *occup; // Vector giving the number of lineages that occupy each location phydbl *ldscape; // Coordinates of the locations int *loc; // Location for each lineage int *loc_beneath; // Gives the location occupied beneath each node in the tree int ldscape_sz; // Landscape size: number of locations int n_dim; // Dimension of the data (e.g., longitude + lattitude -> n_dim = 2) int update_fmat; phydbl sigma; // Dispersal parameter phydbl min_sigma; phydbl max_sigma; phydbl sigma_thresh; // beyond sigma_thresh, there is no dispersal bias. phydbl lbda; // Competition parameter phydbl min_lbda; phydbl max_lbda; phydbl c_lnL; struct __Node **sorted_nd; // Table of nodes sorted wrt their heights. phydbl tau; // overall migration rate parameter phydbl min_tau; phydbl max_tau; }t_geo; /*!********************************************************/ /*!********************************************************/ /*!********************************************************/ /*!********************************************************/ /*!********************************************************/ void Unroot_Tree(char **subtrees); void Init_Tree(t_tree *tree,int n_otu); void Init_Edge_Light(t_edge *b,int num); void Init_Node_Light(t_node *n,int num); void Make_Edge_Dirs(t_edge *b,t_node *a,t_node *d,t_tree *tree); void Init_NNI(nni *a_nni); void Init_Nexus_Format(nexcom **com); void Restrict_To_Coding_Position(align **data,option *io); void Uppercase(char *ch); void Lowercase(char *ch); calign *Compact_Data(align **data,option *io); calign *Compact_Cdata(calign *data,option *io); void Traverse_Prefix_Tree(int site,int seqnum,int *patt_num,int *n_patt,align **data,option *io,pnode *n); pnode *Create_Pnode(int size); void Get_Base_Freqs(calign *data); void Get_AA_Freqs(calign *data); void Swap_Nodes_On_Edges(t_edge *e1,t_edge *e2,int swap,t_tree *tree); void Connect_Edges_To_Nodes_Recur(t_node *a,t_node *d,t_tree *tree); void Connect_One_Edge_To_Two_Nodes(t_node *a,t_node *d,t_edge *b,t_tree *tree); void Update_Dirs(t_tree *tree); void Exit(char *message); void *mCalloc(int nb,size_t size); void *mRealloc(void *p,int nb,size_t size); int Sort_Phydbl_Decrease(const void *a,const void *b); void Qksort_Int(int *A,int *B,int ilo,int ihi); void Qksort(phydbl *A,phydbl *B,int ilo,int ihi); void Qksort_Matrix(phydbl **A,int col,int ilo,int ihi); void Order_Tree_Seq(t_tree *tree,align **data); char *Add_Taxa_To_Constraint_Tree(FILE *fp,calign *cdata); void Check_Constraint_Tree_Taxa_Names(t_tree *tree,calign *cdata); void Order_Tree_CSeq(t_tree *tree,calign *cdata); void Init_Mat(matrix *mat,calign *data); void Copy_Tax_Names_To_Tip_Labels(t_tree *tree,calign *data); void Share_Lk_Struct(t_tree *t_full,t_tree *t_empt); void Share_Spr_Struct(t_tree *t_full,t_tree *t_empt); void Share_Pars_Struct(t_tree *t_full,t_tree *t_empt); int Sort_Edges_NNI_Score(t_tree *tree,t_edge **sorted_edges,int n_elem); int Sort_Edges_Depth(t_tree *tree,t_edge **sorted_edges,int n_elem); void NNI(t_tree *tree,t_edge *b_fcus,int do_swap); void NNI_Pars(t_tree *tree,t_edge *b_fcus,int do_swap); void Swap(t_node *a,t_node *b,t_node *c,t_node *d,t_tree *tree); void Update_All_Partial_Lk(t_edge *b_fcus,t_tree *tree); void Update_SubTree_Partial_Lk(t_edge *b_fcus,t_node *a,t_node *d,t_tree *tree); void Copy_Seq_Names_To_Tip_Labels(t_tree *tree,calign *data); calign *Copy_Cseq(calign *ori,option *io); int Filexists(char *filename); matrix *K80_dist(calign *data,phydbl g_shape); matrix *JC69_Dist(calign *data,t_mod *mod); matrix *Hamming_Dist(calign *data,t_mod *mod); int Is_Invar(int patt_num,int stepsize,int datatype,calign *data); int Is_Ambigu(char *state,int datatype,int stepsize); void Check_Ambiguities(calign *data,int datatype,int stepsize); int Get_State_From_Ui(int ui,int datatype); int Assign_State(char *c,int datatype,int stepsize); char Reciproc_Assign_State(int i_state,int datatype); int Assign_State_With_Ambiguity(char *c,int datatype,int stepsize); void Clean_Tree_Connections(t_tree *tree); void Bootstrap(t_tree *tree); void Br_Len_Involving_Invar(t_tree *tree); void Br_Len_Not_Involving_Invar(t_tree *tree); void Getstring_Stdin(char *s); phydbl Num_Derivatives_One_Param(phydbl (*func)(t_tree *tree), t_tree *tree, phydbl f0, phydbl *param, int which, int n_param, phydbl stepsize, int logt, phydbl *err, int precise, int is_positive); phydbl Num_Derivatives_One_Param_Nonaligned(phydbl (*func)(t_tree *tree), t_tree *tree, phydbl f0, phydbl **param, int which, int n_param, phydbl stepsize, int logt, phydbl *err, int precise, int is_positive); int Num_Derivative_Several_Param(t_tree *tree,phydbl *param,int n_param,phydbl stepsize,int logt,phydbl(*func)(t_tree *tree),phydbl *derivatives, int is_positive); int Num_Derivative_Several_Param_Nonaligned(t_tree *tree, phydbl **param, int n_param, phydbl stepsize, int logt, phydbl (*func)(t_tree *tree), phydbl *derivatives, int is_positive); int Compare_Two_States(char *state1,char *state2,int state_size); void Copy_One_State(char *from,char *to,int state_size); void Copy_Dist(phydbl **cpy,phydbl **orig,int n); t_mod *Copy_Model(t_mod *ori); void Record_Model(t_mod *ori,t_mod *cpy); void Set_Defaults_Input(option *io); void Set_Defaults_Model(t_mod *mod); void Set_Defaults_Optimiz(t_opt *s_opt); void Test_Node_Table_Consistency(t_tree *tree); void Get_Bip(t_node *a,t_node *d,t_tree *tree); void Alloc_Bip(t_tree *tree); int Sort_Phydbl_Increase(const void *a,const void *b); int Sort_String(const void *a,const void *b); int Compare_Bip(t_tree *tree1,t_tree *tree2,int on_existing_edges_only); void Match_Tip_Numbers(t_tree *tree1,t_tree *tree2); void Test_Multiple_Data_Set_Format(option *io); int Are_Compatible(char *statea,char *stateb,int stepsize,int datatype); void Hide_Ambiguities(calign *data); void Copy_Tree(t_tree *ori,t_tree *cpy); void Prune_Subtree(t_node *a,t_node *d,t_edge **target,t_edge **residual,t_tree *tree); void Graft_Subtree(t_edge *target,t_node *link,t_edge *residual,t_tree *tree); void Reassign_Node_Nums(t_node *a,t_node *d,int *curr_ext_node,int *curr_int_node,t_tree *tree); void Reassign_Edge_Nums(t_node *a,t_node *d,int *curr_br,t_tree *tree); void Find_Mutual_Direction(t_node *n1,t_node *n2,short int *dir_n1_to_n2,short int *dir_n2_to_n1); void Update_Dir_To_Tips(t_node *a,t_node *d,t_tree *tree); void Fill_Dir_Table(t_tree *tree); int Get_Subtree_Size(t_node *a,t_node *d); void Fast_Br_Len(t_edge *b,t_tree *tree,int approx); void Init_Eigen_Struct(eigen *this); phydbl Triple_Dist(t_node *a,t_tree *tree,int approx); void Make_Symmetric(phydbl **F,int size); void Round_Down_Freq_Patt(phydbl **F,t_tree *tree); phydbl Get_Sum_Of_Cells(phydbl *F,t_tree *tree); void Divide_Cells(phydbl **F,phydbl div,t_tree *tree); void Divide_Mat_By_Vect(phydbl **F,phydbl *vect,int size); void Multiply_Mat_By_Vect(phydbl **F,phydbl *vect,int size); void Found_In_Subtree(t_node *a,t_node *d,t_node *target,int *match,t_tree *tree); void Get_List_Of_Target_Edges(t_node *a,t_node *d,t_edge **list,int *list_size,t_tree *tree); void Fix_All(t_tree *tree); void Record_Br_Len(t_tree *tree); void Restore_Br_Len(t_tree *tree); void Get_Dist_Btw_Edges(t_node *a,t_node *d,t_tree *tree); void Detect_Polytomies(t_edge *b,phydbl l_thresh,t_tree *tree); void Get_List_Of_Nodes_In_Polytomy(t_node *a,t_node *d,t_node ***list,int *size_list); void Check_Path(t_node *a,t_node *d,t_node *target,t_tree *tree); void Connect_Two_Nodes(t_node *a,t_node *d); void Get_List_Of_Adjacent_Targets(t_node *a,t_node *d,t_node ***node_list,t_edge ***edge_list,int *list_size); void Sort_List_Of_Adjacent_Targets(t_edge ***list,int list_size); t_node *Common_Nodes_Btw_Two_Edges(t_edge *a,t_edge *b); int KH_Test(phydbl *site_lk_M1,phydbl *site_lk_M2,t_tree *tree); void Random_Tree(t_tree *tree); void Reorganize_Edges_Given_Lk_Struct(t_tree *tree); void Random_NNI(int n_moves,t_tree *tree); void Fill_Missing_Dist(matrix *mat); void Fill_Missing_Dist_XY(int x,int y,matrix *mat); phydbl Least_Square_Missing_Dist_XY(int x,int y,phydbl dxy,matrix *mat); void Check_Memory_Amount(t_tree *tree); int Get_State_From_P_Lk(phydbl *p_lk,int pos,t_tree *tree); int Get_State_From_P_Pars(short int *p_pars,int pos,t_tree *tree); void Check_Dirs(t_tree *tree); void Warn_And_Exit(char *s); void Randomize_Sequence_Order(calign *cdata); void Update_Root_Pos(t_tree *tree); void Add_Root(t_edge *target,t_tree *tree); void Update_Ancestors(t_node *a,t_node *d,t_tree *tree); #if (defined PHYTIME || defined SERGEII) t_tree *Generate_Random_Tree_From_Scratch(int n_otu,int rooted); #endif void Random_Lineage_Rates(t_node *a,t_node *d,t_edge *b,phydbl stick_prob,phydbl *rates,int curr_rate,int n_rates,t_tree *tree); t_edge *Find_Edge_With_Label(char *label,t_tree *tree); void Evolve(calign *data,t_mod *mod,t_tree *tree); int Pick_State(int n,phydbl *prob); void Evolve_Recur(t_node *a,t_node *d,t_edge *b,int a_state,int r_class,int site_num,calign *gen_data,t_mod *mod,t_tree *tree); void Site_Diversity(t_tree *tree); void Site_Diversity_Post(t_node *a,t_node *d,t_edge *b,t_tree *tree); void Site_Diversity_Pre(t_node *a,t_node *d,t_edge *b,t_tree *tree); void Subtree_Union(t_node *n,t_edge *b_fcus,t_tree *tree); void Binary_Decomposition(int value,int *bit_vect,int size); void Print_Diversity_Header(FILE *fp,t_tree *tree); phydbl Univariate_Kernel_Density_Estimate(phydbl where,phydbl *x,int sample_size); phydbl Multivariate_Kernel_Density_Estimate(phydbl *where,phydbl **x,int sample_size,int vect_size); phydbl Var(phydbl *x,int n); phydbl Mean(phydbl *x,int n); void Best_Of_NNI_And_SPR(t_tree *tree); int Polint(phydbl *xa,phydbl *ya,int n,phydbl x,phydbl *y,phydbl *dy); void JF(t_tree *tree); t_tree *Dist_And_BioNJ(calign *cdata,t_mod *mod,option *io); void Add_BioNJ_Branch_Lengths(t_tree *tree,calign *cdata,t_mod *mod); char *Bootstrap_From_String(char *s_tree,calign *cdata,t_mod *mod,option *io); char *aLRT_From_String(char *s_tree,calign *cdata,t_mod *mod,option *io); void Prepare_Tree_For_Lk(t_tree *tree); void Find_Common_Tips(t_tree *tree1,t_tree *tree2); phydbl Get_Tree_Size(t_tree *tree); void Dist_To_Root_Pre(t_node *a,t_node *d,t_edge *b,t_tree *tree); void Dist_To_Root(t_node *n_root,t_tree *tree); char *Basename(char *path); t_node *Find_Lca_Pair_Of_Nodes(t_node *n1,t_node *n2,t_tree *tree); t_node *Find_Lca_Clade(t_node **node_list,int node_list_size,t_tree *tree); int Get_List_Of_Ancestors(t_node *ref_node,t_node **list,int *size,t_tree *tree); int Edge_Num_To_Node_Num(int edge_num,t_tree *tree); void Time_To_Branch(t_tree *tree); void Time_To_Branch_Pre(t_node *a,t_node *d,t_tree *tree); void Branch_Lengths_To_Rate_Lengths(t_tree *tree); void Branch_Lengths_To_Rate_Lengths_Pre(t_node *a,t_node *d,t_tree *tree); int Find_Clade(char **tax_name_list,int list_size,t_tree *tree); void Find_Clade_Pre(t_node *a,t_node *d,int *tax_num_list,int list_size,int *num,t_tree *tree); t_edge *Find_Root_Edge(FILE *fp_input_tree,t_tree *tree); void Copy_Tree_Topology_With_Labels(t_tree *ori,t_tree *cpy); void Set_Model_Name(t_mod *mod); void Adjust_Min_Diff_Lk(t_tree *tree); void Translate_Tax_Names(char **tax_names,t_tree *tree); void Skip_Comment(FILE *fp); void Get_Best_Root_Position(t_tree *tree); void Get_Best_Root_Position_Post(t_node *a,t_node *d,int *has_outgrp,t_tree *tree); void Get_Best_Root_Position_Pre(t_node *a,t_node *d,t_tree *tree); void Get_OutIn_Scores(t_node *a,t_node *d); int Check_Sequence_Name(char *s); int Scale_Subtree_Height(t_node *a,phydbl K,phydbl floor,int *n_nodes,t_tree *tree); void Scale_Node_Heights_Post(t_node *a,t_node *d,phydbl K,phydbl floor,int *n_nodes,t_tree *tree); int Scale_Subtree_Rates(t_node *a,phydbl mult,int *n_nodes,t_tree *tree); void Check_Br_Len_Bounds(t_tree *tree); int Scale_Subtree_Rates_Post(t_node *a,t_node *d,phydbl mult,int *n_nodes,t_tree *tree); void Get_Node_Ranks(t_tree *tree); void Get_Node_Ranks_Pre(t_node *a,t_node *d,t_tree *tree); void Log_Br_Len(t_tree *tree); phydbl Diff_Lk_Norm_At_Given_Edge(t_edge *b,t_tree *tree); void Adjust_Variances(t_tree *tree); phydbl Effective_Sample_Size(phydbl first_val,phydbl last_val,phydbl sum,phydbl sumsq,phydbl sumcurnext,int n); void Rescale_Free_Rate_Tree(t_tree *tree); phydbl Rescale_Br_Len_Multiplier_Tree(t_tree *tree); phydbl Unscale_Br_Len_Multiplier_Tree(t_tree *tree); phydbl Reflect(phydbl x,phydbl l,phydbl u); int Are_Equal(phydbl a,phydbl b,phydbl eps); int Check_Topo_Constraints(t_tree *big_tree,t_tree *small_tree); void Prune_Tree(t_tree *big_tree,t_tree *small_tree); void Match_Nodes_In_Small_Tree(t_tree *small_tree,t_tree *big_tree); void Find_Surviving_Edges_In_Small_Tree(t_tree *small_tree,t_tree *big_tree); void Find_Surviving_Edges_In_Small_Tree_Post(t_node *a,t_node *d,t_tree *small_tree,t_tree *big_tree); void Set_Taxa_Id_Ranking(t_tree *tree); void Get_Edge_Binary_Coding_Number(t_tree *tree); void Make_Ancestral_Seq(t_tree *tree); void Make_MutMap(t_tree *tree); int Get_Mutmap_Val(int edge,int site,int mut,t_tree *tree); void Get_Mutmap_Coord(int idx,int *edge,int *site,int *mut,t_tree *tree); void Copy_Edge_Lengths(t_tree *to,t_tree *from); void Init_Scalar_Dbl(scalar_dbl *p); void Init_Scalar_Int(scalar_int *p); void Init_Vect_Dbl(int len,vect_dbl *p); void Init_Vect_Int(int len,vect_int *p); char *To_Lower_String(char *in); phydbl String_To_Dbl(char *string); char *To_Upper_String(char *in); void Connect_CSeqs_To_Nodes(calign *cdata, t_tree *tree); void Switch_Eigen(int state, t_mod *mod); void Joint_Proba_States_Left_Right(phydbl *Pij, phydbl *p_lk_left, phydbl *p_lk_rght, vect_dbl *pi, int scale_left, int scale_rght, phydbl *F, int n, int site, t_tree *tree); void Set_Both_Sides(int yesno, t_tree *tree); void Set_D_States(calign *data, int datatype, int stepsize); void Branch_To_Time(t_tree *tree); void Branch_To_Time_Pre(t_node *a, t_node *d, t_tree *tree); void Path_Length(t_node *dep, t_node *arr, phydbl *len, t_tree *tree); phydbl *Dist_Btw_Tips(t_tree *tree); void Random_SPRs_On_Rooted_Tree(t_tree *tree); void Set_P_Lk_One_Side(phydbl **Pij, phydbl **p_lk, int **sum_scale, t_node *d, t_edge *b, t_tree *tree); void Set_All_P_Lk(t_node **n_v1, t_node **n_v2, phydbl **p_lk , int **sum_scale , int **p_lk_loc, phydbl **Pij1, phydbl **p_lk1, int **sum_scale1, phydbl **Pij2, phydbl **p_lk2, int **sum_scale2, t_node *d, t_edge *b, t_tree *tree); void Optimum_Root_Position_IL_Model(t_tree *tree); void Set_Br_Len_Var(t_tree *tree); void Check_Br_Lens(t_tree *tree); #include "xml.h" #include "free.h" #include "spr.h" #include "lk.h" #include "optimiz.h" #include "models.h" #include "bionj.h" #include "simu.h" #include "eigen.h" #include "pars.h" #include "alrt.h" #include "stats.h" #include "help.h" #include "io.h" #include "make.h" #include "nexus.h" #include "init.h" #include "geo.h" #ifdef MPI #include "mpi_boot.h" #endif #ifdef MG #include "mg.h" #endif #ifdef TIME #include "times.h" #include "rates.h" #endif #ifdef _NOT_NEEDED_A_PRIORI #include "m4.h" #endif #endif