Skip to content

Commit 31f56ce

Browse files
4ment4ment
committed
Fixed some bugs
1 parent 86ac5b6 commit 31f56ce

File tree

6 files changed

+54
-42
lines changed

6 files changed

+54
-42
lines changed

‎src/phyc/heterotachy.c

+1-2
Original file line numberDiff line numberDiff line change
@@ -113,8 +113,7 @@ ClockSearch * new_ClockSearch( SingleTreeLikelihood *tlk, unsigned nThreads ){
113113

114114
search->ic = INFORMATION_CRITERION_AICc;
115115
search->IC = ClockSearch_IC;
116-
//search->ic_sample_size = search->pool->tlks[0]->sp->count;
117-
search->ic_sample_size = search->pool->tlks[0]->sp->nsites;
116+
search->ic_sample_size = search->pool->tlks[0]->sp->count;
118117

119118
return search;
120119
}

‎src/phyc/optimize.c

+1-1
Original file line numberDiff line numberDiff line change
@@ -1542,7 +1542,7 @@ double optimize_singletreelikelihood( SingleTreeLikelihood *stlk ){
15421542
// Scaling of root height and its children if they have short branches and rate at the same time
15431543
// Works very well when a child of the root is stuck to the root and the root cannot descrease its height
15441544

1545-
if( !isCalibrated ){
1545+
if( !isCalibrated && false){
15461546
if( !Node_isleaf( Node_left(Tree_root(stlk->tree)) ) && !Node_isleaf( Node_right(Tree_root(stlk->tree)) ) ){
15471547
opt_set_objective_function(opt_scale_height_rate, _brent_optimize_scale_root_height_rate);
15481548

‎src/phyc/patristic.c

+1-1
Original file line numberDiff line numberDiff line change
@@ -55,7 +55,7 @@ double ** calculate_patristic( Tree *tree ){
5555
b = Node_parent(b);
5656
}
5757
}
58-
matrix[i][j] = d;
58+
matrix[Node_class_id(nodes[i])][Node_class_id(nodes[j])] = matrix[Node_class_id(nodes[j])][Node_class_id(nodes[i])] = d;
5959
}
6060
}
6161
return matrix;

‎src/phyc/phyboot.c

-1
Original file line numberDiff line numberDiff line change
@@ -393,7 +393,6 @@ void SingleTreeLikelihood_resampling_openmp( const SingleTreeLikelihood *tlk, re
393393

394394
SingleTreeLikelihood *tlk2 = new_SingleTreeLikelihood(tree, sm, sp, bm);
395395
OptConfig_copy(&tlk->opt, &tlk2->opt);
396-
tlk2->opt.verbosity = 1;
397396

398397
if( tlk2->opt.topology_optimize ){
399398
TopologyOptimizer *topology = new_TopologyOptimizer( tlk2, TREE_SEARCH_PARSIMONY_SPR );

‎src/phyc/spropt.c

+41-32
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@
3333
#include <omp.h>
3434
#endif
3535

36-
#define DEBUG_TOPOLOGY_SPR 1
36+
//#define DEBUG_TOPOLOGY_SPR 1
3737

3838

3939
void compare_trees(Node *node1, Node *node2 ){
@@ -251,16 +251,18 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
251251
double lnl = stlk->calculate(stlk);
252252
double lnl2 = lnl+10;
253253
data->index_param = -1;
254+
#ifdef DEBUG_TOPOLOGY_SPR
254255
printf("-- lnl %f\n",lnl);
255-
256+
#endif
256257

257258
if( Node_isroot(node1) ){
258259
printf("Cannot be called on the root\n%s\n%s: %d\n", __func__, __FILE__, __LINE__);
259260
exit(1);
260261
}
261262

262-
263+
#ifdef DEBUG_TOPOLOGY_SPR
263264
printf("%f %f %f %f\n", Node_distance(node1), Node_distance(node1->left), Node_distance(node1->right), Node_distance(node2));
265+
#endif
264266
while ( 1 ) {
265267

266268
SingleTreeLikelihood_update_all_nodes(stlk);
@@ -272,8 +274,9 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
272274
if( status == OPT_ERROR ) error("OPT.DISTANCE No SUCCESS!!!!!!!!!!!!\n");
273275

274276
Node_set_distance( node1->left, Parameters_value(param, 0) );
275-
277+
#ifdef DEBUG_TOPOLOGY_SPR
276278
printf("** lnl %f (%f) %s %f\n", -lnl2, lnl,node1->left->name, Parameters_value(param, 0) );
279+
#endif
277280

278281

279282

@@ -286,9 +289,9 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
286289
if( status == OPT_ERROR ) error("OPT.DISTANCE No SUCCESS!!!!!!!!!!!!\n");
287290

288291
Node_set_distance( node1->right, Parameters_value(param, 0) );
289-
292+
#ifdef DEBUG_TOPOLOGY_SPR
290293
printf("** lnl %f (%f) %s %f\n", -lnl2, lnl,node1->right->name, Parameters_value(param, 0) );
291-
294+
#endif
292295

293296
SingleTreeLikelihood_update_all_nodes(stlk);
294297
data->index_param = node1->postorder_idx;
@@ -299,9 +302,9 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
299302
if( status == OPT_ERROR ) error("OPT.DISTANCE No SUCCESS!!!!!!!!!!!!\n");
300303

301304
Node_set_distance( node1, Parameters_value(param, 0) );
302-
305+
#ifdef DEBUG_TOPOLOGY_SPR
303306
printf("** lnl %f (%f) %s %f\n", -lnl2, lnl,node1->name, Parameters_value(param, 0) );
304-
307+
#endif
305308

306309

307310
SingleTreeLikelihood_update_all_nodes(stlk);
@@ -313,9 +316,9 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
313316
if( status == OPT_ERROR ) error("OPT.DISTANCE No SUCCESS!!!!!!!!!!!!\n");
314317

315318
Node_set_distance( node2, Parameters_value(param, 0) );
316-
319+
#ifdef DEBUG_TOPOLOGY_SPR
317320
printf("** lnl %f (%f) %s %f\n", -lnl2, lnl,node2->name, Parameters_value(param, 0) );
318-
321+
#endif
319322
if ( -lnl2 - lnl < 0.01 ){
320323
lnl = -lnl2;
321324
break;
@@ -324,8 +327,9 @@ double optimize_brent_branch_length_3( SingleTreeLikelihood *stlk, Optimizer *op
324327
}
325328

326329
SingleTreeLikelihood_update_all_nodes(stlk);
327-
330+
#ifdef DEBUG_TOPOLOGY_SPR
328331
printf("%f %f %f %f\n", Node_distance(node1), Node_distance(node1->left), Node_distance(node1->right), Node_distance(node2));
332+
#endif
329333
//return stlk->calculate(stlk);
330334
return lnl;
331335
}
@@ -422,9 +426,9 @@ double spr_optimize_bl_parsimony( struct TopologyOptimizer * opt ){
422426

423427
Parsimony *parsimony = new_Parsimony(tlk2->sp, tlk2->tree);
424428
double score = parsimony->calculate(parsimony);
425-
426-
printf("\nParsimony score: %f\n\n", score);
427-
429+
if(tlk->opt.verbosity > 0){
430+
printf("\nParsimony score: %f\n\n", score);
431+
}
428432

429433
Optimizer *opt_bl = new_Optimizer(OPT_BRENT);
430434
BrentData *data_brent = new_BrentData( tlk2 );
@@ -809,7 +813,6 @@ double spr_optimize_bl_parsimony( struct TopologyOptimizer * opt ){
809813
//#define SPR_THREAD 1
810814
//#define _OPENMP 1
811815

812-
#define DEBUG_TOPOLOGY_SPR 1
813816

814817
// neighborhood size 2(n − 3)(2n − 7)
815818
double spr_optimize_bl_parsimony_only( struct TopologyOptimizer * opt ){
@@ -826,9 +829,10 @@ double spr_optimize_bl_parsimony_only( struct TopologyOptimizer * opt ){
826829

827830
int max_distance = 20;//5;
828831

829-
printf("\nParsimony score: %f (%i)\n\n", score, total);
830-
printf("Maximum radius %d\n\n", max_distance);
831-
832+
if(tlk->opt.verbosity > 0){
833+
printf("\nParsimony score: %f (%i)\n\n", score, total);
834+
printf("Maximum radius %d\n\n", max_distance);
835+
}
832836

833837
int nNodes = Tree_node_count(tlk->tree);
834838
bool failed = false;
@@ -1110,9 +1114,9 @@ double spr_optimize_bl_parsimony_only( struct TopologyOptimizer * opt ){
11101114
Parsimony_update_all_nodes(parsimony);
11111115

11121116
spr_score = parsimony->calculate(parsimony);
1113-
1114-
printf("Parsimony %f (%f) d(prune,graft)=%d max %f\n", spr_score, scores[i], ds[i], score);
1115-
1117+
if(tlk->opt.verbosity > 0){
1118+
printf("Parsimony %f (%f) d(prune,graft)=%d max %f\n", spr_score, scores[i], ds[i], score);
1119+
}
11161120
// skip moves that increase the score
11171121
if ( spr_score >= score ) {
11181122

@@ -1139,7 +1143,9 @@ double spr_optimize_bl_parsimony_only( struct TopologyOptimizer * opt ){
11391143
opt->moves++;
11401144
}
11411145
}
1142-
printf("Parsimony score: %f #SPR: %d\n\n", score, accepted);
1146+
if(tlk->opt.verbosity > 0){
1147+
printf("Parsimony score: %f #SPR: %d\n\n", score, accepted);
1148+
}
11431149
}
11441150
else {
11451151
failed = true;
@@ -1183,8 +1189,9 @@ double spr_optimize_bl_parsimony_only_openmp( struct TopologyOptimizer * opt ){
11831189

11841190
int total = (Tree_tip_count(tlk->tree)-3)*2*(2*Tree_tip_count(tlk->tree)-7);
11851191

1186-
printf("\nParsimony score: %f (%i)\n\n", score, total);
1187-
1192+
if(tlk->opt.verbosity > 0){
1193+
printf("\nParsimony score: %f (%i)\n\n", score, total);
1194+
}
11881195

11891196
int nNodes = Tree_node_count(tlk->tree);
11901197
bool failed = false;
@@ -1510,7 +1517,7 @@ double spr_optimize_bl_parsimony_only_openmp( struct TopologyOptimizer * opt ){
15101517
//break;
15111518
}
15121519
else {
1513-
#ifndef DEBUG_TOPOLOGY_SPR
1520+
#ifdef DEBUG_TOPOLOGY_SPR
15141521
printf("Parsimony %f (%f) d(prune,graft)=%d max %f\n", spr_score, scores[i], ds[i], score);
15151522
#endif
15161523

@@ -1526,7 +1533,9 @@ double spr_optimize_bl_parsimony_only_openmp( struct TopologyOptimizer * opt ){
15261533
opt->moves++;
15271534
}
15281535
}
1529-
printf("Parsimony score: %f #SPR: %d\n\n", score, accepted);
1536+
if(tlk->opt.verbosity > 0){
1537+
printf("Parsimony score: %f #SPR: %d\n\n", score, accepted);
1538+
}
15301539
}
15311540
else {
15321541
failed = true;
@@ -1789,9 +1798,9 @@ double spr_optimize_bl_openmp( struct TopologyOptimizer * opt ){
17891798
free_Parameters_soft(oneparameter);
17901799

17911800

1792-
//if(spr_score > lnl){
1801+
#ifdef DEBUG_TOPOLOGY_SPR
17931802
printf("%s %s lnl %f (%f) %s\n\n",graft->name, prune->name, spr_score, lnl, (spr_score > lnl ? "*" : ""));
1794-
//}
1803+
#endif
17951804

17961805
//spr_score = tlk->calculate(tlk);
17971806

@@ -1914,9 +1923,9 @@ double spr_optimize_bl_openmp( struct TopologyOptimizer * opt ){
19141923
SingleTreeLikelihood_update_all_nodes(tlk);
19151924

19161925
spr_score = tlk->calculate(tlk);
1917-
1926+
#ifdef DEBUG_TOPOLOGY_SPR
19181927
printf("Parsimony %f (%f) d(prune,graft)=%d max %f\n", spr_score, scores[i], ds[i], lnl);
1919-
1928+
#endif
19201929
// skip moves that increase the score
19211930
if ( spr_score >= lnl ) {
19221931

@@ -1937,7 +1946,7 @@ double spr_optimize_bl_openmp( struct TopologyOptimizer * opt ){
19371946
//break;
19381947
}
19391948
else {
1940-
#ifndef DEBUG_TOPOLOGY_SPR
1949+
#ifdef DEBUG_TOPOLOGY_SPR
19411950
printf("Parsimony %f (%f) d(prune,graft)=%d max %f\n", spr_score, scores[i], ds[i], score);
19421951
#endif
19431952
SingleTreeLikelihood_update_all_nodes(pool[0]);
@@ -1949,7 +1958,7 @@ double spr_optimize_bl_openmp( struct TopologyOptimizer * opt ){
19491958
lnl = spr_score;
19501959
}
19511960

1952-
#ifndef DEBUG_TOPOLOGY_SPR
1961+
#ifdef DEBUG_TOPOLOGY_SPR
19531962
printf("Parsimony %f\n", score);
19541963
#endif
19551964
opt->moves++;

‎src/physher.c

+10-5
Original file line numberDiff line numberDiff line change
@@ -320,10 +320,11 @@ void run_bootstrap( SingleTreeLikelihood *tlk, const char *output_stem, int nthr
320320
Node **nodes = NULL;
321321
Tree *tree = tlk->tree;
322322
double ci_alpha = 0.05; // should be an option
323-
bool save_patterns = true; // should be an option
323+
bool save_patterns = false; // should be an option
324324
char *jackfile_trees = NULL;
325325
char *jackfile_params = NULL;
326326

327+
int verbosity = tlk->opt.verbosity;
327328
tlk->opt.verbosity = 0;
328329

329330
if( tlk->bm != NULL ){
@@ -428,7 +429,7 @@ void run_bootstrap( SingleTreeLikelihood *tlk, const char *output_stem, int nthr
428429
if( file_exists(buffer->c) ){
429430
Phyboot_confidence_intervals(buffer->c, 1-ci_alpha, jackfile_params);
430431
}
431-
432+
tlk->opt.verbosity = verbosity;
432433

433434
free_StringBuffer(buffer);
434435
if(jackfile_params != NULL )free(jackfile_params);
@@ -569,6 +570,10 @@ void run_local_greedy_likelihood( SingleTreeLikelihood *tlk, const char* output_
569570

570571
if ( sample_size > 0 ) {
571572
ClockSearch_set_ic_sample_size(greedy_local, sample_size );
573+
fprintf(stdout, "Sample size: %d (custom)\n", greedy_local->ic_sample_size);
574+
}
575+
else{
576+
fprintf(stdout, "Sample size: %d (number of site patterns)\n", greedy_local->ic_sample_size);
572577
}
573578

574579

@@ -1459,7 +1464,7 @@ int main(int argc, char* argv[]){
14591464
{ARGS_OPTION_FLAG, 'U', "unrooted", "tree.unrooted", &tree_unrooted, "Input tree is rooted"},
14601465
{ARGS_OPTION_DOUBLE, 's', "scaler", "tree.scaler", &tree_scaler, "Scale input tree"},
14611466
{ARGS_OPTION_STRING, 'O', "treeopt", "tree.topolology.optimize", &topology_optimization_algorithm, "Optimize topology nni or spr (experimental)"},
1462-
{ARGS_OPTION_FLAG, 0, "parsimony", "tree.topolology.optimize.parsimony", &use_parsimony, "Quick optimizaton of tree topology using parsimony before ML optimization"},
1467+
{ARGS_OPTION_BOOLEAN, 0, "parsimony", "tree.topolology.optimize.parsimony", &use_parsimony, "Quick optimizaton of tree topology using parsimony before ML optimization"},
14631468

14641469
{ARGS_OPTION_STRING, 'C', "clock", "clock", &clock, "Clock type: strict, local, discrete"},
14651470
{ARGS_OPTION_FLAG, 0, "forward", "clock.forward", &forward, "Time is forward"},
@@ -2280,10 +2285,10 @@ int main(int argc, char* argv[]){
22802285
tlk->opt.bl.tolx = 0.00000001;
22812286
//tlk->opt.bl.tolx = 0.0001;
22822287
//tlk->opt.bl.method = OPT_CG_PR;
2283-
tlk->opt.freqs.optimize = !frequencies_fixed;
2288+
tlk->opt.freqs.optimize = !frequencies_fixed && sm->m->freqs != NULL;
22842289
//tlk->opt.freqs.method = OPT_CG_PR;
22852290
//if( Parameters_count(mod->rates) > 1 ) tlk->opt.relative_rates.method = OPT_CG_PR; //Brent optmize better in the flub simulated data
2286-
tlk->opt.relative_rates.optimize = !rates_fixed;
2291+
tlk->opt.relative_rates.optimize = !rates_fixed && sm->m->rates != NULL;
22872292
tlk->opt.pinv.optimize = use_pinv && !pinv_fixed;
22882293
tlk->opt.gamma.optimize = use_gamma && !alpha_fixed;
22892294
tlk->opt.topology_optimize = false;

0 commit comments

Comments
 (0)