/* 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 "init.h" void Init_Eigen_Struct(eigen *this) { this->next = NULL; this->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Scalar_Dbl(scalar_dbl *p) { p->v = -1.; p->onoff = ON; p->next = NULL; p->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Scalar_Int(scalar_int *p) { p->v = -1.; p->next = NULL; p->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Vect_Dbl(int len, vect_dbl *p) { p->len = len; p->next = NULL; p->prev = NULL; p->v = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Vect_Int(int len, vect_int *p) { p->len = len; p->next = NULL; p->prev = NULL; p->v = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_String(t_string *ts) { ts->len = -1.; ts->next = NULL; ts->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Triplet_Struct(triplet *t) { t->next = NULL; t->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Efrq(t_efrq *f) { f->next = NULL; f->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Tree(t_tree *tree, int n_otu) { tree->n_otu = n_otu; tree->mat = NULL; tree->n_root = NULL; tree->e_root = NULL; tree->ps_tree = NULL; tree->short_l = NULL; tree->mutmap = NULL; tree->next = NULL; tree->prev = NULL; tree->next = NULL; tree->prev = NULL; tree->mixt_tree = NULL; tree->geo = NULL; tree->is_mixt_tree = NO; tree->tree_num = 0; tree->depth_curr_path = 0; tree->has_bip = NO; tree->n_moves = 0; tree->n_improvements = 0; tree->bl_from_node_stamps = 0; tree->lock_topo = NO; tree->ps_page_number = 0; tree->init_lnL = UNLIKELY; tree->best_lnL = UNLIKELY; tree->old_lnL = UNLIKELY; tree->c_lnL = UNLIKELY; tree->sum_min_sum_scale = .0; tree->n_swap = 0; tree->best_pars = 1E+5; tree->n_pattern = -1; tree->n_root_pos = -1.; tree->write_labels = YES; tree->write_br_lens = YES; tree->print_boot_val = 0; tree->print_alrt_val = 0; tree->num_curr_branch_available = 0; tree->tip_order_score = .0; tree->write_tax_names = YES; tree->update_alias_subpatt = NO; tree->bl_ndigits = 8; tree->n_short_l = 100; tree->norm_scale = 0.0; tree->br_len_recorded = NO; tree->max_spr_depth = 0; tree->apply_lk_scaling = YES; tree->dp = 0; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Edge_Light(t_edge *b, int num) { b->num = num; b->bip_score = 0; b->dist_btw_edges = .0; b->topo_dist_btw_edges = 0; b->has_zero_br_len = NO; b->n_jumps = 0; b->l_var->v = -1.; b->does_exist = YES; b->l->v = -1.; b->bin_cod_num = -1.; b->l->onoff = ON; b->next = NULL; b->prev = NULL; b->next_mixt = NULL; b->prev_mixt = NULL; b->p_lk_left = NULL; b->p_lk_rght = NULL; b->p_lk_loc_left = NULL; b->p_lk_loc_rght = NULL; b->Pij_rr = NULL; b->labels = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Node_Light(t_node *n, int num) { n->num = num; n->tax = -1; n->dist_to_root = .0; n->common = 1; n->ext_node = NULL; n->name = NULL; n->ori_name = NULL; n->c_seq = NULL; n->c_seq_anc = NULL; n->y_rank = 0.; n->y_rank_ori = 0.; n->y_rank_max = 0.; n->y_rank_min = 0.; n->anc = NULL; n->rank = 0; n->match_node = NULL; n->id_rank = 0; n->next = NULL; n->prev = NULL; /* n->next = NULL; */ /* n->prev = NULL; */ } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_NNI(nni *a_nni) { a_nni->left = NULL; a_nni->rght = NULL; a_nni->b = NULL; a_nni->init_l = -1.; a_nni->init_lk = .0; a_nni->score = +1.0; a_nni->best_l = -1.; a_nni->swap_node_v1 = NULL; a_nni->swap_node_v2 = NULL; a_nni->swap_node_v3 = NULL; a_nni->swap_node_v4 = NULL; a_nni->lk0 = UNLIKELY; a_nni->lk1 = UNLIKELY; a_nni->lk2 = UNLIKELY; a_nni->l0 = -1.0; a_nni->l1 = -1.0; a_nni->l2 = -1.0; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Nexus_Format(nexcom **com) { /*****************************/ strcpy(com[0]->name,"dimensions"); com[0]->nparm = 2; com[0]->nxt_token_t = NEXUS_PARM; com[0]->cur_token_t = NEXUS_COM; com[0]->parm[0] = Make_Nexus_Parm(); strcpy(com[0]->parm[0]->name,"ntax"); com[0]->parm[0]->fp = Read_Nexus_Dimensions; com[0]->parm[0]->com = com[0]; com[0]->parm[0]->nxt_token_t = NEXUS_EQUAL; com[0]->parm[0]->cur_token_t = NEXUS_PARM; com[0]->parm[1] = Make_Nexus_Parm(); strcpy(com[0]->parm[1]->name,"nchar"); com[0]->parm[1]->fp = Read_Nexus_Dimensions; com[0]->parm[1]->com = com[0]; com[0]->parm[1]->nxt_token_t = NEXUS_EQUAL; com[0]->parm[1]->cur_token_t = NEXUS_PARM; /*****************************/ strcpy(com[1]->name,"format"); com[1]->nparm = 11; com[1]->nxt_token_t = NEXUS_PARM; com[1]->cur_token_t = NEXUS_COM; com[1]->parm[0] = Make_Nexus_Parm(); strcpy(com[1]->parm[0]->name,"datatype"); com[1]->parm[0]->fp = Read_Nexus_Format; com[1]->parm[0]->com = com[1]; com[1]->parm[0]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[0]->cur_token_t = NEXUS_PARM; com[1]->parm[1] = Make_Nexus_Parm(); strcpy(com[1]->parm[1]->name,"respectcase"); com[1]->parm[1]->fp = Read_Nexus_Format; com[1]->parm[1]->com = com[1]; com[1]->parm[1]->nxt_token_t = NEXUS_PARM; com[1]->parm[1]->cur_token_t = NEXUS_VALUE; com[1]->parm[2] = Make_Nexus_Parm(); strcpy(com[1]->parm[2]->name,"missing"); com[1]->parm[2]->fp = Read_Nexus_Format; com[1]->parm[2]->com = com[1]; com[1]->parm[2]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[2]->cur_token_t = NEXUS_PARM; com[1]->parm[3] = Make_Nexus_Parm(); strcpy(com[1]->parm[3]->name,"gap"); com[1]->parm[3]->fp = Read_Nexus_Format; com[1]->parm[3]->com = com[1]; com[1]->parm[3]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[3]->cur_token_t = NEXUS_PARM; com[1]->parm[4] = Make_Nexus_Parm(); strcpy(com[1]->parm[4]->name,"symbols"); com[1]->parm[4]->fp = Read_Nexus_Format; com[1]->parm[4]->com = com[1]; com[1]->parm[4]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[4]->cur_token_t = NEXUS_PARM; com[1]->parm[5] = Make_Nexus_Parm(); strcpy(com[1]->parm[5]->name,"equate"); com[1]->parm[5]->fp = Read_Nexus_Format; com[1]->parm[5]->com = com[1]; com[1]->parm[5]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[5]->cur_token_t = NEXUS_PARM; com[1]->parm[6] = Make_Nexus_Parm(); strcpy(com[1]->parm[6]->name,"matchchar"); com[1]->parm[6]->fp = Read_Nexus_Format; com[1]->parm[6]->com = com[1]; com[1]->parm[6]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[6]->cur_token_t = NEXUS_PARM; com[1]->parm[7] = Make_Nexus_Parm(); strcpy(com[1]->parm[7]->name,"transpose"); com[1]->parm[7]->fp = Read_Nexus_Format; com[1]->parm[7]->com = com[1]; com[1]->parm[7]->nxt_token_t = NEXUS_PARM; com[1]->parm[7]->cur_token_t = NEXUS_VALUE; com[1]->parm[8] = Make_Nexus_Parm(); strcpy(com[1]->parm[8]->name,"interleave"); com[1]->parm[8]->fp = Read_Nexus_Format; com[1]->parm[8]->com = com[1]; com[1]->parm[8]->nxt_token_t = NEXUS_PARM; com[1]->parm[8]->cur_token_t = NEXUS_VALUE; com[1]->parm[9] = Make_Nexus_Parm(); strcpy(com[1]->parm[9]->name,"items"); com[1]->parm[9]->fp = Read_Nexus_Format; com[1]->parm[9]->com = com[1]; com[1]->parm[9]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[9]->cur_token_t = NEXUS_PARM; com[1]->parm[10] = Make_Nexus_Parm(); strcpy(com[1]->parm[10]->name,"statesformat"); com[1]->parm[10]->fp = Read_Nexus_Format; com[1]->parm[10]->com = com[1]; com[1]->parm[10]->nxt_token_t = NEXUS_EQUAL; com[1]->parm[10]->cur_token_t = NEXUS_PARM; /*****************************/ strcpy(com[2]->name,"eliminate"); com[2]->nparm = 0; com[2]->nxt_token_t = NEXUS_VALUE; com[2]->cur_token_t = NEXUS_COM; /*****************************/ strcpy(com[3]->name,"taxlabels"); com[3]->nparm = 0; com[3]->nxt_token_t = -1; com[3]->cur_token_t = -1; /*****************************/ strcpy(com[4]->name,"charstatelabels"); com[4]->nparm = 0; com[4]->nxt_token_t = -1; com[4]->cur_token_t = -1; /*****************************/ strcpy(com[5]->name,"charlabels"); com[5]->nparm = 0; com[5]->nxt_token_t = -1; com[5]->cur_token_t = -1; /*****************************/ strcpy(com[6]->name,"statelabels"); com[6]->nparm = 0; com[6]->nxt_token_t = -1; com[6]->cur_token_t = -1; /*****************************/ strcpy(com[7]->name,"matrix"); com[7]->nparm = 1; com[7]->nxt_token_t = NEXUS_COM; com[7]->cur_token_t = NEXUS_VALUE; /* This will allow us to skip directly to the matrix reading function */ com[7]->parm[0] = Make_Nexus_Parm(); strcpy(com[7]->parm[0]->name,"matrix"); com[7]->parm[0]->fp = Read_Nexus_Matrix; com[7]->parm[0]->com = com[7]; com[7]->parm[0]->nxt_token_t = NEXUS_COM; com[7]->parm[0]->cur_token_t = -1; /*****************************/ strcpy(com[8]->name,"begin"); com[8]->nparm = 3; com[8]->nxt_token_t = NEXUS_PARM; com[8]->cur_token_t = NEXUS_COM; com[8]->parm[0] = Make_Nexus_Parm(); strcpy(com[8]->parm[0]->name,"data"); com[8]->parm[0]->fp = Read_Nexus_Begin; com[8]->parm[0]->com = com[8]; com[8]->parm[0]->nxt_token_t = NEXUS_COM; com[8]->parm[0]->cur_token_t = NEXUS_PARM; com[8]->parm[1] = Make_Nexus_Parm(); strcpy(com[8]->parm[1]->name,"trees"); com[8]->parm[1]->fp = Read_Nexus_Begin; com[8]->parm[1]->com = com[8]; com[8]->parm[1]->nxt_token_t = NEXUS_COM; com[8]->parm[1]->cur_token_t = NEXUS_PARM; com[8]->parm[2] = Make_Nexus_Parm(); strcpy(com[8]->parm[2]->name,"taxa"); com[8]->parm[2]->fp = Read_Nexus_Taxa; com[8]->parm[2]->com = com[8]; com[8]->parm[2]->nxt_token_t = NEXUS_COM; com[8]->parm[2]->cur_token_t = NEXUS_VALUE; /*****************************/ strcpy(com[9]->name,"end"); com[9]->nparm = 0; com[9]->nxt_token_t = -1; com[9]->cur_token_t = -1; /*****************************/ strcpy(com[10]->name,"translate"); com[10]->nparm = 1; com[10]->nxt_token_t = NEXUS_COM; com[10]->cur_token_t = NEXUS_VALUE; com[10]->parm[0] = Make_Nexus_Parm(); strcpy(com[10]->parm[0]->name,"translate"); com[10]->parm[0]->fp = Read_Nexus_Translate; com[10]->parm[0]->com = com[10]; com[10]->parm[0]->nxt_token_t = NEXUS_COM; com[10]->parm[0]->cur_token_t = -1; /*****************************/ strcpy(com[11]->name,"tree"); com[11]->nparm = 1; com[11]->nxt_token_t = NEXUS_COM; com[11]->cur_token_t = NEXUS_VALUE; com[11]->parm[0] = Make_Nexus_Parm(); strcpy(com[11]->parm[0]->name,"tree"); com[11]->parm[0]->fp = Read_Nexus_Tree; com[11]->parm[0]->com = com[11]; com[11]->parm[0]->nxt_token_t = -1; com[11]->parm[0]->cur_token_t = -1; /*****************************/ } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Mat(matrix *mat, calign *data) { int i; mat->n_otu = data->n_otu; mat->r = mat->n_otu; mat->curr_int = mat->n_otu; mat->method = 1; For(i,data->n_otu) { strcpy(mat->name[i],data->c_seq[i]->name); mat->on_off[i] = 1; } } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Set_Defaults_Input(option* io) { io->fp_in_align = NULL; io->fp_in_tree = NULL; io->fp_in_constraint_tree = NULL; io->fp_out_tree = NULL; io->fp_out_trees = NULL; io->fp_out_boot_tree = NULL; io->fp_out_boot_stats = NULL; io->fp_out_stats = NULL; io->long_tax_names = NULL; io->short_tax_names = NULL; io->lon = NULL; io->lat = NULL; io->z_scores = NULL; io->cstr_tree = NULL; io->next = NULL; io->prev = NULL; io->tree = NULL; io->mod = NULL; strcpy(io->nt_or_cd,"nucleotides"); io->n_data_sets = 1; io->interleaved = 1; io->in_tree = 0; io->out_tree_file_open_mode = 1; io->out_stats_file_open_mode = 1; io->init_len = -1; io->n_otu = -1; io->n_data_set_asked = -1; io->print_boot_trees = 1; io->n_part = 1; io->ratio_test = ABAYES; io->multigene = 0; io->config_multigene = 0; io->curr_interface = 0; io->r_seed = -1; io->collapse_boot = 0; io->random_boot_seq_order = 1; io->print_trace = 0; io->print_site_lnl = 0; io->m4_model = NO; io->rm_ambigu = 0; io->append_run_ID = 0; io->quiet = 0; io->datatype = NT; io->colalias = YES; io->data_file_format = PHYLIP; io->tree_file_format = PHYLIP; io->boot_prog_every = 20; io->mem_question = YES; io->do_alias_subpatt = NO; io->lk_approx = EXACT; io->codpos = -1; io->mutmap = NO; io->state_len = 1; MCMC_Init_MCMC_Struct(NULL,io,io->mcmc); RATES_Init_Rate_Struct(io->rates,NULL,-1); io->rates->model = GUINDON; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Rmat(t_rmat *rmat) { rmat->n_diff_rr = 1; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Set_Defaults_Model(t_mod *mod) { strcpy(mod->modelname->s,"HKY85"); strcpy(mod->custom_mod_string->s,"000000"); mod->next = NULL; mod->prev = NULL; mod->next_mixt = NULL; mod->prev_mixt = NULL; mod->r_mat = NULL; mod->e_frq = NULL; mod->whichmodel = HKY85; mod->ras->n_catg = 4; mod->n_mixt_classes = 0; mod->mod_num = 0; mod->update_eigen = NO; mod->is_mixt_mod = NO; mod->ras->normalise_rr = YES; mod->kappa->v = 4.0; mod->ras->alpha->v = 1.0; mod->lambda->v = 1.0; mod->ras->pinvar->v = 0.0; mod->l_var_sigma = 1.E-2; mod->e_frq_weight->v = 1.0; mod->r_mat_weight->v = 1.0; mod->bootstrap = 0; mod->ras->invar = NO; mod->ns = 4; mod->use_m4mod = NO; mod->ras->gamma_median = NO; mod->m4mod = NULL; /* mod->r_mat->n_diff_rr = 0; */ /* mod->r_mat->rr = NULL; */ /* mod->r_mat->rr_val = NULL; */ /* mod->r_mat->n_rr_per_cat = NULL; */ mod->io = NULL; mod->log_l = NO; mod->ras->free_mixt_rates = NO; mod->gamma_mgf_bl = NO; mod->br_len_multiplier->v = 1.0; mod->ras->parent_class_number = 0; mod->ras->init_r_proba = YES; mod->ras->init_rr = YES; #if !(defined PHYTIME || defined SERGEII) mod->l_min = 1.E-8; mod->l_max = 100.0; #else mod->l_min = 1.E-8; mod->l_max = 2.0; #endif mod->l_var_min = mod->l_min; mod->l_var_max = mod->l_max; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Set_Defaults_Optimiz(t_opt *s_opt) { s_opt->print = YES; s_opt->last_opt = YES; s_opt->opt_subst_param = YES; s_opt->opt_alpha = YES; s_opt->opt_kappa = YES; s_opt->opt_bl = YES; s_opt->opt_lambda = NO; s_opt->opt_pinvar = NO; s_opt->opt_cov_delta = NO; s_opt->opt_cov_alpha = NO; s_opt->opt_cov_free_rates = NO; s_opt->opt_rr = NO; s_opt->init_lk = UNLIKELY; s_opt->n_it_max = 1000; s_opt->opt_topo = YES; s_opt->topo_search = NNI_MOVE; s_opt->random_input_tree = 0; s_opt->n_rand_starts = 5; s_opt->brent_it_max = BRENT_IT_MAX; s_opt->steph_spr = YES; s_opt->user_state_freq = NO; s_opt->min_diff_lk_local = 1.E-04; s_opt->min_diff_lk_global = 1.E-03; s_opt->min_diff_lk_move = 1.E-02; s_opt->p_moves_to_examine = 0.15; s_opt->fast_nni = NO; s_opt->greedy = NO; s_opt->general_pars = NO; s_opt->tree_size_mult = 1; s_opt->opt_five_branch = YES; s_opt->pars_thresh = 5; s_opt->hybrid_thresh = NO; s_opt->quickdirty = NO; s_opt->spr_pars = YES; s_opt->spr_lnL = NO; s_opt->min_depth_path = 0; s_opt->max_depth_path = 20; s_opt->deepest_path = 20; s_opt->max_delta_lnL_spr = 50.; s_opt->br_len_in_spr = 10; s_opt->opt_free_mixt_rates = YES; s_opt->constrained_br_len = NO; s_opt->opt_gamma_br_len = NO; s_opt->first_opt_free_mixt_rates = YES; s_opt->wim_n_rgrft = -1; s_opt->wim_n_globl = -1; s_opt->wim_max_dist = -1; s_opt->wim_n_optim = -1; s_opt->wim_n_best = -1; s_opt->wim_inside_opt = 0; s_opt->opt_rmat_weight = NO; s_opt->opt_efrq_weight = NO; s_opt->skip_tree_traversal = NO; s_opt->serial_free_rates = YES; s_opt->curr_opt_free_rates = NO; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void XML_Init_Attribute(xml_attr *attr) { } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void XML_Init_Node(xml_node *parent, xml_node *new_node, char *name) { if(name) strcpy(new_node->name,name); new_node->parent = parent ? parent : NULL; new_node->next = NULL; new_node->prev = NULL; new_node->child = NULL; new_node->ds->obj = NULL; new_node->ds->next = NULL; if(parent) { if(!parent->child) { parent->child = new_node; } else { xml_node *node = parent->child; while(node->next) node = node->next; node->next = new_node; new_node->prev = node; } } new_node->attr = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void RATES_Init_Rate_Struct(t_rate *rates, t_rate *existing_rates, int n_otu) { int i; if(existing_rates && (existing_rates->model != -1)) { rates->model = existing_rates->model; } else { rates->model = NONE; } if(rates->model == NONE) rates->model_log_rates = NO; else if(rates->model == THORNE) rates->model_log_rates = YES; else if(rates->model == GUINDON) rates->model_log_rates = YES; else if(rates->model == GAMMA) rates->model_log_rates = NO; else if(rates->model == STRICTCLOCK) rates->model_log_rates = NO; else { PhyML_Printf("\n== Please initialize model properly."); PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__); Warn_And_Exit(""); } rates->met_within_gibbs = NO; rates->c_lnL_rates = UNLIKELY; rates->c_lnL_jps = UNLIKELY; rates->adjust_rates = 0; rates->use_rates = 1; rates->lexp = 1.E-3; rates->norm_fact = 1.0; rates->inflate_var = 1.0; rates->nd_t_recorded = NO; rates->br_r_recorded = NO; rates->birth_rate = 1.E-2; rates->birth_rate_min = 1.E-6; rates->birth_rate_max = 1.E+1; if(rates->model_log_rates == YES) { rates->max_rate = LOG(10.); rates->min_rate = -LOG(10.); /* rates->max_rate = MDBL_MAX; */ /* rates->min_rate = -MDBL_MAX; */ } else { rates->max_rate = 10.0; rates->min_rate = 0.0; } /* rates->max_rate = 6.0; */ /* rates->min_rate = 0.0; */ rates->clock_r = 1.E-4; rates->min_clock = 1.E-8; rates->max_clock = 1.E-3; /* rates->clock_r = 3.E-4; */ /* rates->max_clock = 1.E-3; */ /* rates->min_clock = 1.E-5; */ if(rates->model != GAMMA) { rates->nu = 1.E-3; rates->min_nu = 0.0; rates->max_nu = 1.0; } else { rates->nu = 1.E-1; rates->min_nu = 1.E-2; rates->max_nu = 100.0; } /* rates->max_nu = 2.0; */ /* rates->nu = 1.E-4; */ /* rates->max_nu = 1.E-1; */ /* rates->min_nu = 1.E-5; */ rates->min_dt = 0.0; rates->step_rate = 1.E-4; rates->approx = 1; rates->bl_from_rt = 0; rates->update_mean_l = NO; rates->update_cov_l = NO; rates->p_max = 0.01; rates->true_tree_size = 0.0; if(n_otu > 0) { For(i,(2*n_otu-2)*(2*n_otu-2)) rates->cov_l[i] = 0.0; For(i,2*n_otu-2) { rates->n_jps[i] = -1; rates->t_jps[i] = -1; rates->mean_r[i] = 1.0; rates->mean_l[i] = 0.0; } For(i,2*n_otu-1) { if(rates->model_log_rates == YES) { rates->nd_r[i] = 0.0; rates->br_r[i] = 0.0; } else { rates->nd_r[i] = 1.0; rates->br_r[i] = 1.0; } rates->mean_t[i] = 0.0; rates->nd_t[i] = 0.0; rates->true_t[i] = 0.0; if(i < n_otu) { rates->t_has_prior[i] = YES; rates->t_prior_max[i] = 0.0; rates->t_prior_min[i] = 0.0; } else { rates->t_has_prior[i] = NO; rates->t_prior_max[i] = BIG; rates->t_prior_min[i] = -BIG; } rates->br_do_updt[i] = YES; rates->has_survived[i] = NO; rates->t_ranked[i] = i; } } rates->calib = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_One_Spr(t_spr *a_spr) { a_spr->lnL = UNLIKELY; a_spr->pars = 1E+5; a_spr->depth_path = 0; a_spr->dist = 0; a_spr->init_target_l = -1.; a_spr->init_target_v = -1.; a_spr->l0 = -1.; a_spr->l1 = -1.; a_spr->l2 = -1.; a_spr->v0 = -1.; a_spr->v1 = -1.; a_spr->v2 = -1.; a_spr->n_link = NULL; a_spr->n_opp_to_link = NULL; a_spr->b_opp_to_link = NULL; a_spr->b_target = NULL; a_spr->b_init_target = NULL; a_spr->next = NULL; a_spr->prev = NULL; a_spr->next = NULL; a_spr->prev = NULL; } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// void Init_Model(calign *data, t_mod *mod, option *io) { int i,j; phydbl sum,aux; int result; phydbl *dr, *di, *space; mod->ns = io->mod->ns; if(io->datatype == GENERIC) mod->whichmodel = JC69; /* if(!mod->ras->invar) For(i,data->crunch_len) data->invar[i] = 0; */ dr = (phydbl *)mCalloc( mod->ns,sizeof(phydbl)); di = (phydbl *)mCalloc( mod->ns,sizeof(phydbl)); space = (phydbl *)mCalloc(2*mod->ns,sizeof(phydbl)); if(mod->log_l == YES) { mod->l_min = LOG(mod->l_min); mod->l_max = LOG(mod->l_max); } if(mod->ras->init_r_proba == YES) { For(i,mod->ras->n_catg) mod->ras->gamma_r_proba->v[i] = (phydbl)1./mod->ras->n_catg; For(i,mod->ras->n_catg) mod->ras->gamma_r_proba_unscaled->v[i] = (phydbl)(i+1); } else { mod->ras->gamma_r_proba_unscaled->v[0] = mod->ras->gamma_r_proba->v[0]; for(i=1;iras->n_catg;i++) mod->ras->gamma_r_proba_unscaled->v[i] = mod->ras->gamma_r_proba_unscaled->v[i-1] + mod->ras->gamma_r_proba->v[i]; } if(mod->ras->init_rr == YES) { if(mod->ras->n_catg > 1) { For(i,mod->ras->n_catg) mod->ras->gamma_rr->v[i] = (phydbl)i; For(i,mod->ras->n_catg) mod->ras->gamma_rr_unscaled->v[i] = (phydbl)i; } else { mod->ras->gamma_rr->v[0] = 1.0; mod->ras->gamma_rr_unscaled->v[0] = 1.0; } } mod->br_len_multiplier->v = 1.0; For(i,mod->ns) { mod->e_frq->pi->v[i] = data->b_frq[i]; mod->e_frq->pi_unscaled->v[i] = mod->e_frq->pi->v[i] * 100.; } if(io->datatype == NT) { /* Set the substitution parameters to their default values if they are not fixed by the user */ if(mod->s_opt->opt_kappa == YES) { mod->kappa->v = 4.0; mod->lambda->v = 1.0; } /* if(mod->s_opt->opt_rr) */ if(mod->whichmodel == GTR || mod->whichmodel == CUSTOM) { int i; For(i,6) { mod->r_mat->rr->v[i] = 1.0; mod->r_mat->rr_val->v[i] = 1.0; } } } if(mod->s_opt->opt_alpha) mod->ras->alpha->v = 1.0; if(mod->s_opt->opt_pinvar) mod->ras->pinvar->v = 0.2; if(io->datatype == NT) /* Nucleotides */ { /* init for nucleotides */ mod->lambda->v = 1.; /* mod->update_eigen = YES; */ if(mod->whichmodel == JC69) { mod->e_frq->pi->v[0] = mod->e_frq->pi->v[1] = mod->e_frq->pi->v[2] = mod->e_frq->pi->v[3] = .25; mod->kappa->v = 1.; mod->s_opt->opt_state_freq = NO; mod->s_opt->opt_kappa = NO; mod->s_opt->opt_lambda = NO; mod->update_eigen = NO; } if(mod->whichmodel == K80) { mod->e_frq->pi->v[0] = mod->e_frq->pi->v[1] = mod->e_frq->pi->v[2] = mod->e_frq->pi->v[3] = .25; mod->s_opt->opt_state_freq = NO; mod->s_opt->opt_lambda = NO; mod->update_eigen = NO; } if(mod->whichmodel == F81) { mod->kappa->v = 1.; mod->update_eigen = NO; } if(mod->whichmodel == F84) { aux = ((mod->e_frq->pi->v[0]+mod->e_frq->pi->v[2])-(mod->e_frq->pi->v[1]+mod->e_frq->pi->v[3]))/(2.*mod->kappa->v); mod->lambda->v = ((mod->e_frq->pi->v[1]+mod->e_frq->pi->v[3]) + aux)/((mod->e_frq->pi->v[0]+mod->e_frq->pi->v[2]) - aux); mod->update_eigen = NO; } if(mod->whichmodel == TN93) { mod->update_eigen = NO; if(io->mod->s_opt->opt_kappa) io->mod->s_opt->opt_lambda = 1; } if(mod->whichmodel == GTR) { mod->kappa->v = 1.; mod->update_eigen = YES; io->mod->s_opt->opt_rr = YES; } if(mod->whichmodel == CUSTOM) { mod->kappa->v = 1.; mod->update_eigen = YES; /* io->mod->s_opt->opt_rr = YES; */ /* What if the user decided not to optimise the rates? */ } if(mod->whichmodel == GTR) { mod->custom_mod_string->s[0] = '0'; mod->custom_mod_string->s[1] = '1'; mod->custom_mod_string->s[2] = '2'; mod->custom_mod_string->s[3] = '3'; mod->custom_mod_string->s[4] = '4'; mod->custom_mod_string->s[5] = '5'; Translate_Custom_Mod_String(mod); } if(mod->s_opt->user_state_freq == YES) { For(i,4) { mod->e_frq->pi->v[i] = mod->user_b_freq->v[i]; } } if(((mod->whichmodel == GTR) || (mod->whichmodel == CUSTOM) || (mod->whichmodel == HKY85)) && (mod->use_m4mod == NO)) { mod->update_eigen = YES; Update_Eigen(mod); mod->update_eigen = NO; } /* if((mod->whichmodel != GTR) && */ /* (mod->whichmodel != CUSTOM) && */ /* (mod->whichmodel != HKY85)) mod->update_eigen = NO; */ } else if(mod->io->datatype == AA) { /* init for amino-acids */ /* see comments of PMat_Empirical for details */ /* read pi and Q from file */ /* These initialisations are needed when analysing multiple * data sets */ For(i,mod->ns*mod->ns) mod->r_mat->qmat->v[i] = .0; For(i,mod->ns ) mod->e_frq->pi->v[i] = .0; switch(mod->whichmodel) { case DAYHOFF : { Init_Qmat_Dayhoff(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case JTT : { Init_Qmat_JTT(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case MTREV : { Init_Qmat_MtREV(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case LG : { Init_Qmat_LG(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case WAG : { Init_Qmat_WAG(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case DCMUT : { Init_Qmat_DCMut(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case RTREV : { Init_Qmat_RtREV(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case CPREV : { Init_Qmat_CpREV(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case VT : { Init_Qmat_VT(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case BLOSUM62 : { Init_Qmat_Blosum62(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case MTMAM : { Init_Qmat_MtMam(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case MTART : { Init_Qmat_MtArt(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case HIVW : { Init_Qmat_HIVw(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case HIVB : { Init_Qmat_HIVb(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } case CUSTOMAA : { Read_Qmat(mod->r_mat->qmat->v,mod->e_frq->pi->v,mod->fp_aa_rate_mat); /* Print_Qmat_AA(mod->r_mat->qmat->v,mod->e_frq->pi->v); */ if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } default : { Init_Qmat_LG(mod->r_mat->qmat->v,mod->e_frq->pi->v); if(mod->s_opt->opt_state_freq) For(i,mod->ns) mod->e_frq->pi->v[i] = data->b_frq[i]; break; } } /* /\* multiply the nth col of Q by the nth term of pi/100 just as in PAML *\/ */ For(i,mod->ns) For(j,mod->ns) mod->r_mat->qmat->v[i*mod->ns+j] *= mod->e_frq->pi->v[j] / 100.0; /* compute diagonal terms of Q and mean rate mr = l/t */ mod->mr->v= .0; For (i,mod->ns) { sum=.0; For(j, mod->ns) sum += mod->r_mat->qmat->v[i*mod->ns+j]; mod->r_mat->qmat->v[i*mod->ns+i] = -sum; mod->mr->v += mod->e_frq->pi->v[i] * sum; } /* scale imod->nstantaneous rate matrix so that mu=1 */ For (i,mod->ns*mod->ns) mod->r_mat->qmat->v[i] /= mod->mr->v; /* compute eigenvectors/values */ result = 0; For(i,mod->ns*mod->ns) mod->r_mat->qmat_buff->v[i] = mod->r_mat->qmat->v[i]; if(!Eigen(1,mod->r_mat->qmat_buff->v,mod->eigen->size,mod->eigen->e_val, mod->eigen->e_val_im,mod->eigen->r_e_vect, mod->eigen->r_e_vect_im,mod->eigen->space)) { /* compute inverse(Vr) into Vi */ For (i,mod->ns*mod->ns) mod->eigen->l_e_vect[i] = mod->eigen->r_e_vect[i]; if(!Matinv(mod->eigen->l_e_vect,mod->eigen->size,mod->eigen->size,YES)) { PhyML_Printf("\n== Err in file %s at line %d.",__FILE__,__LINE__); Exit("\n"); } /* compute the diagonal terms of EXP(D) */ For(i,mod->ns) mod->eigen->e_val[i] = (phydbl)EXP(mod->eigen->e_val[i]); } else { if (result==-1) PhyML_Printf("\n== Eigenvalues/vectors computation does not converge : computation cancelled"); else if (result==1) PhyML_Printf("\n== Complex eigenvalues/vectors : computation cancelled"); } } else if(mod->io->datatype == GENERIC) { /* Uniform state frequencies */ For(i,mod->ns) mod->e_frq->pi->v[i] = 1./(phydbl)mod->ns; mod->kappa->v = 1; mod->s_opt->opt_state_freq = NO; mod->s_opt->opt_kappa = NO; mod->s_opt->opt_lambda = NO; mod->update_eigen = NO; } else { PhyML_Printf("\n== Datatype not recognized.\n"); PhyML_Printf("\n== Err. in file %s at line %d\n",__FILE__,__LINE__); Warn_And_Exit(""); } if(!mod->use_m4mod) Set_Model_Parameters(mod); Init_Eigen_Struct(mod->eigen); free(dr);free(di);free(space); } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// int Init_Qmat_Dayhoff(phydbl *daa, phydbl *pi) { /* Dayhoff's model data * Dayhoff, M.O., Schwartz, R.M., Orcutt, B.C. (1978) * "A model of evolutionary change in proteins." * Dayhoff, M.O.(ed.) Atlas of Protein Sequence Structur., Vol5, Suppl3. * National Biomedical Research Foundation, Washington DC, pp.345-352. */ int i,j,naa; naa = 20; /* PhyML_Printf("\n\n. REMINDER : THIS IS NOT DAYHOFF !!!\n\n"); */ /* daa[1*20 + 0] = 0.538903; */ /* daa[2*20 + 0] = 0.412504; */ /* daa[2*20 + 1] = 0.736081; */ /* daa[3*20 + 0] = 0.586915; */ /* daa[3*20 + 1] = 0.108051; */ /* daa[3*20 + 2] = 5.446642; */ /* daa[4*20 + 0] = 2.189718; */ /* daa[4*20 + 1] = 0.830604; */ /* daa[4*20 + 2] = 0.573426; */ /* daa[4*20 + 3] = 0.077565; */ /* daa[5*20 + 0] = 1.08213; */ /* daa[5*20 + 1] = 2.950693; */ /* daa[5*20 + 2] = 1.739514; */ /* daa[5*20 + 3] = 0.559035; */ /* daa[5*20 + 4] = 0.111314; */ /* daa[6*20 + 0] = 1.386865; */ /* daa[6*20 + 1] = 0.358434; */ /* daa[6*20 + 2] = 0.541447; */ /* daa[6*20 + 3] = 5.406871; */ /* daa[6*20 + 4] = 0.003738; */ /* daa[6*20 + 5] = 4.124043; */ /* daa[7*20 + 0] = 2.085747; */ /* daa[7*20 + 1] = 0.453132; */ /* daa[7*20 + 2] = 1.815844; */ /* daa[7*20 + 3] = 1.08647; */ /* daa[7*20 + 4] = 0.526418; */ /* daa[7*20 + 5] = 0.347693; */ /* daa[7*20 + 6] = 0.438476; */ /* daa[8*20 + 0] = 0.407898; */ /* daa[8*20 + 1] = 2.689322; */ /* daa[8*20 + 2] = 5.386808; */ /* daa[8*20 + 3] = 0.884563; */ /* daa[8*20 + 4] = 0.658583; */ /* daa[8*20 + 5] = 4.588358; */ /* daa[8*20 + 6] = 0.386218; */ /* daa[8*20 + 7] = 0.368668; */ /* daa[9*20 + 0] = 0.10177; */ /* daa[9*20 + 1] = 0.104875; */ /* daa[9*20 + 2] = 0.16239; */ /* daa[9*20 + 3] = 0.011698; */ /* daa[9*20 + 4] = 0.253282; */ /* daa[9*20 + 5] = 0.083872; */ /* daa[9*20 + 6] = 0.041767; */ /* daa[9*20 + 7] = 0.009778; */ /* daa[9*20 + 8] = 0.1042; */ /* daa[10*20 + 0] = 0.248202; */ /* daa[10*20 + 1] = 0.375742; */ /* daa[10*20 + 2] = 0.093863; */ /* daa[10*20 + 3] = 0.019241; */ /* daa[10*20 + 4] = 0.572688; */ /* daa[10*20 + 5] = 0.703538; */ /* daa[10*20 + 6] = 0.071961; */ /* daa[10*20 + 7] = 0.03006; */ /* daa[10*20 + 8] = 0.418222; */ /* daa[10*20 + 9] = 3.702051; */ /* daa[11*20 + 0] = 0.652019; */ /* daa[11*20 + 1] = 5.940478; */ /* daa[11*20 + 2] = 2.352253; */ /* daa[11*20 + 3] = 0.231001; */ /* daa[11*20 + 4] = 0.027995; */ /* daa[11*20 + 5] = 3.646743; */ /* daa[11*20 + 6] = 1.507981; */ /* daa[11*20 + 7] = 0.331175; */ /* daa[11*20 + 8] = 0.698362; */ /* daa[11*20 + 9] = 0.140326; */ /* daa[11*20 + 10] = 0.171396; */ /* daa[12*20 + 0] = 0.694226; */ /* daa[12*20 + 1] = 0.419899; */ /* daa[12*20 + 2] = 0.326927; */ /* daa[12*20 + 3] = 0.039488; */ /* daa[12*20 + 4] = 0.844827; */ /* daa[12*20 + 5] = 1.394214; */ /* daa[12*20 + 6] = 0.133235; */ /* daa[12*20 + 7] = 0.085075; */ /* daa[12*20 + 8] = 0.347092; */ /* daa[12*20 + 9] = 4.051255; */ /* daa[12*20 + 10] = 6.650794; */ /* daa[12*20 + 11] = 0.617549; */ /* daa[13*20 + 0] = 0.155206; */ /* daa[13*20 + 1] = 0.057971; */ /* daa[13*20 + 2] = 0.09816; */ /* daa[13*20 + 3] = 0.020441; */ /* daa[13*20 + 4] = 0.904305; */ /* daa[13*20 + 5] = 0.052719; */ /* daa[13*20 + 6] = 0.0219; */ /* daa[13*20 + 7] = 0.046668; */ /* daa[13*20 + 8] = 0.890005; */ /* daa[13*20 + 9] = 0.844963; */ /* daa[13*20 + 10] = 2.348881; */ /* daa[13*20 + 11] = 0.028372; */ /* daa[13*20 + 12] = 1.671635; */ /* daa[14*20 + 0] = 1.433475; */ /* daa[14*20 + 1] = 0.328393; */ /* daa[14*20 + 2] = 0.173181; */ /* daa[14*20 + 3] = 0.431874; */ /* daa[14*20 + 4] = 0.09902; */ /* daa[14*20 + 5] = 0.592324; */ /* daa[14*20 + 6] = 0.488352; */ /* daa[14*20 + 7] = 0.23865; */ /* daa[14*20 + 8] = 0.462856; */ /* daa[14*20 + 9] = 0.057048; */ /* daa[14*20 + 10] = 0.233532; */ /* daa[14*20 + 11] = 0.387808; */ /* daa[14*20 + 12] = 0.096377; */ /* daa[14*20 + 13] = 0.079912; */ /* daa[15*20 + 0] = 4.887126; */ /* daa[15*20 + 1] = 0.883923; */ /* daa[15*20 + 2] = 4.627163; */ /* daa[15*20 + 3] = 1.122164; */ /* daa[15*20 + 4] = 3.186667; */ /* daa[15*20 + 5] = 1.085947; */ /* daa[15*20 + 6] = 0.569339; */ /* daa[15*20 + 7] = 1.993432; */ /* daa[15*20 + 8] = 0.867972; */ /* daa[15*20 + 9] = 0.070512; */ /* daa[15*20 + 10] = 0.163009; */ /* daa[15*20 + 11] = 0.718913; */ /* daa[15*20 + 12] = 0.301103; */ /* daa[15*20 + 13] = 0.32579; */ /* daa[15*20 + 14] = 1.449582; */ /* daa[16*20 + 0] = 2.030538; */ /* daa[16*20 + 1] = 0.639463; */ /* daa[16*20 + 2] = 2.076294; */ /* daa[16*20 + 3] = 0.377239; */ /* daa[16*20 + 4] = 1.42848; */ /* daa[16*20 + 5] = 0.979403; */ /* daa[16*20 + 6] = 0.647562; */ /* daa[16*20 + 7] = 0.145556; */ /* daa[16*20 + 8] = 0.493329; */ /* daa[16*20 + 9] = 0.973405; */ /* daa[16*20 + 10] = 0.271824; */ /* daa[16*20 + 11] = 1.20033; */ /* daa[16*20 + 12] = 1.659187; */ /* daa[16*20 + 13] = 0.1217; */ /* daa[16*20 + 14] = 0.571399; */ /* daa[16*20 + 15] = 6.641034; */ /* daa[17*20 + 0] = 0.131405; */ /* daa[17*20 + 1] = 0.552911; */ /* daa[17*20 + 2] = 0.079985; */ /* daa[17*20 + 3] = 0.060514; */ /* daa[17*20 + 4] = 0.633662; */ /* daa[17*20 + 5] = 0.21823; */ /* daa[17*20 + 6] = 0.074988; */ /* daa[17*20 + 7] = 0.169114; */ /* daa[17*20 + 8] = 0.847725; */ /* daa[17*20 + 9] = 0.10627; */ /* daa[17*20 + 10] = 0.622044; */ /* daa[17*20 + 11] = 0.060755; */ /* daa[17*20 + 12] = 0.719575; */ /* daa[17*20 + 13] = 3.14824; */ /* daa[17*20 + 14] = 0.077123; */ /* daa[17*20 + 15] = 0.276716; */ /* daa[17*20 + 16] = 0.148883; */ /* daa[18*20 + 0] = 0.165179; */ /* daa[18*20 + 1] = 0.224883; */ /* daa[18*20 + 2] = 0.528334; */ /* daa[18*20 + 3] = 0.121252; */ /* daa[18*20 + 4] = 1.174118; */ /* daa[18*20 + 5] = 0.177062; */ /* daa[18*20 + 6] = 0.074715; */ /* daa[18*20 + 7] = 0.042356; */ /* daa[18*20 + 8] = 5.911187; */ /* daa[18*20 + 9] = 0.192481; */ /* daa[18*20 + 10] = 0.321454; */ /* daa[18*20 + 11] = 0.090556; */ /* daa[18*20 + 12] = 0.406415; */ /* daa[18*20 + 13] = 10.908861; */ /* daa[18*20 + 14] = 0.070752; */ /* daa[18*20 + 15] = 0.328483; */ /* daa[18*20 + 16] = 0.181539; */ /* daa[18*20 + 17] = 3.823886; */ /* daa[19*20 + 0] = 1.78517; */ /* daa[19*20 + 1] = 0.166975; */ /* daa[19*20 + 2] = 0.106482; */ /* daa[19*20 + 3] = 0.041707; */ /* daa[19*20 + 4] = 1.876812; */ /* daa[19*20 + 5] = 0.22421; */ /* daa[19*20 + 6] = 0.247356; */ /* daa[19*20 + 7] = 0.06688; */ /* daa[19*20 + 8] = 0.1436; */ /* daa[19*20 + 9] = 9.60184; */ /* daa[19*20 + 10] = 1.599119; */ /* daa[19*20 + 11] = 0.17319; */ /* daa[19*20 + 12] = 1.645134; */ /* daa[19*20 + 13] = 0.438571; */ /* daa[19*20 + 14] = 0.252486; */ /* daa[19*20 + 15] = 0.105536; */ /* daa[19*20 + 16] = 1.789097; */ /* daa[19*20 + 17] = 0.147951; */ /* daa[19*20 + 18] = 0.200571; */ /* for (i=0; i 2664 (67% of the original 3933 positions) The species included in the analysis were: Harpiosquilla harpax [NCBI_TaxID 287944] Ixodes uriae [NCBI_TaxID 59655] Heptathela hangzhouensis [NCBI_TaxID 216259] Triops longicaudatus [NCBI_TaxID 58777] Gryllotalpa orientalis [NCBI_TaxID 213494] lepidopsocid RS-2001 [NCBI_TaxID 159971] Locusta migratoria [NCBI_TaxID 7004] Drosophila yakuba [NCBI_TaxID 7245] Ostrinia furnacalis [NCBI_TaxID 93504] Megabalanus volcano [NCBI_TaxID 266495] Periplaneta fuliginosa [NCBI_TaxID 36977] Thermobia domestica [NCBI_TaxID 89055] Aleurochiton aceris [NCBI_TaxID 266942] Schizaphis graminum [NCBI_TaxID 13262] Pteronarcys princeps [NCBI_TaxID 285953] Aleurodicus dugesii [NCBI_TaxID 30099] Pollicipes polymerus [NCBI_TaxID 36137] Gomphiocephalus hodgsoni [NCBI_TaxID 221270] Habronattus oregonensis [NCBI_TaxID 130930] Speleonectes tulumensis [NCBI_TaxID 84346] Hutchinsoniella macracantha [NCBI_TaxID 84335] Haemaphysalis flava [NCBI_TaxID 181088] Scutigera coleoptrata [NCBI_TaxID 29022] Vargula hilgendorfii [NCBI_TaxID 6674] Tricholepidion gertschi [NCBI_TaxID 89825] Varroa destructor [NCBI_TaxID 109461] Bombyx mandarina [NCBI_TaxID 7092] Thyropygus sp. [NCBI_TaxID 174155] Tribolium castaneum [NCBI_TaxID 7070] Pagurus longicarpus [NCBI_TaxID 111067] Limulus polyphemus [NCBI_TaxID 6850] Tetrodontophora bielanensis [NCBI_TaxID 48717] Penaeus monodon [NCBI_TaxID 6687] Daphnia pulex [NCBI_TaxID 6669] Apis mellifera [NCBI_TaxID 7469] Anopheles gambiae [NCBI_TaxID 7165] The topology used for inferring the model was: (((Daph_pulex,Trio_longi),((((((Aleu_aceri,Aleu_duges),Schi_grami),lepi_RS_20), ((((Ostr_furna,Bomb_manda),(Dros_yakub,Anop_gambi)),Apis_melli),Trib_casta)), ((Gryl_orien,Locu_migra),(Pter_princ,Peri_fulig))),(Tric_gerts,Ther_domes)), (Scut_coleo,Thyr_sp),Varg_hilge,Hutc_macra,((((Ixod_uriae,Haem_flava),Varr_destr), (Habr_orego,Hept_hangz)),Limu_polyp),(Poll_polym,Mega_volca),(Gomp_hodgs,Tetr_biela), ((Pagu_longi,Pena_monod),Harp_harpa),Spel_tulum)); Note this is not the ML topoLOGy but the consensus one (based on morphoLOGical data, phyLOGenetic reconstruction using nuclear genes, etc). Where relationships are not clear, a polytomy was introduced (it contains quite a lot of polytomies!). The model was estimated using (the great and helpful) Ziheng Yang's Paml software package. A four-categorized gamma distribution was used to account for heterogeneity (alpha was estimated to be 0.47821). Sites with ambiguity data were taken into account. If you would like the data related to this matrix, please contact fabascal@uvigo.es. Federico Abascal (c)2005. */ int i,j,naa; naa = 20; daa[1*20+ 0] = 0.2; daa[2*20+ 0] = 0.2; daa[2*20+ 1] = 0.2; daa[3*20+ 0] = 0.6; daa[3*20+ 1] = 4.3; daa[3*20+ 2] = 500.2; daa[4*20+ 0] = 253.5; daa[4*20+ 1] = 35.5; daa[4*20+ 2] = 98.2; daa[4*20+ 3] = 10.6; daa[5*20+ 0] = 0.2; daa[5*20+ 1] = 154.0; daa[5*20+ 2] = 261.8; daa[5*20+ 3] = 0.2; daa[5*20+ 4] = 0.2; daa[6*20+ 0] = 0.2; daa[6*20+ 1] = 0.2; daa[6*20+ 2] = 183.0; daa[6*20+ 3] = 861.8; daa[6*20+ 4] = 0.2; daa[6*20+ 5] = 261.6; daa[7*20+ 0] = 199.8; daa[7*20+ 1] = 0.2; daa[7*20+ 2] = 120.5; daa[7*20+ 3] = 12.5; daa[7*20+ 4] = 80.5; daa[7*20+ 5] = 2.6; daa[7*20+ 6] = 43.9; daa[8*20+ 0] = 0.2; daa[8*20+ 1] = 41.3; daa[8*20+ 2] = 179.5; daa[8*20+ 3] = 0.2; daa[8*20+ 4] = 12.4; daa[8*20+ 5] = 313.5; daa[8*20+ 6] = 15.2; daa[8*20+ 7] = 0.2; daa[9*20+ 0] = 25.7; daa[9*20+ 1] = 1.8; daa[9*20+ 2] = 21.3; daa[9*20+ 3] = 6.6; daa[9*20+ 4] = 63.0; daa[9*20+ 5] = 10.5; daa[9*20+ 6] = 6.8; daa[9*20+ 7] = 2.7; daa[9*20+ 8] = 0.2; daa[10*20+ 0] = 3.7; daa[10*20+ 1] = 1.8; daa[10*20+ 2] = 12.6; daa[10*20+ 3] = 1.2; daa[10*20+ 4] = 78.7; daa[10*20+ 5] = 16.3; daa[10*20+ 6] = 1.7; daa[10*20+ 7] = 1.4; daa[10*20+ 8] = 5.5; daa[10*20+ 9] = 514.5; daa[11*20+ 0] = 0.2; daa[11*20+ 1] = 208.6; daa[11*20+ 2] = 467.3; daa[11*20+ 3] = 1.7; daa[11*20+ 4] = 0.2; daa[11*20+ 5] = 349.3; daa[11*20+ 6] = 106.3; daa[11*20+ 7] = 0.2; daa[11*20+ 8] = 0.2; daa[11*20+ 9] = 3.5; daa[11*20+ 10] = 3.8; daa[12*20+ 0] = 120.6; daa[12*20+ 1] = 5.2; daa[12*20+ 2] = 78.8; daa[12*20+ 3] = 0.2; daa[12*20+ 4] = 312.3; daa[12*20+ 5] = 67.3; daa[12*20+ 6] = 0.2; daa[12*20+ 7] = 55.7; daa[12*20+ 8] = 0.2; daa[12*20+ 9] = 514.8; daa[12*20+ 10] = 885.5; daa[12*20+ 11] = 105.6; daa[13*20+ 0] = 13.1; daa[13*20+ 1] = 4.7; daa[13*20+ 2] = 19.7; daa[13*20+ 3] = 0.2; daa[13*20+ 4] = 184.1; daa[13*20+ 5] = 0.2; daa[13*20+ 6] = 0.2; daa[13*20+ 7] = 0.8; daa[13*20+ 8] = 13.8; daa[13*20+ 9] = 117.9; daa[13*20+ 10] = 262.6; daa[13*20+ 11] = 10.7; daa[13*20+ 12] = 321.6; daa[14*20+ 0] = 49.3; daa[14*20+ 1] = 0.2; daa[14*20+ 2] = 16.5; daa[14*20+ 3] = 0.2; daa[14*20+ 4] = 0.2; daa[14*20+ 5] = 39.3; daa[14*20+ 6] = 7.9; daa[14*20+ 7] = 0.2; daa[14*20+ 8] = 0.8; daa[14*20+ 9] = 0.2; daa[14*20+ 10] = 12.2; daa[14*20+ 11] = 16.8; daa[14*20+ 12] = 5.3; daa[14*20+ 13] = 14.6; daa[15*20+ 0] = 673.0; daa[15*20+ 1] = 2.7; daa[15*20+ 2] = 398.4; daa[15*20+ 3] = 44.4; daa[15*20+ 4] = 664.2; daa[15*20+ 5] = 52.4; daa[15*20+ 6] = 31.5; daa[15*20+ 7] = 226.0; daa[15*20+ 8] = 10.6; daa[15*20+ 9] = 7.2; daa[15*20+ 10] = 8.2; daa[15*20+ 11] = 144.2; daa[15*20+ 12] = 111.7; daa[15*20+ 13] = 36.1; daa[15*20+ 14] = 86.5; daa[16*20+ 0] = 243.9; daa[16*20+ 1] = 0.2; daa[16*20+ 2] = 165.9; daa[16*20+ 3] = 0.2; daa[16*20+ 4] = 182.8; daa[16*20+ 5] = 43.7; daa[16*20+ 6] = 43.4; daa[16*20+ 7] = 0.2; daa[16*20+ 8] = 18.6; daa[16*20+ 9] = 203.7; daa[16*20+ 10] = 47.8; daa[16*20+ 11] = 69.5; daa[16*20+ 12] = 288.6; daa[16*20+ 13] = 13.5; daa[16*20+ 14] = 46.8; daa[16*20+ 15] = 660.4; daa[17*20+ 0] = 0.2; daa[17*20+ 1] = 0.2; daa[17*20+ 2] = 7.7; daa[17*20+ 3] = 0.2; daa[17*20+ 4] = 21.6; daa[17*20+ 5] = 6.7; daa[17*20+ 6] = 11.0; daa[17*20+ 7] = 1.9; daa[17*20+ 8] = 0.2; daa[17*20+ 9] = 0.2; daa[17*20+ 10] = 21.1; daa[17*20+ 11] = 16.0; daa[17*20+ 12] = 70.7; daa[17*20+ 13] = 53.7; daa[17*20+ 14] = 0.2; daa[17*20+ 15] = 2.4; daa[17*20+ 16] = 0.2; daa[18*20+ 0] = 1.2; daa[18*20+ 1] = 3.9; daa[18*20+ 2] = 251.2; daa[18*20+ 3] = 0.2; daa[18*20+ 4] = 72.0; daa[18*20+ 5] = 86.7; daa[18*20+ 6] = 7.7; daa[18*20+ 7] = 8.6; daa[18*20+ 8] = 191.4; daa[18*20+ 9] = 12.3; daa[18*20+ 10] = 19.8; daa[18*20+ 11] = 117.1; daa[18*20+ 12] = 70.9; daa[18*20+ 13] = 791.6; daa[18*20+ 14] = 18.4; daa[18*20+ 15] = 30.5; daa[18*20+ 16] = 46.0; daa[18*20+ 17] = 37.7; daa[19*20+ 0] = 339.9; daa[19*20+ 1] = 0.2; daa[19*20+ 2] = 22.6; daa[19*20+ 3] = 0.2; daa[19*20+ 4] = 350.4; daa[19*20+ 5] = 0.2; daa[19*20+ 6] = 13.6; daa[19*20+ 7] = 2.6; daa[19*20+ 8] = 0.2; daa[19*20+ 9] = 1854.5; daa[19*20+ 10] = 84.7; daa[19*20+ 11] = 26.1; daa[19*20+ 12] = 281.3; daa[19*20+ 13] = 51.9; daa[19*20+ 14] = 31.7; daa[19*20+ 15] = 60.6; daa[19*20+ 16] = 544.1; daa[19*20+ 17] = 0.2; daa[19*20+ 18] = 1.6; /* MtArt.old: esta es la MtArt que hice con 26 secuencias (2 outgroups) con una topoLOGia incorrecta daa[1*20+ 0] = 0.2; daa[2*20+ 0] = 0.2; daa[2*20+ 1] = 8.0; daa[3*20+ 0] = 0.2; daa[3*20+ 1] = 0.2; daa[3*20+ 2] = 441.7; daa[4*20+ 0] = 287.9; daa[4*20+ 1] = 48.4; daa[4*20+ 2] = 82.4; daa[4*20+ 3] = 0.2; daa[5*20+ 0] = 0.2; daa[5*20+ 1] = 149.9; daa[5*20+ 2] = 278.6; daa[5*20+ 3] = 0.2; daa[5*20+ 4] = 21.7; daa[6*20+ 0] = 6.6; daa[6*20+ 1] = 0.2; daa[6*20+ 2] = 213.9; daa[6*20+ 3] = 760.8; daa[6*20+ 4] = 0.2; daa[6*20+ 5] = 292.9; daa[7*20+ 0] = 228.2; daa[7*20+ 1] = 0.2; daa[7*20+ 2] = 97.1; daa[7*20+ 3] = 10.4; daa[7*20+ 4] = 98.4; daa[7*20+ 5] = 4.0; daa[7*20+ 6] = 48.7; daa[8*20+ 0] = 0.2; daa[8*20+ 1] = 56.7; daa[8*20+ 2] = 156.4; daa[8*20+ 3] = 24.5; daa[8*20+ 4] = 15.5; daa[8*20+ 5] = 328.6; daa[8*20+ 6] = 7.0; daa[8*20+ 7] = 8.4; daa[9*20+ 0] = 26.4; daa[9*20+ 1] = 1.6; daa[9*20+ 2] = 40.1; daa[9*20+ 3] = 0.2; daa[9*20+ 4] = 22.1; daa[9*20+ 5] = 13.8; daa[9*20+ 6] = 0.2; daa[9*20+ 7] = 3.6; daa[9*20+ 8] = 0.2; daa[10*20+ 0] = 3.4; daa[10*20+ 1] = 0.6; daa[10*20+ 2] = 13.8; daa[10*20+ 3] = 0.7; daa[10*20+ 4] = 76.9; daa[10*20+ 5] = 12.1; daa[10*20+ 6] = 5.4; daa[10*20+ 7] = 2.5; daa[10*20+ 8] = 2.9; daa[10*20+ 9] = 542.6; daa[11*20+ 0] = 0.2; daa[11*20+ 1] = 240.2; daa[11*20+ 2] = 602.8; daa[11*20+ 3] = 35.5; daa[11*20+ 4] = 0.2; daa[11*20+ 5] = 357.6; daa[11*20+ 6] = 62.6; daa[11*20+ 7] = 0.2; daa[11*20+ 8] = 3.3; daa[11*20+ 9] = 0.2; daa[11*20+ 10] = 17.5; daa[12*20+ 0] = 119.0; daa[12*20+ 1] = 0.2; daa[12*20+ 2] = 91.4; daa[12*20+ 3] = 6.4; daa[12*20+ 4] = 332.3; daa[12*20+ 5] = 65.4; daa[12*20+ 6] = 0.2; daa[12*20+ 7] = 60.4; daa[12*20+ 8] = 2.4; daa[12*20+ 9] = 492.5; daa[12*20+ 10] = 815.8; daa[12*20+ 11] = 67.3; daa[13*20+ 0] = 8.2; daa[13*20+ 1] = 6.4; daa[13*20+ 2] = 31.5; daa[13*20+ 3] = 3.4; daa[13*20+ 4] = 174.4; daa[13*20+ 5] = 5.7; daa[13*20+ 6] = 5.7; daa[13*20+ 7] = 2.1; daa[13*20+ 8] = 11.0; daa[13*20+ 9] = 94.4; daa[13*20+ 10] = 243.3; daa[13*20+ 11] = 12.3; daa[13*20+ 12] = 357.8; daa[14*20+ 0] = 62.5; daa[14*20+ 1] = 0.4; daa[14*20+ 2] = 17.5; daa[14*20+ 3] = 0.2; daa[14*20+ 4] = 0.2; daa[14*20+ 5] = 48.6; daa[14*20+ 6] = 17.7; daa[14*20+ 7] = 2.7; daa[14*20+ 8] = 0.2; daa[14*20+ 9] = 0.2; daa[14*20+ 10] = 11.2; daa[14*20+ 11] = 21.7; daa[14*20+ 12] = 5.2; daa[14*20+ 13] = 12.6; daa[15*20+ 0] = 659.0; daa[15*20+ 1] = 5.2; daa[15*20+ 2] = 469.8; daa[15*20+ 3] = 52.3; daa[15*20+ 4] = 570.7; daa[15*20+ 5] = 47.8; daa[15*20+ 6] = 37.3; daa[15*20+ 7] = 227.8; daa[15*20+ 8] = 12.7; daa[15*20+ 9] = 12.3; daa[15*20+ 10] = 7.4; daa[15*20+ 11] = 189.0; daa[15*20+ 12] = 155.3; daa[15*20+ 13] = 43.8; daa[15*20+ 14] = 103.4; daa[16*20+ 0] = 276.4; daa[16*20+ 1] = 1.6; daa[16*20+ 2] = 175.6; daa[16*20+ 3] = 0.2; daa[16*20+ 4] = 96.2; daa[16*20+ 5] = 71.4; daa[16*20+ 6] = 37.4; daa[16*20+ 7] = 0.2; daa[16*20+ 8] = 14.2; daa[16*20+ 9] = 212.5; daa[16*20+ 10] = 38.5; daa[16*20+ 11] = 97.4; daa[16*20+ 12] = 254.7; daa[16*20+ 13] = 2.1; daa[16*20+ 14] = 41.6; daa[16*20+ 15] = 670.6; daa[17*20+ 0] = 6.2; daa[17*20+ 1] = 0.2; daa[17*20+ 2] = 0.2; daa[17*20+ 3] = 5.6; daa[17*20+ 4] = 0.2; daa[17*20+ 5] = 0.2; daa[17*20+ 6] = 3.1; daa[17*20+ 7] = 0.4; daa[17*20+ 8] = 0.2; daa[17*20+ 9] = 15.2; daa[17*20+ 10] = 11.5; daa[17*20+ 11] = 32.6; daa[17*20+ 12] = 82.4; daa[17*20+ 13] = 81.9; daa[17*20+ 14] = 0.2; daa[17*20+ 15] = 9.7; daa[17*20+ 16] = 0.2; daa[18*20+ 0] = 1.6; daa[18*20+ 1] = 7.7; daa[18*20+ 2] = 242.5; daa[18*20+ 3] = 0.2; daa[18*20+ 4] = 88.0; daa[18*20+ 5] = 93.1; daa[18*20+ 6] = 0.2; daa[18*20+ 7] = 6.0; daa[18*20+ 8] = 113.7; daa[18*20+ 9] = 22.1; daa[18*20+ 10] = 17.2; daa[18*20+ 11] = 138.5; daa[18*20+ 12] = 37.6; daa[18*20+ 13] = 770.2; daa[18*20+ 14] = 5.3; daa[18*20+ 15] = 25.0; daa[18*20+ 16] = 55.5; daa[18*20+ 17] = 69.3; daa[19*20+ 0] = 307.8; daa[19*20+ 1] = 2.2; daa[19*20+ 2] = 6.9; daa[19*20+ 3] = 0.2; daa[19*20+ 4] = 405.7; daa[19*20+ 5] = 0.8; daa[19*20+ 6] = 20.2; daa[19*20+ 7] = 5.7; daa[19*20+ 8] = 0.2; daa[19*20+ 9] = 1687.9; daa[19*20+ 10] = 49.4; daa[19*20+ 11] = 23.4; daa[19*20+ 12] = 329.9; daa[19*20+ 13] = 86.3; daa[19*20+ 14] = 27.3; daa[19*20+ 15] = 95.0; daa[19*20+ 16] = 443.0; daa[19*20+ 17] = 2.4; daa[19*20+ 18] = 0.2; 3*/ for (i=0; iio->datatype == NT) m4mod->n_o = 4; else if(mod->io->datatype == AA) m4mod->n_o = 20; else { PhyML_Printf("\n== Not implemented yet."); PhyML_Printf("\n== Err in file %s at line %d\n",__FILE__,__LINE__); Warn_And_Exit(""); } mod->ns = m4mod->n_o * m4mod->n_h; For(i,m4mod->n_o) m4mod->o_fq[i] = mod->e_frq->pi->v[i]; /*! At that stage, the mod->pi vector as been initialized under a standard non covarion type of model. Use these frequencies as they have been set according to the nucleotide substitution model chosen (e.g., 1/4 for JC69). !*/ For(i,(int)(m4mod->n_h)) m4mod->multipl[i] = 1.; ct = 0; For(i,m4mod->n_o-1) { for(j=i+1;jn_o;j++) { m4mod->o_rr[ct] = MAX(mod->r_mat->qmat->v[i*m4mod->n_o+j],1.E-5); ct++; } } For(i,(int)(m4mod->n_h*(m4mod->n_h-1)/2)) m4mod->h_rr[i] = 1.; fq = (phydbl)(1./m4mod->n_h); if(mod->s_opt->opt_cov_delta) m4mod->delta = 1.0; if(mod->s_opt->opt_cov_alpha) m4mod->alpha = 1.0; For(i,m4mod->n_h) m4mod->h_fq[i] = fq; For(i,m4mod->n_h) m4mod->h_fq_unscaled[i] = 1.0; For(i,m4mod->n_h) m4mod->multipl[i] = (phydbl)i; For(i,m4mod->n_h) m4mod->multipl_unscaled[i] = (phydbl)i; Switch_Eigen(YES,mod); M4_Update_Qmat(m4mod,mod); } ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////