1 | // =============================================================== // |
---|
2 | // // |
---|
3 | // File : PARS_dtree.cxx // |
---|
4 | // Purpose : // |
---|
5 | // // |
---|
6 | // Institute of Microbiology (Technical University Munich) // |
---|
7 | // http://www.arb-home.de/ // |
---|
8 | // // |
---|
9 | // =============================================================== // |
---|
10 | |
---|
11 | #include "pars_dtree.hxx" |
---|
12 | #include "pars_main.hxx" |
---|
13 | #include "pars_debug.hxx" |
---|
14 | #include "ap_tree_nlen.hxx" |
---|
15 | |
---|
16 | #include <AP_seq_dna.hxx> |
---|
17 | #include <AP_seq_protein.hxx> |
---|
18 | #include <AP_filter.hxx> |
---|
19 | |
---|
20 | #include <ColumnStat.hxx> |
---|
21 | #include <awt_sel_boxes.hxx> |
---|
22 | #include <awt_filter.hxx> |
---|
23 | |
---|
24 | #include <gui_aliview.hxx> |
---|
25 | |
---|
26 | #include <aw_preset.hxx> |
---|
27 | #include <aw_awar.hxx> |
---|
28 | #include <aw_msg.hxx> |
---|
29 | #include <arb_progress.h> |
---|
30 | #include <aw_root.hxx> |
---|
31 | #include <aw_question.hxx> |
---|
32 | |
---|
33 | #include "PerfMeter.h" |
---|
34 | |
---|
35 | static void AWT_graphic_parsimony_root_changed(void *cd, AP_tree *old, AP_tree *newroot) { |
---|
36 | AWT_graphic_tree *agt = (AWT_graphic_tree*)cd; |
---|
37 | |
---|
38 | if (old == agt->displayed_root) agt->displayed_root = newroot; |
---|
39 | } |
---|
40 | |
---|
41 | static AliView *pars_generate_aliview(WeightedFilter *pars_weighted_filter) { |
---|
42 | GBDATA *gb_main = pars_weighted_filter->get_gb_main(); |
---|
43 | char *ali_name; |
---|
44 | { |
---|
45 | GB_transaction ta(gb_main); |
---|
46 | ali_name = GBT_read_string(gb_main, AWAR_ALIGNMENT); |
---|
47 | } |
---|
48 | GB_ERROR error = NULL; |
---|
49 | AliView *aliview = pars_weighted_filter->create_aliview(ali_name, error); |
---|
50 | if (!aliview) aw_popup_exit(error); |
---|
51 | free(ali_name); |
---|
52 | return aliview; |
---|
53 | } |
---|
54 | |
---|
55 | void PARS_tree_init(AWT_graphic_tree *agt) { |
---|
56 | ap_assert(agt->get_root_node()); |
---|
57 | ap_assert(agt == ap_main->get_tree_root()); |
---|
58 | |
---|
59 | GB_transaction ta(GLOBAL_gb_main); |
---|
60 | |
---|
61 | const char *use = ap_main->get_aliname(); |
---|
62 | long ali_len = GBT_get_alignment_len(GLOBAL_gb_main, use); |
---|
63 | if (ali_len <= 1) { |
---|
64 | aw_popup_exit("No valid alignment selected! Try again"); |
---|
65 | } |
---|
66 | |
---|
67 | agt->tree_static->set_root_changed_callback(AWT_graphic_parsimony_root_changed, agt); |
---|
68 | } |
---|
69 | |
---|
70 | static double funktion_quadratisch(double wert, double *param_list, int param_anz) { |
---|
71 | if (param_anz != 3) { |
---|
72 | ap_assert(0); // wrong number of parameters |
---|
73 | return 0; |
---|
74 | } |
---|
75 | return wert * wert * param_list[0] + wert * param_list[1] + param_list[2]; |
---|
76 | } |
---|
77 | |
---|
78 | |
---|
79 | void ArbParsimony::kernighan_optimize_tree(AP_tree *at) { |
---|
80 | GB_push_transaction(GLOBAL_gb_main); |
---|
81 | |
---|
82 | long prevCombineCount = AP_sequence::combine_count(); |
---|
83 | |
---|
84 | AP_FLOAT pars_start = get_root_node()->costs(); |
---|
85 | const AP_FLOAT pars_org = pars_start; |
---|
86 | |
---|
87 | int rek_deep_max = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/maxdepth"); |
---|
88 | |
---|
89 | AP_KL_FLAG funktype = (AP_KL_FLAG)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/function_type"); |
---|
90 | |
---|
91 | int param_anz; |
---|
92 | double param_list[3]; |
---|
93 | double f_startx, f_maxy, f_maxx, f_max_deep; |
---|
94 | f_max_deep = (double)rek_deep_max; // x2 |
---|
95 | f_startx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/start"); |
---|
96 | f_maxy = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxy"); |
---|
97 | f_maxx = (double)*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/maxx"); |
---|
98 | |
---|
99 | double (*funktion)(double wert, double *param_list, int param_anz); |
---|
100 | switch (funktype) { |
---|
101 | default: |
---|
102 | case AP_QUADRAT_START: |
---|
103 | funktion = funktion_quadratisch; |
---|
104 | param_anz = 3; |
---|
105 | param_list[2] = f_startx; |
---|
106 | param_list[0] = (f_startx - f_maxy) / (f_maxx * f_maxx); |
---|
107 | param_list[1] = -2.0 * param_list[0] * f_maxx; |
---|
108 | break; |
---|
109 | case AP_QUADRAT_MAX: // parameter liste fuer quadratische gleichung (y =ax^2 +bx +c) |
---|
110 | funktion = funktion_quadratisch; |
---|
111 | param_anz = 3; |
---|
112 | param_list[0] = - f_maxy / ((f_max_deep - f_maxx) * (f_max_deep - f_maxx)); |
---|
113 | param_list[1] = -2.0 * param_list[0] * f_maxx; |
---|
114 | param_list[2] = f_maxy + param_list[0] * f_maxx * f_maxx; |
---|
115 | break; |
---|
116 | } |
---|
117 | |
---|
118 | |
---|
119 | AP_KL_FLAG searchflag=(AP_KL_FLAG)0; |
---|
120 | if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/dynamic/enable")) { |
---|
121 | searchflag = AP_DYNAMIK; |
---|
122 | } |
---|
123 | if (*GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/enable")) { |
---|
124 | searchflag = (AP_KL_FLAG)(searchflag|AP_STATIC); |
---|
125 | } |
---|
126 | |
---|
127 | int rek_breite[8]; |
---|
128 | rek_breite[0] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth0"); |
---|
129 | rek_breite[1] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth1"); |
---|
130 | rek_breite[2] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth2"); |
---|
131 | rek_breite[3] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth3"); |
---|
132 | rek_breite[4] = *GBT_read_int(GLOBAL_gb_main, "genetic/kh/static/depth4"); |
---|
133 | int rek_breite_anz = 5; |
---|
134 | |
---|
135 | int anzahl = (int)(*GBT_read_float(GLOBAL_gb_main, "genetic/kh/nodes")*at->count_leafs()); |
---|
136 | AP_tree **list = at->getRandomNodes(anzahl); |
---|
137 | |
---|
138 | arb_progress progress(anzahl); |
---|
139 | |
---|
140 | progress.subtitle(GBS_global_string("Old Parsimony: %.1f", pars_start)); |
---|
141 | |
---|
142 | GB_pop_transaction(GLOBAL_gb_main); |
---|
143 | |
---|
144 | for (int i=0; i<anzahl && !progress.aborted(); i++) { |
---|
145 | AP_tree_nlen *tree_elem = DOWNCAST(AP_tree_nlen*, list[i]); // @@@ pass 'at' as AP_tree_nlen |
---|
146 | |
---|
147 | bool in_folded_group = tree_elem->gr.hidden || |
---|
148 | (tree_elem->father && tree_elem->get_father()->gr.hidden); |
---|
149 | |
---|
150 | if (!in_folded_group) { |
---|
151 | bool better_tree_found = false; |
---|
152 | ap_main->push(); |
---|
153 | display_clear(funktion, param_list, param_anz, (int)pars_start, (int)rek_deep_max); |
---|
154 | |
---|
155 | tree_elem->kernighan_rek(0, |
---|
156 | rek_breite, rek_breite_anz, rek_deep_max, |
---|
157 | funktion, param_list, param_anz, |
---|
158 | pars_start, pars_start, pars_org, |
---|
159 | searchflag, &better_tree_found); |
---|
160 | |
---|
161 | if (better_tree_found) { |
---|
162 | ap_main->clear(); |
---|
163 | pars_start = get_root_node()->costs(); |
---|
164 | progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", pars_start, pars_org-pars_start)); |
---|
165 | } |
---|
166 | else { |
---|
167 | ap_main->pop(); |
---|
168 | } |
---|
169 | } |
---|
170 | progress.inc(); |
---|
171 | } |
---|
172 | delete list; |
---|
173 | printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount); |
---|
174 | } |
---|
175 | |
---|
176 | |
---|
177 | |
---|
178 | void ArbParsimony::optimize_tree(AP_tree *at, arb_progress& progress) { |
---|
179 | AP_tree *oldrootleft = get_root_node()->get_leftson(); |
---|
180 | AP_tree *oldrootright = get_root_node()->get_rightson(); |
---|
181 | const AP_FLOAT org_pars = get_root_node()->costs(); |
---|
182 | AP_FLOAT prev_pars = org_pars; |
---|
183 | |
---|
184 | OptiPerfMeter perf("optimize_tree", org_pars); |
---|
185 | |
---|
186 | progress.subtitle(GBS_global_string("Old parsimony: %.1f", org_pars)); |
---|
187 | |
---|
188 | while (!progress.aborted()) { |
---|
189 | AP_FLOAT nni_pars = DOWNCAST(AP_tree_nlen*, at)->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false); |
---|
190 | |
---|
191 | if (nni_pars == prev_pars) { // NNI did not reduce costs -> kern-lin |
---|
192 | kernighan_optimize_tree(at); |
---|
193 | AP_FLOAT ker_pars = get_root_node()->costs(); |
---|
194 | if (ker_pars == prev_pars) break; // kern-lin did not improve tree -> done |
---|
195 | prev_pars = ker_pars; |
---|
196 | } |
---|
197 | else { |
---|
198 | prev_pars = nni_pars; |
---|
199 | } |
---|
200 | progress.subtitle(GBS_global_string("New parsimony: %.1f (gain: %.1f)", prev_pars, org_pars-prev_pars)); |
---|
201 | } |
---|
202 | |
---|
203 | if (oldrootleft->father == oldrootright) oldrootleft->set_root(); |
---|
204 | else oldrootright->set_root(); |
---|
205 | |
---|
206 | get_root_node()->costs(); |
---|
207 | |
---|
208 | perf.dump(stdout, prev_pars); |
---|
209 | } |
---|
210 | |
---|
211 | AWT_graphic_parsimony::AWT_graphic_parsimony(ArbParsimony& parsimony_, GBDATA *gb_main_, AD_map_viewer_cb map_viewer_cb_) |
---|
212 | : AWT_graphic_tree(parsimony_.get_awroot(), gb_main_, map_viewer_cb_), |
---|
213 | parsimony(parsimony_) |
---|
214 | {} |
---|
215 | |
---|
216 | void ArbParsimony::generate_tree(WeightedFilter *pars_weighted_filter) { |
---|
217 | AliView *aliview = pars_generate_aliview(pars_weighted_filter); |
---|
218 | AP_sequence *seq_templ = 0; |
---|
219 | |
---|
220 | GBDATA *gb_main = aliview->get_gb_main(); |
---|
221 | { |
---|
222 | GB_transaction ta(gb_main); |
---|
223 | bool is_aa = GBT_is_alignment_protein(gb_main, aliview->get_aliname()); |
---|
224 | |
---|
225 | if (is_aa) seq_templ = new AP_sequence_protein(aliview); |
---|
226 | else seq_templ = new AP_sequence_parsimony(aliview); |
---|
227 | } |
---|
228 | |
---|
229 | tree = new AWT_graphic_parsimony(*this, aliview->get_gb_main(), PARS_map_viewer); |
---|
230 | tree->init(new AP_TreeNlenNodeFactory, aliview, seq_templ, true, false); |
---|
231 | ap_main->set_tree_root(tree); |
---|
232 | } |
---|
233 | |
---|
234 | AW_gc_manager AWT_graphic_parsimony::init_devices(AW_window *aww, AW_device *device, AWT_canvas* ntw) { |
---|
235 | AW_init_color_group_defaults("arb_pars"); |
---|
236 | |
---|
237 | AW_gc_manager gc_manager = |
---|
238 | AW_manage_GC(aww, |
---|
239 | ntw->get_gc_base_name(), |
---|
240 | device, AWT_GC_CURSOR, AWT_GC_MAX, /* AWT_GC_CURSOR+7, */ AW_GCM_DATA_AREA, |
---|
241 | makeWindowCallback(AWT_resize_cb, ntw), |
---|
242 | true, // uses color groups |
---|
243 | "#AAAA55", |
---|
244 | |
---|
245 | // Important note : |
---|
246 | // Many gc indices are shared between ABR_NTREE and ARB_PARSIMONY |
---|
247 | // e.g. the tree drawing routines use same gc's for drawing both trees |
---|
248 | // (check AWT_dtree.cxx AWT_graphic_tree::init_devices) |
---|
249 | |
---|
250 | "Cursor$#FFFFFF", |
---|
251 | "Branch remarks$#DBE994", |
---|
252 | "+-Bootstrap$#DBE994", "-B.(limited)$white", |
---|
253 | "--unused$#ff0000", |
---|
254 | "Marked$#FFFF00", |
---|
255 | "Some marked$#eeee88", |
---|
256 | "Not marked$black", |
---|
257 | "Zombies etc.$#cc5924", |
---|
258 | |
---|
259 | "--unused", "--unused", // these reserve the numbers which are used for probe colors in ARB_NTREE |
---|
260 | "--unused", "--unused", // (this is necessary because ARB_PARS and ARB_NTREE use the same tree painting routines) |
---|
261 | "--unused", "--unused", |
---|
262 | "--unused", "--unused", |
---|
263 | |
---|
264 | NULL); |
---|
265 | return gc_manager; |
---|
266 | } |
---|
267 | |
---|
268 | void AWT_graphic_parsimony::show(AW_device *device) { |
---|
269 | AP_tree_nlen *root_node = parsimony.get_root_node(); |
---|
270 | AW_awar *awar_pars = aw_root->awar(AWAR_PARSIMONY); |
---|
271 | AW_awar *awar_best = aw_root->awar(AWAR_BEST_PARSIMONY); |
---|
272 | |
---|
273 | long parsval = root_node ? root_node->costs() : 0; |
---|
274 | awar_pars->write_int(parsval); |
---|
275 | |
---|
276 | long best_pars = awar_best->read_int(); |
---|
277 | if (parsval < best_pars || 0==best_pars) awar_best->write_int(parsval); |
---|
278 | |
---|
279 | AWT_graphic_tree::show(device); |
---|
280 | } |
---|
281 | |
---|
282 | void AWT_graphic_parsimony::handle_command(AW_device *device, AWT_graphic_event& event) { |
---|
283 | ClickedTarget clicked(this, event.best_click()); |
---|
284 | bool recalc_branchlengths_on_structure_change = true; |
---|
285 | |
---|
286 | switch (event.cmd()) { |
---|
287 | // @@@ something is designed completely wrong here! |
---|
288 | // why do all commands close TA and reopen when done? |
---|
289 | |
---|
290 | case AWT_MODE_NNI: |
---|
291 | if (event.type()==AW_Mouse_Press) { |
---|
292 | GB_pop_transaction(gb_main); |
---|
293 | switch (event.button()) { |
---|
294 | case AW_BUTTON_LEFT: { |
---|
295 | if (clicked.node()) { |
---|
296 | arb_progress progress("NNI optimize subtree"); |
---|
297 | AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, clicked.node()); |
---|
298 | atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false); |
---|
299 | exports.save = 1; |
---|
300 | ASSERT_VALID_TREE(get_root_node()); |
---|
301 | } |
---|
302 | break; |
---|
303 | } |
---|
304 | case AW_BUTTON_RIGHT: { |
---|
305 | arb_progress progress("NNI optimize tree"); |
---|
306 | long prevCombineCount = AP_sequence::combine_count(); |
---|
307 | |
---|
308 | AP_tree_nlen *atn = DOWNCAST(AP_tree_nlen*, get_root_node()); |
---|
309 | atn->nn_interchange_rek(-1, AP_BL_NNI_ONLY, false); |
---|
310 | printf("Combines: %li\n", AP_sequence::combine_count()-prevCombineCount); |
---|
311 | |
---|
312 | exports.save = 1; |
---|
313 | ASSERT_VALID_TREE(get_root_node()); |
---|
314 | break; |
---|
315 | } |
---|
316 | |
---|
317 | default: break; |
---|
318 | } |
---|
319 | GB_begin_transaction(gb_main); |
---|
320 | } |
---|
321 | break; |
---|
322 | case AWT_MODE_KERNINGHAN: |
---|
323 | if (event.type()==AW_Mouse_Press) { |
---|
324 | GB_pop_transaction(gb_main); |
---|
325 | switch (event.button()) { |
---|
326 | case AW_BUTTON_LEFT: |
---|
327 | if (clicked.node()) { |
---|
328 | arb_progress progress("Kernighan-Lin optimize subtree"); |
---|
329 | parsimony.kernighan_optimize_tree(clicked.node()); |
---|
330 | this->exports.save = 1; |
---|
331 | ASSERT_VALID_TREE(get_root_node()); |
---|
332 | } |
---|
333 | break; |
---|
334 | case AW_BUTTON_RIGHT: { |
---|
335 | arb_progress progress("Kernighan-Lin optimize tree"); |
---|
336 | parsimony.kernighan_optimize_tree(get_root_node()); |
---|
337 | this->exports.save = 1; |
---|
338 | ASSERT_VALID_TREE(get_root_node()); |
---|
339 | break; |
---|
340 | } |
---|
341 | default: break; |
---|
342 | } |
---|
343 | GB_begin_transaction(gb_main); |
---|
344 | } |
---|
345 | break; |
---|
346 | case AWT_MODE_OPTIMIZE: |
---|
347 | if (event.type()==AW_Mouse_Press) { |
---|
348 | GB_pop_transaction(gb_main); |
---|
349 | switch (event.button()) { |
---|
350 | case AW_BUTTON_LEFT: |
---|
351 | if (clicked.node()) { |
---|
352 | arb_progress progress("Optimizing subtree"); |
---|
353 | parsimony.optimize_tree(clicked.node(), progress); |
---|
354 | this->exports.save = 1; |
---|
355 | ASSERT_VALID_TREE(get_root_node()); |
---|
356 | } |
---|
357 | break; |
---|
358 | case AW_BUTTON_RIGHT: { |
---|
359 | arb_progress progress("Optimizing tree"); |
---|
360 | |
---|
361 | parsimony.optimize_tree(get_root_node(), progress); |
---|
362 | this->exports.save = 1; |
---|
363 | ASSERT_VALID_TREE(get_root_node()); |
---|
364 | break; |
---|
365 | } |
---|
366 | default: break; |
---|
367 | } |
---|
368 | GB_begin_transaction(gb_main); |
---|
369 | } |
---|
370 | break; |
---|
371 | |
---|
372 | default: |
---|
373 | recalc_branchlengths_on_structure_change = false; |
---|
374 | // fall-through (modes listed below trigger branchlength calculation) |
---|
375 | case AWT_MODE_MOVE: |
---|
376 | AWT_graphic_tree::handle_command(device, event); |
---|
377 | break; |
---|
378 | } |
---|
379 | |
---|
380 | if (exports.save == 1 && recalc_branchlengths_on_structure_change) { |
---|
381 | arb_progress progress("Recalculating branch lengths"); |
---|
382 | rootEdge()->calc_branchlengths(); |
---|
383 | reorder_tree(BIG_BRANCHES_TO_TOP); // beautify after recalc_branch_lengths |
---|
384 | } |
---|
385 | } |
---|
386 | |
---|
387 | |
---|
388 | // -------------------------------------------------------------------------------- |
---|
389 | |
---|
390 | #ifdef UNIT_TESTS |
---|
391 | #include <arb_diff.h> |
---|
392 | #ifndef TEST_UNIT_H |
---|
393 | #include <test_unit.h> |
---|
394 | #endif |
---|
395 | |
---|
396 | void fake_AW_init_color_groups(); |
---|
397 | |
---|
398 | struct fake_agt : public AWT_graphic_tree, virtual Noncopyable { |
---|
399 | AP_sequence_parsimony *templ; |
---|
400 | |
---|
401 | fake_agt() |
---|
402 | : AWT_graphic_tree(NULL, GLOBAL_gb_main, NULL), |
---|
403 | templ(NULL) |
---|
404 | { |
---|
405 | } |
---|
406 | ~fake_agt() { |
---|
407 | delete templ; |
---|
408 | } |
---|
409 | void init(AliView *aliview) { |
---|
410 | fake_AW_init_color_groups(); // acts like no species has a color |
---|
411 | delete templ; |
---|
412 | templ = aliview->has_data() ? new AP_sequence_parsimony(aliview) : NULL; |
---|
413 | AWT_graphic_tree::init(new AP_TreeNlenNodeFactory, aliview, templ, true, false); |
---|
414 | } |
---|
415 | }; |
---|
416 | |
---|
417 | class PARSIMONY_testenv : virtual Noncopyable { |
---|
418 | GB_shell shell; |
---|
419 | AP_main apMain; |
---|
420 | fake_agt *agt; |
---|
421 | |
---|
422 | void common_init(const char *dbname) { |
---|
423 | GLOBAL_gb_main = NULL; |
---|
424 | apMain.open(dbname); |
---|
425 | |
---|
426 | TEST_EXPECT_NULL(ap_main); |
---|
427 | ap_main = &apMain; |
---|
428 | |
---|
429 | agt = new fake_agt; |
---|
430 | apMain.set_tree_root(agt); |
---|
431 | } |
---|
432 | |
---|
433 | public: |
---|
434 | PARSIMONY_testenv(const char *dbname) { |
---|
435 | common_init(dbname); |
---|
436 | agt->init(new AliView(GLOBAL_gb_main)); |
---|
437 | } |
---|
438 | PARSIMONY_testenv(const char *dbname, const char *aliName) { |
---|
439 | common_init(dbname); |
---|
440 | GB_transaction ta(GLOBAL_gb_main); |
---|
441 | size_t aliLength = GBT_get_alignment_len(GLOBAL_gb_main, aliName); |
---|
442 | |
---|
443 | AP_filter filter(aliLength); |
---|
444 | if (!filter.is_invalid()) { |
---|
445 | AP_weights weights(&filter); |
---|
446 | agt->init(new AliView(GLOBAL_gb_main, filter, weights, aliName)); |
---|
447 | } |
---|
448 | } |
---|
449 | ~PARSIMONY_testenv() { |
---|
450 | TEST_EXPECT_EQUAL(ap_main, &apMain); |
---|
451 | ap_main = NULL; |
---|
452 | |
---|
453 | delete agt; |
---|
454 | GB_close(GLOBAL_gb_main); |
---|
455 | GLOBAL_gb_main = NULL; |
---|
456 | } |
---|
457 | |
---|
458 | GB_ERROR load_tree(const char *tree_name) { |
---|
459 | GB_transaction ta(GLOBAL_gb_main); // @@@ do inside AWT_graphic_tree::load? |
---|
460 | return agt->load(GLOBAL_gb_main, tree_name, 0, 0); |
---|
461 | } |
---|
462 | AP_tree_nlen *tree_root() { return apMain.get_root_node(); } |
---|
463 | |
---|
464 | void push() { apMain.push(); } |
---|
465 | void pop() { apMain.pop(); } |
---|
466 | }; |
---|
467 | |
---|
468 | void TEST_tree_modifications() { |
---|
469 | PARSIMONY_testenv env("TEST_trees.arb"); |
---|
470 | TEST_EXPECT_NO_ERROR(env.load_tree("tree_test")); |
---|
471 | |
---|
472 | { |
---|
473 | AP_tree_nlen *root = env.tree_root(); |
---|
474 | TEST_REJECT_NULL(root); |
---|
475 | |
---|
476 | AP_tree_edge::initialize(root); // builds edges |
---|
477 | TEST_ASSERT_VALID_TREE(root); |
---|
478 | |
---|
479 | root->compute_tree(); |
---|
480 | |
---|
481 | // first check initial state: |
---|
482 | { |
---|
483 | AP_tree_members& root_info = root->gr; |
---|
484 | |
---|
485 | TEST_EXPECT_EQUAL(root_info.grouped, 0); |
---|
486 | TEST_EXPECT_EQUAL(root_info.hidden, 0); |
---|
487 | TEST_EXPECT_EQUAL(root_info.has_marked_children, 1); |
---|
488 | TEST_EXPECT_EQUAL(root_info.leaf_sum, 15); |
---|
489 | |
---|
490 | TEST_EXPECT_SIMILAR(root_info.max_tree_depth, 1.624975, 0.000001); |
---|
491 | TEST_EXPECT_SIMILAR(root_info.min_tree_depth, 0.341681, 0.000001); |
---|
492 | |
---|
493 | GB_transaction ta(GLOBAL_gb_main); |
---|
494 | GBT_mark_all(GLOBAL_gb_main, 0); // unmark all species |
---|
495 | root->compute_tree(); |
---|
496 | TEST_EXPECT_EQUAL(root_info.has_marked_children, 0); |
---|
497 | } |
---|
498 | |
---|
499 | |
---|
500 | #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" |
---|
501 | #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" |
---|
502 | #define B2_TOP "(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010,CloPaste:0.179):0.131" |
---|
503 | #define B2_BOT "(CloPaste:0.179,(CloCarni:0.120,(CloButy2:0.009,CloButyr:0.000):0.564):0.010):0.131" |
---|
504 | |
---|
505 | |
---|
506 | #define B3_LEFT_TOP_SONS "(((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053,CelBiazo:0.059)" |
---|
507 | #define B3_TOP_SONS B3_LEFT_TOP_SONS ":0.207,CytAquat:0.711" |
---|
508 | #define B3_TOP_SONS_CCR "((CorAquat:0.187,CorGluta:0.522):0.053,CelBiazo:0.059):0.207,CytAquat:0.711" // CCR = CurCitre removed |
---|
509 | #define B3_TOP "(" B3_TOP_SONS "):0.081" |
---|
510 | #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" |
---|
511 | |
---|
512 | |
---|
513 | const char *top_topo = "((" B1_TOP "," B2_TOP "):0.081," B3_TOP ");"; |
---|
514 | const char *edge_topo = "((" B1_TOP "," B2_BOT "):0.081," B3_BOT ");"; |
---|
515 | const char *bottom_topo = "(" B3_BOT ",(" B2_BOT "," B1_BOT "):0.081);"; |
---|
516 | |
---|
517 | const char *radial_topo = |
---|
518 | "(((CloPaste:0.179,((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010):0.131," |
---|
519 | "((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," |
---|
520 | "((CelBiazo:0.059,((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053):0.207,CytAquat:0.711):0.081);"; |
---|
521 | const char *radial_topo2 = |
---|
522 | "(((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," |
---|
523 | "(CytAquat:0.711," B3_LEFT_TOP_SONS ":0.207):0.081);"; |
---|
524 | |
---|
525 | // expect that no mode reproduces another mode: |
---|
526 | TEST_EXPECT_DIFFERENT(top_topo, edge_topo); |
---|
527 | TEST_EXPECT_DIFFERENT(top_topo, bottom_topo); |
---|
528 | TEST_EXPECT_DIFFERENT(top_topo, radial_topo); |
---|
529 | TEST_EXPECT_DIFFERENT(top_topo, radial_topo2); |
---|
530 | TEST_EXPECT_DIFFERENT(edge_topo, bottom_topo); |
---|
531 | TEST_EXPECT_DIFFERENT(edge_topo, radial_topo); |
---|
532 | TEST_EXPECT_DIFFERENT(edge_topo, radial_topo2); |
---|
533 | TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo); |
---|
534 | TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo2); |
---|
535 | TEST_EXPECT_DIFFERENT(radial_topo, radial_topo2); |
---|
536 | |
---|
537 | env.push(); // 1st stack level (=top_topo) |
---|
538 | |
---|
539 | TEST_ASSERT_VALID_TREE(root); |
---|
540 | |
---|
541 | TEST_EXPECT_NEWICK(nLENGTH, root, top_topo); |
---|
542 | // test reorder_tree: |
---|
543 | root->reorder_tree(BIG_BRANCHES_TO_EDGE); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo); env.push(); // 2nd stack level (=edge_topo) |
---|
544 | root->reorder_tree(BIG_BRANCHES_TO_BOTTOM); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo); env.push(); // 3rd stack level (=bottom_topo) |
---|
545 | root->reorder_tree(BIG_BRANCHES_TO_CENTER); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo); |
---|
546 | root->reorder_tree(BIG_BRANCHES_ALTERNATING); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo2); |
---|
547 | root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo); |
---|
548 | |
---|
549 | TEST_ASSERT_VALID_TREE(root); |
---|
550 | |
---|
551 | // test set root: |
---|
552 | AP_tree_nlen *CloTyrob = root->findLeafNamed("CloTyrob"); |
---|
553 | TEST_REJECT_NULL(CloTyrob); |
---|
554 | |
---|
555 | ARB_edge rootEdge(root->get_leftson(), root->get_rightson()); |
---|
556 | CloTyrob->set_root(); |
---|
557 | |
---|
558 | TEST_ASSERT_VALID_TREE(root); |
---|
559 | |
---|
560 | const char *rootAtCloTyrob_topo = |
---|
561 | "(CloTyrob:0.004," |
---|
562 | "(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017," |
---|
563 | "((((" B3_TOP_SONS "):0.162," B2_TOP "):0.124,CloBifer:0.388):0.057,CloInnoc:0.371):0.274):0.004);"; |
---|
564 | |
---|
565 | TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo); |
---|
566 | env.push(); // 4th stack level (=rootAtCloTyrob_topo) |
---|
567 | |
---|
568 | TEST_ASSERT_VALID_TREE(root); |
---|
569 | |
---|
570 | AP_tree_nlen *CelBiazoFather = root->findLeafNamed("CelBiazo")->get_father(); |
---|
571 | TEST_REJECT_NULL(CelBiazoFather); |
---|
572 | CelBiazoFather->set_root(); |
---|
573 | |
---|
574 | const char *rootAtCelBiazoFather_topo = "(" B3_LEFT_TOP_SONS ":0.104,((" B1_TOP "," B2_TOP "):0.162,CytAquat:0.711):0.104);"; |
---|
575 | TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCelBiazoFather_topo); |
---|
576 | |
---|
577 | TEST_ASSERT_VALID_TREE(root); |
---|
578 | |
---|
579 | ARB_edge oldRootEdge(rootEdge.source(), rootEdge.dest()); |
---|
580 | DOWNCAST(AP_tree_nlen*,oldRootEdge.son())->set_root(); |
---|
581 | |
---|
582 | const char *rootSetBack_topo = top_topo; |
---|
583 | TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo); |
---|
584 | env.push(); // 5th stack level (=rootSetBack_topo) |
---|
585 | |
---|
586 | TEST_ASSERT_VALID_TREE(root); |
---|
587 | |
---|
588 | // test remove: |
---|
589 | AP_tree_nlen *CurCitre = root->findLeafNamed("CurCitre"); |
---|
590 | TEST_REJECT_NULL(CurCitre); |
---|
591 | TEST_REJECT_NULL(CurCitre->get_father()); |
---|
592 | |
---|
593 | CurCitre->remove(); |
---|
594 | const char *CurCitre_removed_topo = "((" B1_TOP "," B2_TOP "):0.081,(" B3_TOP_SONS_CCR "):0.081);"; |
---|
595 | // ------------------------------------------------------------------- ^^^ = B3_TOP_SONS minus CurCitre |
---|
596 | TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo); |
---|
597 | |
---|
598 | TEST_ASSERT_VALID_TREE(root); |
---|
599 | TEST_ASSERT_VALID_TREE(CurCitre); |
---|
600 | |
---|
601 | TEST_EXPECT_EQUAL(root->gr.leaf_sum, 15); // out of date |
---|
602 | root->compute_tree(); |
---|
603 | TEST_EXPECT_EQUAL(root->gr.leaf_sum, 14); |
---|
604 | |
---|
605 | env.push(); // 6th stack level (=CurCitre_removed_topo) |
---|
606 | |
---|
607 | TEST_ASSERT_VALID_TREE(root); |
---|
608 | |
---|
609 | // test insert: |
---|
610 | AP_tree_nlen *CloCarni = root->findLeafNamed("CloCarni"); |
---|
611 | TEST_REJECT_NULL(CloCarni); |
---|
612 | CurCitre->insert(CloCarni); // this creates two extra edges (not destroyed by destroy() below) and one extra node |
---|
613 | |
---|
614 | 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);"; |
---|
615 | TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_inserted_topo); |
---|
616 | |
---|
617 | AP_tree_nlen *node_del_manually = CurCitre->get_father(); |
---|
618 | AP_tree_edge *edge1_del_manually = CurCitre->edgeTo(node_del_manually); |
---|
619 | AP_tree_edge *edge2_del_manually = CurCitre->get_brother()->edgeTo(node_del_manually); |
---|
620 | |
---|
621 | TEST_ASSERT_VALID_TREE(root); |
---|
622 | |
---|
623 | // now check pops: |
---|
624 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo); |
---|
625 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo); |
---|
626 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo); |
---|
627 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo); |
---|
628 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo); |
---|
629 | env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo); |
---|
630 | |
---|
631 | TEST_ASSERT_VALID_TREE(root); |
---|
632 | |
---|
633 | AP_tree_edge::destroy(root); |
---|
634 | |
---|
635 | // delete memory allocated by insert() above and lost due to pop()s |
---|
636 | delete edge1_del_manually; |
---|
637 | delete edge2_del_manually; |
---|
638 | |
---|
639 | node_del_manually->forget_origin(); |
---|
640 | node_del_manually->father = NULL; |
---|
641 | node_del_manually->leftson = NULL; |
---|
642 | node_del_manually->rightson = NULL; |
---|
643 | delete node_del_manually; |
---|
644 | } |
---|
645 | } |
---|
646 | |
---|
647 | // @@@ Tests wanted: |
---|
648 | // - NNI |
---|
649 | // - tree optimize |
---|
650 | // - ... |
---|
651 | |
---|
652 | void TEST_calc_bootstraps() { |
---|
653 | PARSIMONY_testenv env("TEST_trees.arb", "ali_5s"); |
---|
654 | TEST_EXPECT_NO_ERROR(env.load_tree("tree_test")); |
---|
655 | |
---|
656 | 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%');"; |
---|
657 | 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%');"; |
---|
658 | 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%');"; |
---|
659 | |
---|
660 | { |
---|
661 | AP_tree_nlen *root = env.tree_root(); |
---|
662 | TEST_REJECT_NULL(root); |
---|
663 | |
---|
664 | AP_tree_edge::initialize(root); // builds edges |
---|
665 | TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that) |
---|
666 | |
---|
667 | AP_tree_edge *root_edge = rootEdge(); |
---|
668 | TEST_REJECT_NULL(root_edge); |
---|
669 | |
---|
670 | root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_origi_topo); |
---|
671 | |
---|
672 | 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); |
---|
673 | 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); |
---|
674 | |
---|
675 | TEST_EXPECT_EQUAL(env.tree_root(), root); |
---|
676 | AP_tree_edge::destroy(root); |
---|
677 | } |
---|
678 | } |
---|
679 | |
---|
680 | void TEST_tree_remove_add_all() { |
---|
681 | // reproduces crash as described in #527 |
---|
682 | PARSIMONY_testenv env("TEST_trees.arb", "ali_5s"); |
---|
683 | TEST_EXPECT_NO_ERROR(env.load_tree("tree_nj")); |
---|
684 | |
---|
685 | const int LEAFS = 6; |
---|
686 | AP_tree_nlen *leaf[LEAFS]; |
---|
687 | const char *name[LEAFS] = { |
---|
688 | "CloButy2", |
---|
689 | "CloButyr", |
---|
690 | "CytAquat", |
---|
691 | "CorAquat", |
---|
692 | "CurCitre", |
---|
693 | "CorGluta", |
---|
694 | }; |
---|
695 | |
---|
696 | AP_tree_nlen *root = env.tree_root(); |
---|
697 | TEST_REJECT_NULL(root); |
---|
698 | |
---|
699 | for (int i = 0; i<LEAFS; ++i) { |
---|
700 | leaf[i] = root->findLeafNamed(name[i]); |
---|
701 | TEST_REJECT_NULL(leaf[i]); |
---|
702 | } |
---|
703 | |
---|
704 | AP_tree_edge::initialize(root); // builds edges |
---|
705 | TEST_EXPECT_EQUAL(root, rootNode()); // need tree-access via global 'ap_main' (too much code is based on that) |
---|
706 | |
---|
707 | AP_tree_root *troot = leaf[0]->get_tree_root(); |
---|
708 | TEST_REJECT_NULL(troot); |
---|
709 | |
---|
710 | // Note: following loop leaks father nodes and edges |
---|
711 | // suppressed in valgrind via ../SOURCE_TOOLS/arb.supp@TEST_tree_remove_add_all |
---|
712 | for (int i = 0; i<LEAFS-1; ++i) { // removing the second to last leaf, "removes" both remaining leafs |
---|
713 | TEST_ASSERT_VALID_TREE(root); |
---|
714 | leaf[i]->remove(); |
---|
715 | TEST_ASSERT_VALID_TREE(leaf[i]); |
---|
716 | } |
---|
717 | leaf[LEAFS-1]->father = NULL; // correct final leaf (not removed regularily) |
---|
718 | |
---|
719 | leaf[0]->initial_insert(leaf[1], troot); |
---|
720 | for (int i = 2; i<LEAFS; ++i) { |
---|
721 | TEST_ASSERT_VALID_TREE(leaf[i-1]); |
---|
722 | TEST_ASSERT_VALID_TREE(leaf[i]); |
---|
723 | leaf[i]->insert(leaf[i-1]); |
---|
724 | } |
---|
725 | } |
---|
726 | |
---|
727 | #endif // UNIT_TESTS |
---|
728 | |
---|
729 | // -------------------------------------------------------------------------------- |
---|