8
#include <aw_device.hxx>
9
#include <aw_window.hxx>
10
#include <aw_preset.hxx>
11
#include <awt_canvas.hxx>
12
#include <awt_tree.hxx>
13
#include <awt_seq_dna.hxx>
14
#include <awt_seq_protein.hxx>
16
#include <awt_csp.hxx>
18
#include <awt_dtree.hxx>
19
#include <awt_sel_boxes.hxx>
1
// =============================================================== //
3
// File : PARS_dtree.cxx //
6
// Institute of Microbiology (Technical University Munich) //
7
// http://www.arb-home.de/ //
9
// =============================================================== //
20
11
#include "pars_dtree.hxx"
22
#include "AP_buffer.hxx"
23
#include "parsimony.hxx"
24
#include "ap_tree_nlen.hxx"
25
12
#include "pars_main.hxx"
26
13
#include "pars_debug.hxx"
28
extern AWT_csp *awt_csp;
30
char *AWT_graphic_parsimony_root_changed(void *cd, AP_tree *old, AP_tree *newroot)
14
#include "ap_tree_nlen.hxx"
16
#include <AP_seq_dna.hxx>
17
#include <AP_seq_protein.hxx>
18
#include <AP_filter.hxx>
20
#include <ColumnStat.hxx>
21
#include <awt_sel_boxes.hxx>
22
#include <awt_filter.hxx>
24
#include <gui_aliview.hxx>
26
#include <aw_preset.hxx>
27
#include <aw_awar.hxx>
29
#include <arb_progress.h>
30
#include <aw_root.hxx>
31
#include <aw_question.hxx>
33
static void AWT_graphic_parsimony_root_changed(void *cd, AP_tree *old, AP_tree *newroot) {
32
34
AWT_graphic_tree *agt = (AWT_graphic_tree*)cd;
33
if (old == agt->tree_root_display) agt->tree_root_display = newroot;
34
if (old == agt->tree_root) agt->tree_root = newroot;
35
if (old == GLOBAL_NT->tree->tree_root) GLOBAL_NT->tree->tree_root = newroot;
36
if (old == agt->displayed_root) agt->displayed_root = newroot;
40
/**************************
43
Filter & weights setup
44
loads sequences at the leafs and
45
initalize filters & weights for the tree
46
( AP_tree_nlen expected )
48
**************************/
49
void NT_tree_init(AWT_graphic_tree *agt, adfiltercbstruct *pars_global_filter) {
51
AP_tree *tree = agt->tree_root;
52
GB_transaction dummy(GLOBAL_gb_main);
39
static AliView *pars_generate_aliview(WeightedFilter *pars_weighted_filter) {
40
GBDATA *gb_main = pars_weighted_filter->get_gb_main();
43
GB_transaction ta(gb_main);
44
ali_name = GBT_read_string(gb_main, AWAR_ALIGNMENT);
57
char *use = GBT_read_string(GLOBAL_gb_main,AWAR_ALIGNMENT);
59
long ali_len = GBT_get_alignment_len(GLOBAL_gb_main,use);
46
GB_ERROR error = NULL;
47
AliView *aliview = pars_weighted_filter->create_aliview(ali_name, error);
48
if (!aliview) aw_popup_exit(error);
53
void PARS_tree_init(AWT_graphic_tree *agt) {
54
ap_assert(agt->get_root_node());
55
ap_assert(agt == ap_main->get_tree_root());
57
GB_transaction ta(GLOBAL_gb_main);
59
const char *use = ap_main->get_aliname();
60
long ali_len = GBT_get_alignment_len(GLOBAL_gb_main, use);
61
62
aw_popup_exit("No valid alignment selected! Try again");
65
GB_BOOL is_aa = GBT_is_alignment_protein(GLOBAL_gb_main,use);
67
// filter & weights setup
69
if (!tree->tree_root->sequence_template) {
70
AP_tree_root *tr = tree->tree_root;
73
sproto = (AP_sequence *)new AP_sequence_protein(tr);
75
sproto = (AP_sequence *)new AP_sequence_parsimony(tr);
78
tr->sequence_template = sproto;
79
tr->filter = awt_get_filter(agt->aw_root, pars_global_filter);
80
tr->weights = new AP_weights();
85
for (i=0;i<ali_len;i++){
86
if (awt_csp->rates[i]>0.0000001){
87
awt_csp->weights[i] *= (int)(2.0/ awt_csp->rates[i]);
90
tr->weights->init(awt_csp->weights , tr->filter);
92
tr->weights->init(tr->filter);
94
tree->load_sequences_rek(use,GB_FALSE,GB_TRUE); // load with sequences
96
tree->tree_root->root_changed_cd = (void*)agt;
97
tree->tree_root->root_changed = AWT_graphic_parsimony_root_changed;
65
agt->tree_static->set_root_changed_callback(AWT_graphic_parsimony_root_changed, agt);
102
static int ap_global_abort_flag;
104
double funktion_quadratisch(double x,double *param_list,int param_anz) {
106
double wert = (double)x;
68
static double funktion_quadratisch(double wert, double *param_list, int param_anz) {
107
69
if (param_anz != 3) {
108
AW_ERROR("funktion_quadratisch: Falsche Parameteranzahl !");
70
ap_assert(0); // wrong number of parameters
111
ergebnis = wert * wert * param_list[0] + wert * param_list[1] + param_list[2];
73
return wert * wert * param_list[0] + wert * param_list[1] + param_list[2];
116
void PARS_kernighan_cb(AP_tree *tree) {
77
void ArbParsimony::kernighan_optimize_tree(AP_tree *at) {
118
78
GB_push_transaction(GLOBAL_gb_main);
120
AP_sequence::global_combineCount = 0;
122
AP_FLOAT pars_start, pars_prev;
123
pars_prev = pars_start = GLOBAL_NT->tree->tree_root->costs();
125
int rek_deep_max = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/maxdepth");
127
AP_KL_FLAG funktype = (AP_KL_FLAG)*GBT_read_int(GLOBAL_gb_main,"genetic/kh/function_type");
80
long prevCombineCount = AP_sequence::combine_count();
82
AP_FLOAT pars_start = get_root_node()->costs();
83
const AP_FLOAT pars_org = pars_start;
85
int rek_deep_max = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/maxdepth");
87
AP_KL_FLAG funktype = (AP_KL_FLAG)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/function_type");
130
90
double param_list[3];
131
double f_startx,f_maxy,f_maxx,f_max_deep;
91
double f_startx, f_maxy, f_maxx, f_max_deep;
132
92
f_max_deep = (double)rek_deep_max; // x2
133
f_startx = (double)*GBT_read_int(GLOBAL_gb_main,"genetic/kh/dynamic/start");
134
f_maxy = (double)*GBT_read_int(GLOBAL_gb_main,"genetic/kh/dynamic/maxy");
135
f_maxx = (double)*GBT_read_int(GLOBAL_gb_main,"genetic/kh/dynamic/maxx");
93
f_startx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/start");
94
f_maxy = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxy");
95
f_maxx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxx");
137
double (*funktion)(double wert,double *param_list,int param_anz);
97
double (*funktion)(double wert, double *param_list, int param_anz);
138
98
switch (funktype) {
140
100
case AP_QUADRAT_START:
157
117
AP_KL_FLAG searchflag=(AP_KL_FLAG)0;
158
if (*GBT_read_int(GLOBAL_gb_main,"genetic/kh/dynamic/enable")){
118
if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/enable")) {
159
119
searchflag = AP_DYNAMIK;
161
if (*GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/enable")){
121
if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/enable")) {
162
122
searchflag = (AP_KL_FLAG)(searchflag|AP_STATIC);
165
125
int rek_breite[8];
166
rek_breite[0] = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/depth0");
167
rek_breite[1] = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/depth1");
168
rek_breite[2] = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/depth2");
169
rek_breite[3] = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/depth3");
170
rek_breite[4] = *GBT_read_int(GLOBAL_gb_main,"genetic/kh/static/depth4");
126
rek_breite[0] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth0");
127
rek_breite[1] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth1");
128
rek_breite[2] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth2");
129
rek_breite[3] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth3");
130
rek_breite[4] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth4");
171
131
int rek_breite_anz = 5;
173
int anzahl = (int)(*GBT_read_float(GLOBAL_gb_main,"genetic/kh/nodes")*tree->arb_tree_leafsum2());
175
list = tree->getRandomNodes(anzahl);
180
aw_openstatus("KL Optimizer");
183
sprintf(buffer,"Old Parsimony: %f",pars_start);
186
int abort_flag = AP_FALSE;
133
int anzahl = (int)(*GBT_read_float(GLOBAL_gb_main, "genetic/kh/nodes")*at->count_leafs());
134
AP_tree **list = at->getRandomNodes(anzahl);
136
arb_progress progress(anzahl);
138
progress.subtitle(GBS_global_string("Old Parsimony: %.1f", pars_start));
188
140
GB_pop_transaction(GLOBAL_gb_main);
190
for (i=0;i<anzahl && ! abort_flag; i++) {
191
abort_flag |= aw_status(i/(double)anzahl);
192
if (abort_flag) break;
194
AP_tree_nlen *tree_elem = (AP_tree_nlen *)list[i];
196
if (tree_elem->gr.hidden ||
197
(tree_elem->father && tree_elem->father->gr.hidden)){
198
continue; // within a folded group
201
AP_BOOL better_tree_found = AP_FALSE;
142
for (int i=0; i<anzahl && !progress.aborted(); i++) {
143
AP_tree_nlen *tree_elem = DOWNCAST(AP_tree_nlen*, list[i]); // @@@ pass 'at' as AP_tree_nlen
145
bool in_folded_group = tree_elem->gr.hidden ||
146
(tree_elem->father && tree_elem->get_father()->gr.hidden);
148
if (!in_folded_group) {
149
bool better_tree_found = false;
203
display_clear(funktion,param_list,param_anz,(int)pars_start,(int)rek_deep_max);
151
display_clear(funktion, param_list, param_anz, (int)pars_start, (int)rek_deep_max);
205
153
tree_elem->kernighan_rek(0,
206
rek_breite,rek_breite_anz,rek_deep_max,
207
funktion, param_list,param_anz,
208
pars_start, pars_start, pars_prev,
209
searchflag,&better_tree_found);
154
rek_breite, rek_breite_anz, rek_deep_max,
155
funktion, param_list, param_anz,
156
pars_start, pars_start, pars_org,
157
searchflag, &better_tree_found);
211
159
if (better_tree_found) {
212
160
ap_main->clear();
213
pars_start = GLOBAL_NT->tree->tree_root->costs();
215
sprintf(buffer,"New Parsimony: %f",pars_start);
216
abort_flag |= aw_status(buffer);
161
pars_start = get_root_node()->costs();
162
progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", pars_start, pars_org-pars_start));
224
ap_global_abort_flag |= abort_flag;
225
printf("Combines: %li\n", AP_sequence::global_combineCount);
171
printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount);
229
void PARS_optimizer_cb(AP_tree *tree) {
230
AP_tree *oldrootleft = GLOBAL_NT->tree->tree_root->leftson;
231
AP_tree *oldrootright = GLOBAL_NT->tree->tree_root->rightson;
233
for (ap_global_abort_flag = 0;!ap_global_abort_flag;){
234
AP_FLOAT old_pars = GLOBAL_NT->tree->tree_root->costs();
236
((AP_tree_nlen *)tree)->nn_interchange_rek(AP_TRUE,ap_global_abort_flag,-1,AP_BL_NNI_ONLY, GB_TRUE); // only not hidden
237
if (ap_global_abort_flag) break;
238
PARS_kernighan_cb(tree);
239
if (old_pars == GLOBAL_NT->tree->tree_root->costs()) {
240
ap_global_abort_flag = 1;
176
void ArbParsimony::optimize_tree(AP_tree *at, arb_progress& progress) {
177
AP_tree *oldrootleft = get_root_node()->get_leftson();
178
AP_tree *oldrootright = get_root_node()->get_rightson();
179
const AP_FLOAT org_pars = get_root_node()->costs();
180
AP_FLOAT prev_pars = org_pars;
182
progress.subtitle(GBS_global_string("Old parsimony: %.1f", org_pars));
184
while (!progress.aborted()) {
185
AP_FLOAT nni_pars = DOWNCAST(AP_tree_nlen*, at)->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
187
if (nni_pars == prev_pars) { // NNI did not reduce costs -> kern-lin
188
kernighan_optimize_tree(at);
189
AP_FLOAT ker_pars = get_root_node()->costs();
190
if (ker_pars == prev_pars) break; // kern-lin did not improve tree -> done
191
prev_pars = ker_pars;
194
prev_pars = nni_pars;
196
progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", prev_pars, org_pars-prev_pars));
243
199
if (oldrootleft->father == oldrootright) oldrootleft->set_root();
244
200
else oldrootright->set_root();
245
GLOBAL_NT->tree->tree_root->costs();
249
AWT_graphic_parsimony::AWT_graphic_parsimony(AW_root *root, GBDATA *gb_maini):AWT_graphic_tree(root,gb_maini)
252
AWT_graphic_tree *PARS_generate_tree(AW_root *root) {
253
AWT_graphic_parsimony *apdt = new AWT_graphic_parsimony(root,GLOBAL_gb_main);
254
AP_tree_nlen *aptnl = new AP_tree_nlen(0);
256
apdt->init((AP_tree *)aptnl);
257
ap_main->tree_root = &apdt->tree_root;
264
AWT_graphic_parsimony::init_devices(AW_window *aww, AW_device *device, AWT_canvas* ntw, AW_CL cd2)
202
get_root_node()->costs();
205
AWT_graphic_parsimony::AWT_graphic_parsimony(ArbParsimony& parsimony_, GBDATA *gb_main_, AD_map_viewer_cb map_viewer_cb_)
206
: AWT_graphic_tree(parsimony_.get_awroot(), gb_main_, map_viewer_cb_),
207
parsimony(parsimony_)
210
void ArbParsimony::generate_tree(WeightedFilter *pars_weighted_filter) {
211
AliView *aliview = pars_generate_aliview(pars_weighted_filter);
212
AP_sequence *seq_templ = 0;
214
GBDATA *gb_main = aliview->get_gb_main();
216
GB_transaction ta(gb_main);
217
bool is_aa = GBT_is_alignment_protein(gb_main, aliview->get_aliname());
219
if (is_aa) seq_templ = new AP_sequence_protein(aliview);
220
else seq_templ = new AP_sequence_parsimony(aliview);
223
tree = new AWT_graphic_parsimony(*this, aliview->get_gb_main(), PARS_map_viewer);
224
tree->init(new AP_TreeNlenNodeFactory, aliview, seq_templ, true, false);
225
ap_main->set_tree_root(tree);
228
AW_gc_manager AWT_graphic_parsimony::init_devices(AW_window *aww, AW_device *device, AWT_canvas* ntw) {
266
229
AW_init_color_group_defaults("arb_pars");
268
AW_gc_manager preset_window =
269
AW_manage_GC(aww,device,AWT_GC_CURSOR, AWT_GC_MAX, /*AWT_GC_CURSOR+7,*/ AW_GCM_DATA_AREA,
270
(AW_CB)AWT_resize_cb, (AW_CL)ntw, cd2,
231
AW_gc_manager gc_manager =
233
ntw->get_gc_base_name(),
234
device, AWT_GC_CURSOR, AWT_GC_MAX, /* AWT_GC_CURSOR+7, */ AW_GCM_DATA_AREA,
235
makeWindowCallback(AWT_resize_cb, ntw),
271
236
true, // uses color groups
291
256
"--unused", "--unused",
294
return preset_window;
297
void AWT_graphic_parsimony::show(AW_device *device)
300
if (GLOBAL_NT->tree->tree_root) parsval = (long)GLOBAL_NT->tree->tree_root->costs();
301
GLOBAL_NT->awr->awar(AWAR_PARSIMONY)->write_int( parsval);
302
long best = GLOBAL_NT->awr->awar(AWAR_BEST_PARSIMONY)->read_int();
303
if (parsval < best || 0==best) {
304
GLOBAL_NT->awr->awar(AWAR_BEST_PARSIMONY)->write_int( parsval);
306
this->AWT_graphic_tree::show(device);
310
void AWT_graphic_parsimony::command(AW_device *device, AWT_COMMAND_MODE cmd, int button, AW_key_mod key_modifier, AW_key_code key_code, char key_char,
311
AW_event_type type, AW_pos x, AW_pos y,
312
AW_clicked_line *cl, AW_clicked_text *ct)
314
static int bl_drag_flag;
319
bool compute_tree = false;
320
bool recalc_branch_lengths = false;
321
bool beautify_tree = false;
323
//GBDATA *gb_main = this->tree_static->gb_main;
325
case AWT_MODE_MOVE: // two point commands !!!!!
326
if(button==AWT_M_MIDDLE){
331
if( !(cl && cl->exists) ){
335
/*** check security level @@@ ***/
336
at = (AP_tree *)cl->client_data1;
337
if(at && at->father){
341
this->old_rot_cl = *cl;
346
if( bl_drag_flag && this->rot_at && this->rot_at->father){
347
this->rot_show_line(device);
353
this->rot_show_line(device);
356
case AW_Mouse_Release:
357
if( bl_drag_flag && this->rot_at && this->rot_at->father){
358
this->rot_show_line(device);
360
if (cl->exists) dest = (AP_tree *)cl->client_data1;
361
AP_tree *source = rot_at;
362
if (!(source && dest)) {
363
aw_message("Please Drag Line to a valid branch");
366
if (source == dest) {
367
aw_message("Please drag mouse from source to destination");
370
// if ( dest->is_son(source)) {
371
// aw_message("This operation is only allowed with two independent subtrees");
374
const char *error = 0;
377
// case AWT_MODE_MOVE:
381
error = source->cantMoveTo(dest);
383
source->moveTo(dest,cl->length);
384
recalc_branch_lengths = true;
388
error = source->move_group_info(dest);
391
error = "????? 45338";
395
// error = "????? 45338";
398
if (error) aw_message(error);
399
this->exports.refresh = 1;
400
this->exports.save = 1;
401
this->exports.resize = 1;
402
this->tree_root->test_tree();
403
//this->tree_root->compute_tree(gb_main);
262
void AWT_graphic_parsimony::show(AW_device *device) {
263
AP_tree_nlen *root_node = parsimony.get_root_node();
264
AW_awar *awar_pars = aw_root->awar(AWAR_PARSIMONY);
265
AW_awar *awar_best = aw_root->awar(AWAR_BEST_PARSIMONY);
267
long parsval = root_node ? root_node->costs() : 0;
268
awar_pars->write_int(parsval);
270
long best_pars = awar_best->read_int();
271
if (parsval < best_pars || 0==best_pars) awar_best->write_int(parsval);
273
AWT_graphic_tree::show(device);
276
void AWT_graphic_parsimony::handle_command(AW_device *device, AWT_graphic_event& event) {
277
ClickedTarget clicked(this, event.best_click());
278
bool recalc_branchlengths_on_structure_change = true;
280
switch (event.cmd()) {
281
// @@@ something is designed completely wrong here!
282
// why do all commands close TA and reopen when done?
413
284
case AWT_MODE_NNI:
414
if(type==AW_Mouse_Press){
285
if (event.type()==AW_Mouse_Press) {
415
286
GB_pop_transaction(gb_main);
418
if (!cl->exists) break;
419
at = (AP_tree *)cl->client_data1;
420
ap_global_abort_flag = AP_FALSE;
421
((AP_tree_nlen *)at)->nn_interchange_rek(AP_TRUE,ap_global_abort_flag,-1);
422
this->exports.refresh = 1;
423
this->exports.save = 1;
424
this->tree_root->test_tree();
425
recalc_branch_lengths = true;
428
AP_sequence::global_combineCount = 0;
429
ap_global_abort_flag = AP_FALSE;
430
((AP_tree_nlen *)this->tree_root)->nn_interchange_rek(AP_TRUE,ap_global_abort_flag,-1);
431
printf("Combines: %li\n", AP_sequence::global_combineCount);
432
this->exports.refresh = 1;
433
this->exports.save = 1;
434
this->tree_root->test_tree();
435
recalc_branch_lengths = true;
287
switch (event.button()) {
288
case AW_BUTTON_LEFT: {
289
if (clicked.node()) {
290
arb_progress progress("NNI optimize subtree");
291
AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, clicked.node());
292
atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
294
ASSERT_VALID_TREE(get_root_node());
298
case AW_BUTTON_RIGHT: {
299
arb_progress progress("NNI optimize tree");
300
long prevCombineCount = AP_sequence::combine_count();
302
AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, get_root_node());
303
atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false);
304
printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount);
307
ASSERT_VALID_TREE(get_root_node());
438
313
GB_begin_transaction(gb_main);
441
316
case AWT_MODE_KERNINGHAN:
442
if(type==AW_Mouse_Press){
317
if (event.type()==AW_Mouse_Press) {
443
318
GB_pop_transaction(gb_main);
446
if (!cl->exists) break;
447
at = (AP_tree *)cl->client_data1;
448
PARS_kernighan_cb(at);
449
this->exports.refresh = 1;
450
this->exports.save = 1;
451
this->tree_root->test_tree();
452
recalc_branch_lengths = true;
455
PARS_kernighan_cb(this->tree_root);
456
this->exports.refresh = 1;
457
this->exports.save = 1;
458
this->tree_root->test_tree();
459
recalc_branch_lengths = true;
319
switch (event.button()) {
321
if (clicked.node()) {
322
arb_progress progress("Kernighan-Lin optimize subtree");
323
parsimony.kernighan_optimize_tree(clicked.node());
324
this->exports.save = 1;
325
ASSERT_VALID_TREE(get_root_node());
328
case AW_BUTTON_RIGHT: {
329
arb_progress progress("Kernighan-Lin optimize tree");
330
parsimony.kernighan_optimize_tree(get_root_node());
331
this->exports.save = 1;
332
ASSERT_VALID_TREE(get_root_node());
462
337
GB_begin_transaction(gb_main);
465
340
case AWT_MODE_OPTIMIZE:
466
if(type==AW_Mouse_Press){
341
if (event.type()==AW_Mouse_Press) {
467
342
GB_pop_transaction(gb_main);
470
if (!cl->exists) break;
471
at = (AP_tree *)cl->client_data1;
472
if (at) PARS_optimizer_cb(at);
473
this->exports.refresh = 1;
474
this->exports.save = 1;
475
this->tree_root->test_tree();
476
recalc_branch_lengths = true;
479
PARS_optimizer_cb(this->tree_root);
480
this->exports.refresh = 1;
481
this->exports.save = 1;
482
this->tree_root->test_tree();
483
recalc_branch_lengths = true;
343
switch (event.button()) {
345
if (clicked.node()) {
346
arb_progress progress("Optimizing subtree");
347
parsimony.optimize_tree(clicked.node(), progress);
348
this->exports.save = 1;
349
ASSERT_VALID_TREE(get_root_node());
352
case AW_BUTTON_RIGHT: {
353
arb_progress progress("Optimizing tree");
355
parsimony.optimize_tree(get_root_node(), progress);
356
this->exports.save = 1;
357
ASSERT_VALID_TREE(get_root_node());
486
362
GB_begin_transaction(gb_main);
492
AWT_graphic_tree::command(device,cmd,button, key_modifier, key_code, key_char, type, x, y, cl, ct);
367
recalc_branchlengths_on_structure_change = false;
368
// fall-through (modes listed below trigger branchlength calculation)
370
AWT_graphic_tree::handle_command(device, event);
496
if (recalc_branch_lengths) {
497
int abort_flag = AP_FALSE;
498
rootEdge()->nni_rek(AP_FALSE,abort_flag,-1, GB_FALSE,AP_BL_BL_ONLY);
500
beautify_tree = true; // beautify after recalc_branch_lengths
504
this->resort_tree(0);
505
this->exports.save = 1;
507
compute_tree = true; // compute_tree after beautify_tree
511
this->tree_root->compute_tree(gb_main);
512
this->exports.refresh = 1;
374
if (exports.save == 1 && recalc_branchlengths_on_structure_change) {
375
arb_progress progress("Recalculating branch lengths");
376
rootEdge()->calc_branchlengths();
377
reorder_tree(BIG_BRANCHES_TO_TOP); // beautify after recalc_branch_lengths
382
// --------------------------------------------------------------------------------
385
#include <arb_diff.h>
387
#include <test_unit.h>
390
void fake_AW_init_color_groups();
392
struct fake_agt : public AWT_graphic_tree, virtual Noncopyable {
393
AP_sequence_parsimony *templ;
396
: AWT_graphic_tree(NULL, GLOBAL_gb_main, NULL),
403
void init(AliView *aliview) {
404
fake_AW_init_color_groups(); // acts like no species has a color
406
templ = aliview->has_data() ? new AP_sequence_parsimony(aliview) : NULL;
407
AWT_graphic_tree::init(new AP_TreeNlenNodeFactory, aliview, templ, true, false);
411
class PARSIMONY_testenv : virtual Noncopyable {
416
void common_init(const char *dbname) {
417
GLOBAL_gb_main = NULL;
420
TEST_EXPECT_NULL(ap_main);
424
apMain.set_tree_root(agt);
428
PARSIMONY_testenv(const char *dbname) {
430
agt->init(new AliView(GLOBAL_gb_main));
432
PARSIMONY_testenv(const char *dbname, const char *aliName) {
434
GB_transaction ta(GLOBAL_gb_main);
435
size_t aliLength = GBT_get_alignment_len(GLOBAL_gb_main, aliName);
437
AP_filter filter(aliLength);
438
if (!filter.is_invalid()) {
439
AP_weights weights(&filter);
440
agt->init(new AliView(GLOBAL_gb_main, filter, weights, aliName));
443
~PARSIMONY_testenv() {
444
TEST_EXPECT_EQUAL(ap_main, &apMain);
448
GB_close(GLOBAL_gb_main);
449
GLOBAL_gb_main = NULL;
452
GB_ERROR load_tree(const char *tree_name) {
453
GB_transaction ta(GLOBAL_gb_main); // @@@ do inside AWT_graphic_tree::load?
454
return agt->load(GLOBAL_gb_main, tree_name, 0, 0);
456
AP_tree_nlen *tree_root() { return apMain.get_root_node(); }
458
void push() { apMain.push(); }
459
void pop() { apMain.pop(); }
462
void TEST_tree_modifications() {
463
PARSIMONY_testenv env("TEST_trees.arb");
464
TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
467
AP_tree_nlen *root = env.tree_root();
468
TEST_REJECT_NULL(root);
470
AP_tree_edge::initialize(root); // builds edges
471
TEST_ASSERT_VALID_TREE(root);
473
root->compute_tree();
475
// first check initial state:
477
AP_tree_members& root_info = root->gr;
479
TEST_EXPECT_EQUAL(root_info.grouped, 0);
480
TEST_EXPECT_EQUAL(root_info.hidden, 0);
481
TEST_EXPECT_EQUAL(root_info.has_marked_children, 1);
482
TEST_EXPECT_EQUAL(root_info.leaf_sum, 15);
484
TEST_EXPECT_SIMILAR(root_info.max_tree_depth, 1.624975, 0.000001);
485
TEST_EXPECT_SIMILAR(root_info.min_tree_depth, 0.341681, 0.000001);
487
GB_transaction ta(GLOBAL_gb_main);
488
GBT_mark_all(GLOBAL_gb_main, 0); // unmark all species
489
root->compute_tree();
490
TEST_EXPECT_EQUAL(root_info.has_marked_children, 0);
494
#define B1_TOP "(((((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,CloTyrob:0.009):0.274,CloInnoc:0.371):0.057,CloBifer:0.388):0.124"
495
#define B1_BOT "(CloBifer:0.388,(CloInnoc:0.371,(CloTyrob:0.009,(CloTyro2:0.017,(CloTyro3:1.046,CloTyro4:0.061):0.026):0.017):0.274):0.057):0.124"
496
#define B2_TOP "(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010,CloPaste:0.179):0.131"
497
#define B2_BOT "(CloPaste:0.179,(CloCarni:0.120,(CloButy2:0.009,CloButyr:0.000):0.564):0.010):0.131"
500
#define B3_LEFT_TOP_SONS "(((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053,CelBiazo:0.059)"
501
#define B3_TOP_SONS B3_LEFT_TOP_SONS ":0.207,CytAquat:0.711"
502
#define B3_TOP_SONS_CCR "((CorAquat:0.187,CorGluta:0.522):0.053,CelBiazo:0.059):0.207,CytAquat:0.711" // CCR = CurCitre removed
503
#define B3_TOP "(" B3_TOP_SONS "):0.081"
504
#define B3_BOT "(CytAquat:0.711,(CelBiazo:0.059,(CorGluta:0.522,(CorAquat:0.084,CurCitre:0.058):0.103):0.053):0.207):0.081"
507
const char *top_topo = "((" B1_TOP "," B2_TOP "):0.081," B3_TOP ");";
508
const char *edge_topo = "((" B1_TOP "," B2_BOT "):0.081," B3_BOT ");";
509
const char *bottom_topo = "(" B3_BOT ",(" B2_BOT "," B1_BOT "):0.081);";
511
const char *radial_topo =
512
"(((CloPaste:0.179,((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010):0.131,"
513
"((CloInnoc:0.371,((CloTyro2:0.017,(CloTyro3:1.046,CloTyro4:0.061):0.026):0.017,CloTyrob:0.009):0.274):0.057,CloBifer:0.388):0.124):0.081,"
514
"((CelBiazo:0.059,((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053):0.207,CytAquat:0.711):0.081);";
515
const char *radial_topo2 =
516
"(((CloBifer:0.388,(CloInnoc:0.371,(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,CloTyrob:0.009):0.274):0.057):0.124," B2_TOP "):0.081,"
517
"(CytAquat:0.711," B3_LEFT_TOP_SONS ":0.207):0.081);";
519
// expect that no mode reproduces another mode:
520
TEST_EXPECT_DIFFERENT(top_topo, edge_topo);
521
TEST_EXPECT_DIFFERENT(top_topo, bottom_topo);
522
TEST_EXPECT_DIFFERENT(top_topo, radial_topo);
523
TEST_EXPECT_DIFFERENT(top_topo, radial_topo2);
524
TEST_EXPECT_DIFFERENT(edge_topo, bottom_topo);
525
TEST_EXPECT_DIFFERENT(edge_topo, radial_topo);
526
TEST_EXPECT_DIFFERENT(edge_topo, radial_topo2);
527
TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo);
528
TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo2);
529
TEST_EXPECT_DIFFERENT(radial_topo, radial_topo2);
531
env.push(); // 1st stack level (=top_topo)
533
TEST_ASSERT_VALID_TREE(root);
535
TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
536
// test reorder_tree:
537
root->reorder_tree(BIG_BRANCHES_TO_EDGE); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo); env.push(); // 2nd stack level (=edge_topo)
538
root->reorder_tree(BIG_BRANCHES_TO_BOTTOM); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo); env.push(); // 3rd stack level (=bottom_topo)
539
root->reorder_tree(BIG_BRANCHES_TO_CENTER); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo);
540
root->reorder_tree(BIG_BRANCHES_ALTERNATING); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo2);
541
root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
543
TEST_ASSERT_VALID_TREE(root);
546
AP_tree_nlen *CloTyrob = root->findLeafNamed("CloTyrob");
547
TEST_REJECT_NULL(CloTyrob);
549
ARB_edge rootEdge(root->get_leftson(), root->get_rightson());
550
CloTyrob->set_root();
552
TEST_ASSERT_VALID_TREE(root);
554
const char *rootAtCloTyrob_topo =
556
"(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,"
557
"((((" B3_TOP_SONS "):0.162," B2_TOP "):0.124,CloBifer:0.388):0.057,CloInnoc:0.371):0.274):0.004);";
559
TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
560
env.push(); // 4th stack level (=rootAtCloTyrob_topo)
562
TEST_ASSERT_VALID_TREE(root);
564
AP_tree_nlen *CelBiazoFather = root->findLeafNamed("CelBiazo")->get_father();
565
TEST_REJECT_NULL(CelBiazoFather);
566
CelBiazoFather->set_root();
568
const char *rootAtCelBiazoFather_topo = "(" B3_LEFT_TOP_SONS ":0.104,((" B1_TOP "," B2_TOP "):0.162,CytAquat:0.711):0.104);";
569
TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCelBiazoFather_topo);
571
TEST_ASSERT_VALID_TREE(root);
573
ARB_edge oldRootEdge(rootEdge.source(), rootEdge.dest());
574
DOWNCAST(AP_tree_nlen*,oldRootEdge.son())->set_root();
576
const char *rootSetBack_topo = top_topo;
577
TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
578
env.push(); // 5th stack level (=rootSetBack_topo)
580
TEST_ASSERT_VALID_TREE(root);
583
AP_tree_nlen *CurCitre = root->findLeafNamed("CurCitre");
584
TEST_REJECT_NULL(CurCitre);
585
TEST_REJECT_NULL(CurCitre->get_father());
588
const char *CurCitre_removed_topo = "((" B1_TOP "," B2_TOP "):0.081,(" B3_TOP_SONS_CCR "):0.081);";
589
// ------------------------------------------------------------------- ^^^ = B3_TOP_SONS minus CurCitre
590
TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
592
TEST_ASSERT_VALID_TREE(root);
593
TEST_ASSERT_VALID_TREE(CurCitre);
595
TEST_EXPECT_EQUAL(root->gr.leaf_sum, 15); // out of date
596
root->compute_tree();
597
TEST_EXPECT_EQUAL(root->gr.leaf_sum, 14);
599
env.push(); // 6th stack level (=CurCitre_removed_topo)
601
TEST_ASSERT_VALID_TREE(root);
604
AP_tree_nlen *CloCarni = root->findLeafNamed("CloCarni");
605
TEST_REJECT_NULL(CloCarni);
606
CurCitre->insert(CloCarni); // this creates two extra edges (not destroyed by destroy() below) and one extra node
608
const char *CurCitre_inserted_topo = "((" B1_TOP ",(((CloButy2:0.009,CloButyr:0.000):0.564,(CurCitre:0.060,CloCarni:0.060):0.060):0.010,CloPaste:0.179):0.131):0.081,(" B3_TOP_SONS_CCR "):0.081);";
609
TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_inserted_topo);
611
AP_tree_nlen *node_del_manually = CurCitre->get_father();
612
AP_tree_edge *edge1_del_manually = CurCitre->edgeTo(node_del_manually);
613
AP_tree_edge *edge2_del_manually = CurCitre->get_brother()->edgeTo(node_del_manually);
615
TEST_ASSERT_VALID_TREE(root);
618
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
619
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
620
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
621
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo);
622
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo);
623
env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
625
TEST_ASSERT_VALID_TREE(root);
627
AP_tree_edge::destroy(root);
629
// delete memory allocated by insert() above and lost due to pop()s
630
delete edge1_del_manually;
631
delete edge2_del_manually;
633
node_del_manually->forget_origin();
634
node_del_manually->father = NULL;
635
node_del_manually->leftson = NULL;
636
node_del_manually->rightson = NULL;
637
delete node_del_manually;
646
void TEST_calc_bootstraps() {
647
PARSIMONY_testenv env("TEST_trees.arb", "ali_5s");
648
TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
650
const char *bs_origi_topo = "(((((((CloTyro3,CloTyro4)'40%',CloTyro2)'0%',CloTyrob)'97%',CloInnoc)'0%',CloBifer)'53%',(((CloButy2,CloButyr)'100%',CloCarni)'33%',CloPaste)'97%')'100%',((((CorAquat,CurCitre)'100%',CorGluta)'17%',CelBiazo)'40%',CytAquat)'100%');";
651
const char *bs_limit_topo = "(((((((CloTyro3,CloTyro4)'87%',CloTyro2)'0%',CloTyrob)'100%',CloInnoc)'87%',CloBifer)'83%',(((CloButy2,CloButyr)'99%',CloCarni)'17%',CloPaste)'56%')'61%',((((CorAquat,CurCitre)'78%',CorGluta)'0%',CelBiazo)'59%',CytAquat)'61%');";
652
const char *bs_estim_topo = "(((((((CloTyro3,CloTyro4)'75%',CloTyro2)'0%',CloTyrob)'100%',CloInnoc)'75%',CloBifer)'78%',(((CloButy2,CloButyr)'99%',CloCarni)'13%',CloPaste)'32%')'53%',((((CorAquat,CurCitre)'74%',CorGluta)'0%',CelBiazo)'56%',CytAquat)'53%');";
655
AP_tree_nlen *root = env.tree_root();
656
TEST_REJECT_NULL(root);
658
AP_tree_edge::initialize(root); // builds edges
659
TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that)
661
AP_tree_edge *root_edge = rootEdge();
662
TEST_REJECT_NULL(root_edge);
664
root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_origi_topo);
666
root_edge->nni_rek(-1, false, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_LIMIT), NULL); root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_limit_topo);
667
root_edge->nni_rek(-1, false, AP_BL_MODE(AP_BL_BL_ONLY|AP_BL_BOOTSTRAP_ESTIMATE), NULL); root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_estim_topo);
669
TEST_EXPECT_EQUAL(env.tree_root(), root);
670
AP_tree_edge::destroy(root);
674
void TEST_tree_remove_add_all() {
675
// reproduces crash as described in #527
676
PARSIMONY_testenv env("TEST_trees.arb", "ali_5s");
677
TEST_EXPECT_NO_ERROR(env.load_tree("tree_nj"));
680
AP_tree_nlen *leaf[LEAFS];
681
const char *name[LEAFS] = {
690
AP_tree_nlen *root = env.tree_root();
691
TEST_REJECT_NULL(root);
693
for (int i = 0; i<LEAFS; ++i) {
694
leaf[i] = root->findLeafNamed(name[i]);
695
TEST_REJECT_NULL(leaf[i]);
698
AP_tree_edge::initialize(root); // builds edges
699
TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that)
701
AP_tree_root *troot = leaf[0]->get_tree_root();
702
TEST_REJECT_NULL(troot);
704
// Note: following loop leaks father nodes and edges
705
for (int i = 0; i<LEAFS-1; ++i) { // removing the second to last leaf, "removes" both remaining leafs
706
TEST_ASSERT_VALID_TREE(root);
708
TEST_ASSERT_VALID_TREE(leaf[i]);
710
leaf[LEAFS-1]->father = NULL; // correct final leaf (not removed regularily)
712
leaf[0]->initial_insert(leaf[1], troot);
713
for (int i = 2; i<LEAFS; ++i) {
714
TEST_ASSERT_VALID_TREE(leaf[i-1]);
715
TEST_ASSERT_VALID_TREE(leaf[i]);
716
leaf[i]->insert(leaf[i-1]);
722
// --------------------------------------------------------------------------------
added
removed
PARSIMONY/PARS_dtree.cxx
