Context Navigation

source: tags/ms_r16q3/PARSIMONY/PARS_dtree.cxx

Visit:

Last change on this file was 15064, checked in by westram, 8 years ago
remove parameter `base_gc` from `AW_manage_GC` (always zero)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.6 KB

Line
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 "PerfMeter.h"
12	#include "pars_dtree.hxx"
13	#include "pars_main.hxx"
14	#include "pars_awars.h"
15	#include "ap_tree_nlen.hxx"
16	#include "ap_main.hxx"
17
18	#include <AP_TreeColors.hxx>
19	#include <AP_seq_dna.hxx>
20	#include <AP_seq_protein.hxx>
21	#include <AP_filter.hxx>
22
23	#include <TreeCallbacks.hxx>
24
25	#include <ColumnStat.hxx>
26	#include <awt_sel_boxes.hxx>
27	#include <awt_filter.hxx>
28
29	#include <gui_aliview.hxx>
30
31	#include <aw_preset.hxx>
32	#include <aw_awar.hxx>
33	#include <aw_msg.hxx>
34	#include <arb_progress.h>
35	#include <aw_root.hxx>
36	#include <aw_question.hxx>
37
38	static void AWT_graphic_parsimony_root_changed(void cd, AP_tree old, AP_tree *newroot) {
39	AWT_graphic_tree agt = (AWT_graphic_tree)cd;
40	UNCOVERED();
41
42	if (old == agt->get_logical_root()) agt->set_logical_root_to(newroot);
43	}
44
45	static AliView pars_generate_aliview(WeightedFilter pars_weighted_filter) {
46	GBDATA *gb_main = pars_weighted_filter->get_gb_main();
47	char *ali_name;
48	{
49	GB_transaction ta(gb_main);
50	ali_name = GBT_read_string(gb_main, AWAR_ALIGNMENT);
51	}
52	GB_ERROR error = NULL;
53	AliView *aliview = pars_weighted_filter->create_aliview(ali_name, error);
54	if (!aliview) aw_popup_exit(error);
55	free(ali_name);
56	return aliview;
57	}
58
59	void PARS_tree_init(AWT_graphic_parsimony *agt) {
60	ap_assert(agt->get_root_node());
61	ap_assert(agt == ap_main->get_graphic_tree());
62
63	GBDATA *gb_main = ap_main->get_gb_main();
64	GB_transaction ta(gb_main);
65
66	const char *use = ap_main->get_aliname();
67	long ali_len = GBT_get_alignment_len(gb_main, use);
68	if (ali_len <= 1) {
69	aw_popup_exit("No valid alignment selected! Try again");
70	}
71
72	agt->get_tree_root()->set_root_changed_callback(AWT_graphic_parsimony_root_changed, agt);
73	}
74
75	QuadraticThreshold::QuadraticThreshold(KL_DYNAMIC_THRESHOLD_TYPE type, double startx, double maxy, double maxx, double maxDepth, AP_FLOAT pars_start) {
76	// set a, b, and c for quadratic equation y = ax^2 + bx + c
77	switch (type) {
78	default:
79	ap_assert(0);
80	// fall-through
81	case AP_QUADRAT_START:
82	c = startx;
83	a = (startx - maxy) / (maxx * maxx);
84	b = -2.0 * a * maxx;
85	break;
86
87	case AP_QUADRAT_MAX: // unused (experimental)
88	a = - maxy / ((maxDepth - maxx) * (maxDepth - maxx));
89	b = -2.0 * a * maxx;
90	c = maxy + a * maxx * maxx;
91	break;
92	}
93	c += pars_start;
94	}
95
96	void ArbParsimony::kernighan_optimize_tree(AP_tree_nlen at, const KL_Settings& settings, const AP_FLOAT pars_global_start, bool dumpPerf) {
97	AP_tree_nlen *oldrootleft = get_root_node()->get_leftson();
98	AP_tree_nlen *oldrootright = get_root_node()->get_rightson();
99	AP_FLOAT pars_curr = get_root_node()->costs();
100	const AP_FLOAT pars_org = pars_curr;
101
102	OptiPerfMeter performance("KL-recursion", pars_curr);
103
104	// setup KL recursion parameters:
105	KL_params KL;
106	{
107	KL.max_rec_depth = settings.maxdepth; ap_assert(KL.max_rec_depth>0);
108	KL.inc_rec_depth = settings.incdepth;
109	KL.thresFunctor = QuadraticThreshold(settings.Dynamic.type, settings.Dynamic.start, settings.Dynamic.maxy, settings.Dynamic.maxx, KL.max_rec_depth, pars_curr);
110	KL.rec_type = KL_RECURSION_TYPE((settings.Dynamic.enabledAP_DYNAMIK)\|(settings.Static.enabledAP_STATIC));
111
112	for (int x = 0; x<CUSTOM_DEPTHS; ++x) {
113	KL.rec_width[x] = settings.Static.depth[x];
114	}
115	KL.stopAtFoldedGroups = settings.whichEdges&SKIP_FOLDED_EDGES;
116	}
117
118	AP_tree_edge *startEdge = NULL;
119	AP_tree_nlen *skipNode = NULL;
120	bool skipStartEdge = true;
121
122	if (!at->father) {
123	startEdge = rootEdge();
124	skipStartEdge = false;
125	}
126	else if (!at->father->father) {
127	startEdge = rootEdge();
128	skipNode = startEdge->otherNode(at);
129	}
130	else {
131	skipNode = at->get_father();
132	startEdge = at->edgeTo(skipNode);
133	}
134	ap_assert(startEdge);
135
136	EdgeChain chain(startEdge, EdgeSpec(SKIP_LEAF_EDGES\|settings.whichEdges), true, skipNode, skipStartEdge);
137	arb_progress progress(chain.size());
138
139	if (pars_global_start) {
140	progress.subtitle(GBS_global_string("best=%.1f (gain=%.1f)", pars_curr, *pars_global_start-pars_curr));
141	}
142	else {
143	progress.subtitle(GBS_global_string("best=%.1f", pars_curr));
144	}
145
146	if (skipStartEdge) startEdge->set_visited(true); // avoid traversal beyond startEdge
147
148	while (chain && !progress.aborted()) {
149	AP_tree_edge edge = chain; ++chain;
150
151	ap_assert(!edge->is_leaf_edge());
152	ap_assert(implicated(KL.stopAtFoldedGroups, !edge->next_to_folded_group()));
153
154	ap_main->remember();
155
156	bool better_tree_found = edge->kl_rec(KL, 0, pars_curr);
157
158	if (better_tree_found) {
159	ap_main->accept();
160	AP_FLOAT pars_new = get_root_node()->costs();
161	KL.thresFunctor.change_parsimony_start(pars_new-pars_curr);
162	pars_curr = pars_new;
163	if (pars_global_start) {
164	progress.subtitle(GBS_global_string("best=%.1f (gain=%.1f, KL=%.1f)", pars_curr, *pars_global_start-pars_curr, pars_org-pars_curr));
165	}
166	else {
167	progress.subtitle(GBS_global_string("best=%.1f (gain=%.1f)", pars_curr, pars_org-pars_curr));
168	}
169	}
170	else {
171	ap_main->revert();
172	}
173	progress.inc();
174	}
175
176	if (skipStartEdge) startEdge->set_visited(false); // reallow traversal beyond startEdge
177
178	if (dumpPerf) performance.dump(stdout, pars_curr);
179
180	if (oldrootleft->father == oldrootright) {
181	oldrootleft->set_root();
182	}
183	else {
184	oldrootright->set_root();
185	}
186	}
187
188
189
190	void ArbParsimony::optimize_tree(AP_tree_nlen *at, const KL_Settings& settings, arb_progress& progress) {
191	AP_tree_nlen *oldrootleft = get_root_node()->get_leftson();
192	AP_tree_nlen *oldrootright = get_root_node()->get_rightson();
193	const AP_FLOAT org_pars = get_root_node()->costs();
194	AP_FLOAT prev_pars = org_pars;
195
196	OptiPerfMeter overallPerf("global optimization", org_pars);
197
198	progress.subtitle(GBS_global_string("best=%.1f", org_pars));
199
200	// define available heuristics
201	enum Heuristic {
202	RECURSIVE_NNI,
203	CUSTOM_KL,
204
205	NO_FURTHER_HEURISTIC,
206	HEURISTIC_COUNT = NO_FURTHER_HEURISTIC,
207	} heuristic = RECURSIVE_NNI;
208
209	struct H_Settings {
210	const char *name; // name shown in OptiPerfMeter
211	const KL_Settings *kl; // ==NULL -> NNI; else KL with these settings
212	bool repeat; // retry same heuristic when tree improved
213	Heuristic onImprove; // continue with this heuristic if improved (repeated or not)
214	Heuristic onFailure; // continue with this heuristic if NOT improved
215	} heuristic_setting[HEURISTIC_COUNT] = {
216	{ "recursive NNIs", NULL, true, CUSTOM_KL, CUSTOM_KL },
217	{ "KL-optimizer", &settings, false, RECURSIVE_NNI, NO_FURTHER_HEURISTIC },
218	};
219
220	AP_FLOAT heu_start_pars = prev_pars;
221	OptiPerfMeter *heuPerf = NULL;
222
223	#if defined(ASSERTION_USED)
224	bool at_is_root = at == rootNode();
225	#endif
226
227	while (heuristic != NO_FURTHER_HEURISTIC && !progress.aborted()) {
228	const H_Settings& hset = heuristic_setting[heuristic];
229	if (!heuPerf) {
230	ap_assert(heu_start_pars == prev_pars);
231	heuPerf = new OptiPerfMeter(hset.name, heu_start_pars);
232	}
233
234	AP_FLOAT this_pars = -1;
235	if (hset.kl) {
236	kernighan_optimize_tree(at, *hset.kl, &org_pars, false);
237	this_pars = get_root_node()->costs();
238	}
239	else {
240	this_pars = at->nn_interchange_rec(settings.whichEdges, AP_BL_NNI_ONLY);
241	}
242	ap_assert(this_pars>=0); // ensure this_pars was set
243	ap_assert(this_pars<=prev_pars); // otherwise heuristic worsened the tree
244
245	ap_assert(at_is_root == (at == rootNode()));
246
247	bool dumpOverall = false;
248	Heuristic nextHeuristic = heuristic;
249	if (this_pars<prev_pars) { // found better tree
250	prev_pars = this_pars;
251	progress.subtitle(GBS_global_string("best=%.1f (gain=%.1f)", this_pars, org_pars-this_pars));
252	if (!hset.repeat) {
253	nextHeuristic = hset.onImprove;
254	dumpOverall = heuristic == CUSTOM_KL;
255	}
256	}
257	else { // last step did not find better tree
258	nextHeuristic = this_pars<heu_start_pars ? hset.onImprove : hset.onFailure;
259	}
260
261	if (nextHeuristic != heuristic) {
262	heuristic = nextHeuristic;
263	heu_start_pars = this_pars;
264
265	heuPerf->dump(stdout, this_pars);
266	delete heuPerf; heuPerf = NULL;
267	}
268	if (dumpOverall) overallPerf.dumpCustom(stdout, this_pars, "overall (so far)");
269	}
270
271	if (oldrootleft->father == oldrootright) {
272	oldrootleft->set_root();
273	}
274	else {
275	oldrootright->set_root();
276	}
277
278	overallPerf.dump(stdout, prev_pars);
279	}
280
281	AWT_graphic_parsimony::AWT_graphic_parsimony(ArbParsimony& parsimony_, GBDATA *gb_main_, AD_map_viewer_cb map_viewer_cb_)
282	: AWT_graphic_tree(AW_root::SINGLETON, gb_main_, map_viewer_cb_),
283	parsimony(parsimony_)
284	{}
285
286	AP_tree_root AWT_graphic_parsimony::create_tree_root(AliView aliview, AP_sequence *seq_prototype, bool insert_delete_cbs) {
287	return new AP_pars_root(aliview, seq_prototype, insert_delete_cbs);
288	}
289
290
291	void ArbParsimony::generate_tree(WeightedFilter *pars_weighted_filter) {
292	AliView *aliview = pars_generate_aliview(pars_weighted_filter);
293	AP_sequence *seq_templ = 0;
294
295	GBDATA *gb_main = aliview->get_gb_main();
296	{
297	GB_transaction ta(gb_main);
298	bool is_aa = GBT_is_alignment_protein(gb_main, aliview->get_aliname());
299
300	if (is_aa) seq_templ = new AP_sequence_protein(aliview);
301	else seq_templ = new AP_sequence_parsimony(aliview);
302	}
303
304	AWT_graphic_parsimony new_tree = new AWT_graphic_parsimony(this, aliview->get_gb_main(), PARS_map_viewer);
305	new_tree->init(aliview, seq_templ, true, false);
306	set_tree(new_tree);
307	}
308
309	AW_gc_manager AWT_graphic_parsimony::init_devices(AW_window aww, AW_device device, AWT_canvas ntw) {
310	AW_init_color_group_defaults("arb_ntree");
311
312	AW_gc_manager *gc_manager =
313	AW_manage_GC(aww,
314	ntw->get_gc_base_name(),
315	device, AWT_GC_MAX, AW_GCM_DATA_AREA,
316	makeGcChangedCallback(TREE_GC_changed_cb, ntw),
317	"#AAAA55",
318
319	// Important note :
320	// Many gc indices are shared between ABR_NTREE and ARB_PARSIMONY
321	// e.g. the tree drawing routines use same gc's for drawing both trees
322	// (keep in sync with ../SL/TREEDISP/TreeDisplay.cxx@init_devices)
323
324	"Cursor$#FFFFFF",
325	"Branch remarks$#DBE994",
326	"+-Bootstrap$#DBE994", "-B.(limited)$white",
327	"!1", // reserve GC-number used for "IRS group box" in arb_ntree
328	"Marked$#FFFF00",
329	"Some marked$#eeee88",
330	"Not marked$black",
331	"Zombies etc.$#cc5924",
332
333	"!14", // reserve 14 GC-numbers which are used for probe colors in ARB_NTREE
334	// (namely 'None', 'All' and 'P1' up to 'P12')
335
336	"&color_groups", // use color groups
337
338	NULL);
339	return gc_manager;
340	}
341
342	void AWT_graphic_parsimony::show(AW_device *device) {
343	AP_tree_nlen *root_node = parsimony.get_root_node();
344	AW_awar *awar_pars = aw_root->awar(AWAR_PARSIMONY);
345	AW_awar *awar_best = aw_root->awar(AWAR_BEST_PARSIMONY);
346
347	long parsval = root_node ? root_node->costs() : 0;
348	awar_pars->write_int(parsval);
349
350	long best_pars = awar_best->read_int();
351	if (parsval < best_pars \|\| 0==best_pars) awar_best->write_int(parsval);
352
353	AWT_graphic_tree::show(device);
354	}
355
356	void AWT_graphic_parsimony::handle_command(AW_device *device, AWT_graphic_event& event) {
357	ClickedTarget clicked(this, event.best_click());
358	bool recalc_branchlengths_on_structure_change = true;
359
360	switch (event.cmd()) {
361	case AWT_MODE_NNI:
362	case AWT_MODE_KERNINGHAN:
363	case AWT_MODE_OPTIMIZE: {
364	const char *what = NULL;
365	const char *where = NULL;
366
367	switch (event.cmd()) {
368	case AWT_MODE_NNI: what = "Recursive NNI on"; break;
369	case AWT_MODE_KERNINGHAN: what = "K.L. optimize"; break;
370	case AWT_MODE_OPTIMIZE: what = "Global optimize"; break;
371	default: break;
372	}
373
374	AP_tree_nlen *startNode = NULL;
375	bool repeatOpti = true;
376
377	if (event.type() == AW_Mouse_Press) {
378	switch (event.button()) {
379	case AW_BUTTON_LEFT:
380	repeatOpti = false;
381	// fall-through
382	case AW_BUTTON_RIGHT:
383	startNode = DOWNCAST(AP_tree_nlen*, clicked.node());
384	where = (startNode == get_root_node()) ? "tree" : "subtree";
385	break;
386
387	default:
388	break;
389	}
390	}
391
392	if (what && where) {
393	arb_progress progress(GBS_global_string("%s %s", what, where));
394
395	AP_FLOAT start_pars = get_root_node()->costs();
396	AP_FLOAT curr_pars = start_pars;
397
398	KL_Settings KL(aw_root);
399
400	do {
401	AP_FLOAT prev_pars = curr_pars;
402
403	switch (event.cmd()) {
404	case AWT_MODE_NNI:
405	startNode->nn_interchange_rec(KL.whichEdges, AP_BL_NNI_ONLY);
406	break;
407	case AWT_MODE_KERNINGHAN:
408	parsimony.kernighan_optimize_tree(startNode, KL, &start_pars, true);
409	break;
410	case AWT_MODE_OPTIMIZE:
411	parsimony.optimize_tree(startNode, KL, progress);
412	repeatOpti = false; // never loop here (optimize_tree already loops until no further improvement)
413	break;
414	default:
415	repeatOpti = false;
416	break;
417	}
418
419	curr_pars = get_root_node()->costs();
420	repeatOpti = repeatOpti && curr_pars<prev_pars;
421	} while (repeatOpti);
422
423	exports.save = 1;
424	ASSERT_VALID_TREE(get_root_node());
425	}
426	break;
427	}
428
429	default:
430	recalc_branchlengths_on_structure_change = false;
431	// fall-through (unlisted modes trigger branchlength calculation internally when needed)
432	case AWT_MODE_MOVE:
433	AWT_graphic_tree::handle_command(device, event);
434	break;
435	}
436
437	if (exports.save == 1 && recalc_branchlengths_on_structure_change) {
438	arb_progress progress("Recalculating branch lengths");
439	rootEdge()->calc_branchlengths();
440	reorder_tree(BIG_BRANCHES_TO_TOP); // beautify after recalc_branch_lengths
441	}
442	}
443
444
445	// --------------------------------------------------------------------------------
446
447	#ifdef UNIT_TESTS
448	#include <arb_diff.h>
449	#include <test_unit.h>
450	#include "test_env.h"
451
452	template<typename SEQTYPE>
453	PARSIMONY_testenv<SEQTYPE>::PARSIMONY_testenv(const char dbname, const char aliName)
454	: parsimony(),
455	klSettings(NULL)
456	{
457	common_init(dbname);
458	GBDATA *gb_main = ap_main->get_gb_main();
459	GB_transaction ta(gb_main);
460	size_t aliLength = GBT_get_alignment_len(gb_main, aliName);
461
462	klSettings = new KL_Settings();
463
464	AP_filter filter(aliLength);
465	if (!filter.is_invalid()) {
466	AP_weights weights(&filter);
467	agt->init(new AliView(gb_main, filter, weights, aliName));
468	}
469	}
470
471	template PARSIMONY_testenv<AP_sequence_protein>::PARSIMONY_testenv(const char dbname, const char aliName); // explicit instanciation (otherwise link error in unittest)
472	template PARSIMONY_testenv<AP_sequence_parsimony>::PARSIMONY_testenv(const char dbname, const char aliName); // explicit instanciation (as above, but only happens with gcc 4.6.3/NDEBUG)
473
474
475	void TEST_basic_tree_modifications() {
476	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb");
477	TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
478
479	{
480	AP_tree_nlen *root = env.root_node();
481
482	// first check initial state:
483	{
484	AP_tree_members& root_info = root->gr;
485
486	TEST_EXPECT_EQUAL(root_info.grouped, 0);
487	TEST_EXPECT_EQUAL(root_info.hidden, 0);
488	TEST_EXPECT_EQUAL(root_info.mark_sum, 6);
489	TEST_EXPECT_EQUAL(root_info.leaf_sum, 15);
490
491	TEST_EXPECT_SIMILAR(root_info.max_tree_depth, 1.624975, 0.000001);
492	TEST_EXPECT_SIMILAR(root_info.min_tree_depth, 0.341681, 0.000001);
493
494	GB_transaction ta(env.gbmain());
495	GBT_mark_all(env.gbmain(), 0); // unmark all species
496	root->compute_tree();
497	TEST_EXPECT_EQUAL(root_info.mark_sum, 0);
498	}
499
500
501	#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"
502	#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"
503	#define B2_TOP "(((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010,CloPaste:0.179):0.131"
504	#define B2_BOT "(CloPaste:0.179,(CloCarni:0.120,(CloButy2:0.009,CloButyr:0.000):0.564):0.010):0.131"
505
506
507	#define B3_LEFT_TOP_SONS "(((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053,CelBiazo:0.059)"
508	#define B3_TOP_SONS B3_LEFT_TOP_SONS ":0.207,CytAquat:0.711"
509	#define B3_TOP_SONS_CCR "((CorAquat:0.187,CorGluta:0.522):0.053,CelBiazo:0.059):0.207,CytAquat:0.711" // CCR = CurCitre removed
510	#define B3_TOP "(" B3_TOP_SONS "):0.081"
511	#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"
512
513
514	const char *top_topo = "((" B1_TOP "," B2_TOP "):0.081," B3_TOP ");";
515	const char *edge_topo = "((" B1_TOP "," B2_BOT "):0.081," B3_BOT ");";
516	const char *bottom_topo = "(" B3_BOT ",(" B2_BOT "," B1_BOT "):0.081);";
517
518	const char *radial_topo =
519	"(((CloPaste:0.179,((CloButy2:0.009,CloButyr:0.000):0.564,CloCarni:0.120):0.010):0.131,"
520	"((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,"
521	"((CelBiazo:0.059,((CorAquat:0.084,CurCitre:0.058):0.103,CorGluta:0.522):0.053):0.207,CytAquat:0.711):0.081);";
522	const char *radial_topo2 =
523	"(((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,"
524	"(CytAquat:0.711," B3_LEFT_TOP_SONS ":0.207):0.081);";
525
526	// expect that no mode reproduces another mode:
527	TEST_EXPECT_DIFFERENT(top_topo, edge_topo);
528	TEST_EXPECT_DIFFERENT(top_topo, bottom_topo);
529	TEST_EXPECT_DIFFERENT(top_topo, radial_topo);
530	TEST_EXPECT_DIFFERENT(top_topo, radial_topo2);
531	TEST_EXPECT_DIFFERENT(edge_topo, bottom_topo);
532	TEST_EXPECT_DIFFERENT(edge_topo, radial_topo);
533	TEST_EXPECT_DIFFERENT(edge_topo, radial_topo2);
534	TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo);
535	TEST_EXPECT_DIFFERENT(bottom_topo, radial_topo2);
536	TEST_EXPECT_DIFFERENT(radial_topo, radial_topo2);
537
538	env.push(); // 1st stack level (=top_topo)
539
540	TEST_EXPECT_VALID_TREE(root);
541
542	TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
543	// test reorder_tree:
544	root->reorder_tree(BIG_BRANCHES_TO_EDGE); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo); env.push(); // 2nd stack level (=edge_topo)
545	root->reorder_tree(BIG_BRANCHES_TO_BOTTOM); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo); env.push(); // 3rd stack level (=bottom_topo)
546	root->reorder_tree(BIG_BRANCHES_TO_CENTER); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo);
547	root->reorder_tree(BIG_BRANCHES_ALTERNATING); TEST_EXPECT_NEWICK(nLENGTH, root, radial_topo2);
548	root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo);
549
550	TEST_EXPECT_VALID_TREE(root);
551
552	// test set root:
553	AP_tree_nlen *CloTyrob = root->findLeafNamed("CloTyrob");
554	TEST_REJECT_NULL(CloTyrob);
555
556	ARB_edge rootEdge(root->get_leftson(), root->get_rightson());
557	CloTyrob->set_root();
558
559	TEST_EXPECT_VALID_TREE(root);
560
561	const char *rootAtCloTyrob_topo =
562	"(CloTyrob:0.004,"
563	"(((CloTyro3:1.046,CloTyro4:0.061):0.026,CloTyro2:0.017):0.017,"
564	"((((" B3_TOP_SONS "):0.162," B2_TOP "):0.124,CloBifer:0.388):0.057,CloInnoc:0.371):0.274):0.004);";
565
566	TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo);
567	env.push(); // 4th stack level (=rootAtCloTyrob_topo)
568
569	TEST_EXPECT_VALID_TREE(root);
570
571	AP_tree_nlen *CelBiazoFather = root->findLeafNamed("CelBiazo")->get_father();
572	TEST_REJECT_NULL(CelBiazoFather);
573	CelBiazoFather->set_root();
574
575	const char *rootAtCelBiazoFather_topo = "(" B3_LEFT_TOP_SONS ":0.104,((" B1_TOP "," B2_TOP "):0.162,CytAquat:0.711):0.104);";
576	TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCelBiazoFather_topo);
577
578	TEST_EXPECT_VALID_TREE(root);
579
580	ARB_edge oldRootEdge(rootEdge.source(), rootEdge.dest());
581	DOWNCAST(AP_tree_nlen*,oldRootEdge.son())->set_root();
582
583	const char *rootSetBack_topo = top_topo;
584	TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo);
585	env.push(); // 5th stack level (=rootSetBack_topo)
586
587	TEST_EXPECT_VALID_TREE(root);
588
589	// test remove:
590	AP_tree_nlen *CurCitre = root->findLeafNamed("CurCitre");
591	TEST_REJECT_NULL(CurCitre);
592	TEST_REJECT_NULL(CurCitre->get_father());
593
594	CurCitre->REMOVE();
595	const char *CurCitre_removed_topo = "((" B1_TOP "," B2_TOP "):0.081,(" B3_TOP_SONS_CCR "):0.081);";
596	// ------------------------------------------------------------------- ^^^ = B3_TOP_SONS minus CurCitre
597	TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo);
598
599	TEST_EXPECT_VALID_TREE(root);
600	TEST_EXPECT_VALID_TREE(CurCitre);
601
602	TEST_EXPECT_EQUAL(root->gr.leaf_sum, 15); // out of date
603	root->compute_tree();
604	TEST_EXPECT_EQUAL(root->gr.leaf_sum, 14);
605
606	env.push(); // 6th stack level (=CurCitre_removed_topo)
607
608	TEST_EXPECT_VALID_TREE(root);
609
610	// test insert:
611	AP_tree_nlen *CloCarni = root->findLeafNamed("CloCarni");
612	TEST_REJECT_NULL(CloCarni);
613	CurCitre->insert(CloCarni); // this creates two extra edges (not destroyed by destroy() below) and one extra node
614
615	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);";
616	TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_inserted_topo);
617
618	TEST_EXPECT_VALID_TREE(root);
619
620	// now check pops:
621	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, CurCitre_removed_topo); TEST_EXPECT_VALID_TREE(root);
622	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootSetBack_topo); TEST_EXPECT_VALID_TREE(root);
623	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, rootAtCloTyrob_topo); TEST_EXPECT_VALID_TREE(root);
624	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, bottom_topo); TEST_EXPECT_VALID_TREE(root);
625	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, edge_topo); TEST_EXPECT_VALID_TREE(root);
626	env.pop(); TEST_EXPECT_NEWICK(nLENGTH, root, top_topo); TEST_EXPECT_VALID_TREE(root);
627	}
628	}
629
630	void TEST_calc_bootstraps() {
631	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", "ali_5s");
632	TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
633
634	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%');";
635	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%');";
636	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%');";
637
638	{
639	AP_tree_nlen *root = env.root_node();
640	AP_tree_edge *root_edge = rootEdge();
641
642	TEST_EXPECT(root && root_edge);
643
644	root->reorder_tree(BIG_BRANCHES_TO_TOP); TEST_EXPECT_NEWICK(nREMARK, root, bs_origi_topo);
645
646	TEST_EXPECT_EQUAL(env.combines_performed(), 0);
647
648	root_edge->nni_rec(ANY_EDGE, AP_BL_MODE(AP_BL_BL_ONLY\|AP_BL_BOOTSTRAP_LIMIT), NULL, true);
649	root->reorder_tree(BIG_BRANCHES_TO_TOP);
650	TEST_EXPECT_NEWICK(nREMARK, root, bs_limit_topo);
651	TEST_EXPECT_EQUAL(env.combines_performed(), 170);
652
653	root_edge->nni_rec(ANY_EDGE, AP_BL_MODE(AP_BL_BL_ONLY\|AP_BL_BOOTSTRAP_ESTIMATE), NULL, true);
654	root->reorder_tree(BIG_BRANCHES_TO_TOP);
655	TEST_EXPECT_NEWICK(nREMARK, root, bs_estim_topo);
656	TEST_EXPECT_EQUAL(env.combines_performed(), 156);
657
658	TEST_EXPECT_EQUAL(env.root_node(), root);
659	}
660
661	}
662
663	void TEST_tree_remove_add_all() {
664	// reproduces crash as described in #527
665	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", "ali_5s");
666	TEST_EXPECT_NO_ERROR(env.load_tree("tree_nj"));
667
668	const int LEAFS = 6;
669	AP_tree_nlen *leaf[LEAFS];
670	const char *name[LEAFS] = {
671	"CloButy2",
672	"CloButyr",
673	"CytAquat",
674	"CorAquat",
675	"CurCitre",
676	"CorGluta",
677	};
678
679	AP_tree_nlen *root = env.root_node();
680
681	for (int i = 0; i<LEAFS; ++i) {
682	leaf[i] = root->findLeafNamed(name[i]);
683	TEST_REJECT_NULL(leaf[i]);
684	}
685
686	TEST_EXPECT_VALID_TREE(root);
687
688	AP_pars_root *troot = leaf[0]->get_tree_root();
689	TEST_REJECT_NULL(troot);
690
691	for (int i = 0; i<LEAFS-1; ++i) {
692	// Note: removing the second to last leaf, "removes" both remaining
693	// leafs (but only destroys their father node)
694
695	TEST_EXPECT_VALID_TREE(root);
696	leaf[i]->REMOVE();
697	TEST_EXPECT_VALID_TREE(leaf[i]);
698	}
699
700	leaf[0]->initial_insert(leaf[1], troot);
701	for (int i = 2; i<LEAFS; ++i) {
702	TEST_EXPECT_VALID_TREE(leaf[i-1]);
703	TEST_EXPECT_VALID_TREE(leaf[i]);
704	leaf[i]->insert(leaf[i-1]);
705	}
706	}
707
708	#endif // UNIT_TESTS
709
710	// --------------------------------------------------------------------------------

Note: See TracBrowser for help on using the repository browser.

Download in other formats: