Context Navigation

source: tags/ms_ra2q34/PARSIMONY/PARS_main.cxx

Visit:

Last change on this file was 17980, checked in by westram, 5 years ago
use security objects instead of manually calling GB_push/pop_my_security. eliminate GB_push/pop_my_security. fix scoping using a TA object.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 125.1 KB

Line
1	// =============================================================== //
2	// //
3	// File : PARS_main.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_main.hxx"
13	#include "pars_klprops.hxx"
14	#include "pars_awars.h"
15	#include "ap_tree_nlen.hxx"
16	#include "ap_main.hxx"
17
18	#include <ColumnStat.hxx>
19	#include <gui_aliview.hxx>
20	#include <macros.hxx>
21	#include <nds.h>
22	#include <TreeCallbacks.hxx>
23
24	#include <aw_awars.hxx>
25	#include <aw_preset.hxx>
26	#include <aw_msg.hxx>
27	#include <aw_root.hxx>
28	#include <aw_question.hxx>
29
30	#include <awt.hxx>
31	#include <awt_sel_boxes.hxx>
32	#include <awt_filter.hxx>
33	#include <awt_config_manager.hxx>
34
35	#include <arb_progress.h>
36	#include <arb_misc.h>
37	#include <arb_defs.h>
38	#include <arb_global_defs.h>
39
40	#include <ad_cb.h>
41
42	#include <list>
43	#include <map>
44
45	#if defined(DEBUG)
46	# define TESTMENU
47	#endif // DEBUG
48
49	using namespace std;
50
51	AW_HEADER_MAIN
52
53	#define AWAR_COLUMNSTAT_BASE "tmp/pars/colstat"
54	#define AWAR_COLUMNSTAT_NAME AWAR_COLUMNSTAT_BASE "/name"
55
56	#define AWT_TREE_PARS(ntw) DOWNCAST(AWT_graphic_parsimony*, (ntw)->gfx)
57
58	static ArbParsimony *GLOBAL_PARS = NULp;
59
60	inline AWT_graphic_parsimony *global_tree() { return GLOBAL_PARS->get_tree(); }
61	inline AP_pars_root *global_tree_root() { return global_tree()->get_tree_root(); }
62
63	// waaah more globals :(
64	AP_main *ap_main; // @@@ move into ArbParsimony? or eliminate ArbParsimony
65
66	void ArbParsimony::set_tree(AWT_graphic_parsimony *tree_) {
67	ap_assert(!tree); // only call once
68	tree = tree_;
69	ap_main->set_tree_root(tree);
70	}
71
72	static void set_keep_ghostnodes() {
73	// avoid that saving tree to DB does delete removed nodes
74	// (hack to fix #528)
75	// see ../ARBDB/adtree.cxx@keep_ghostnodes
76	GBDATA *gb_tree = ap_main->get_tree_root()->get_gb_tree();
77	GB_transaction ta(gb_tree);
78	GBDATA *gb_keep = GB_searchOrCreate_int(gb_tree, "keep_ghostnodes", 1);
79	ASSERT_NO_ERROR(GB_set_temporary(gb_keep));
80	}
81	static void delete_kept_ghostnodes() {
82	if (ap_main->get_graphic_tree()) {
83	GBDATA *gb_tree = ap_main->get_tree_root()->get_gb_tree();
84	GB_transaction ta(gb_tree);
85
86	GBDATA *gb_keep = GB_entry(gb_tree, "keep_ghostnodes");
87	if (gb_keep) { // e.g. wrong for quick-add species
88	GB_ERROR error = GB_delete(gb_keep);
89	if (!error) {
90	if (ap_main->get_tree_root()->was_saved()) {
91	// if tree was saved, DB may contain ghostnodes
92	// -> save again to delete them
93	error = global_tree()->save_to_DB(GB_get_root(gb_tree), NULp);
94	}
95	}
96	if (error) aw_message(error);
97	}
98	}
99	}
100
101	__ATTR__NORETURN static void pars_exit(AW_window *aww) {
102	AW_root *aw_root = aww->get_root();
103	shutdown_macro_recording(aw_root);
104
105	ap_main->accept_all();
106	delete_kept_ghostnodes();
107
108	aw_root->unlink_awars_from_DB(ap_main->get_gb_main());
109	#if defined(DEBUG)
110	AWT_browser_forget_db(ap_main->get_gb_main());
111	#endif // DEBUG
112	delete ap_main; // closes DB
113	ap_main = NULp;
114
115	exit(EXIT_SUCCESS);
116	}
117
118	static void AP_user_push_cb(AW_window *aww) {
119	ap_main->remember_user_state();
120	aww->get_root()->awar(AWAR_STACKPOINTER)->write_int(ap_main->get_user_push_counter());
121	}
122
123	static void AP_user_pop_cb(AW_window aww, TREE_canvas ntw) {
124	if (ap_main->get_user_push_counter()<=0) {
125	aw_message("No tree on stack.");
126	return;
127	}
128
129	AWT_auto_refresh allowed_on(ntw);
130	ap_main->revert_user_state();
131	ntw->request_save();
132
133	aww->get_root()->awar(AWAR_STACKPOINTER)->write_int(ap_main->get_user_push_counter());
134	if (ap_main->get_user_push_counter() <= 0) { // last tree was popped => push again
135	AP_user_push_cb(aww);
136	}
137	}
138
139	class InsertData {
140	bool abort_flag;
141	arb_progress progress;
142
143	public:
144
145	bool quick_add_flag;
146	InsertData(bool quick, long spec_count)
147	: abort_flag(false),
148	progress(GBS_global_string("Inserting %li species", spec_count), spec_count),
149	quick_add_flag(quick)
150	{}
151
152	bool aborted() const { return abort_flag; }
153	void set_aborted(bool aborted_) { abort_flag = aborted_; }
154
155	void inc() {
156	progress.inc();
157	abort_flag = progress.aborted();
158	}
159
160	arb_progress& get_progress() { return progress; }
161	};
162
163
164	static int sort_sequences_by_length(const char, long leaf0_ptr, const char, long leaf1_ptr) { // @@@ any chance to make this typesafe?
165	AP_tree_nlen leaf0 = (AP_tree_nlen)leaf0_ptr;
166	AP_tree_nlen leaf1 = (AP_tree_nlen)leaf1_ptr;
167
168	AP_FLOAT len0 = leaf0->get_seq()->weighted_base_count();
169	AP_FLOAT len1 = leaf1->get_seq()->weighted_base_count();
170
171	// longest sequence first
172	if (len0<len1) return 1;
173	if (len0>len1) return -1;
174
175	// if length equal -> determine order by species name (just to have a defined order!)
176	int cmp = strcmp(leaf1->name, leaf0->name);
177	ap_assert(cmp != 0);
178	return cmp;
179	}
180
181	static long transform_gbd_to_leaf(const char key, long val, void ) {
182	if (!val) return val;
183
184	#if defined(WARN_TODO)
185	#warning use create_linked_leaf() when impl
186	#endif
187
188	GBDATA gb_node = (GBDATA )val;
189	AP_pars_root *troot = ap_main->get_tree_root();
190	AP_tree_nlen leaf = DOWNCAST(AP_tree_nlen, troot->makeNode());
191
192	leaf->forget_origin(); // new leaf is not part of tree yet
193
194	leaf->gb_node = gb_node;
195	leaf->name = ARB_strdup(key);
196	leaf->markAsLeaf();
197
198	leaf->set_seq(troot->get_seqTemplate()->dup());
199	GB_ERROR error = leaf->get_seq()->bind_to_species(gb_node);
200	if (!error) {
201	if (leaf->get_seq()->weighted_base_count() < MIN_SEQUENCE_LENGTH) {
202	error = GBS_global_string("Species %s has too short sequence (%f, minimum is %i)",
203	key,
204	leaf->get_seq()->weighted_base_count(),
205	MIN_SEQUENCE_LENGTH);
206	}
207	}
208	if (error) {
209	GBT_message(gb_node, error);
210	destroy(leaf, troot); leaf = NULp;
211	}
212	return (long)leaf;
213	}
214
215	typedef vector<AP_tree_nlen*> InsertedSpecies;
216
217	static long toInserted(const char , long val, void cd_toInsert) {
218	InsertedSpecies toInsert = (InsertedSpecies)cd_toInsert;
219	AP_tree_nlen node = (AP_tree_nlen)val;
220
221	toInsert->push_back(node);
222	return 0;
223	}
224
225	inline int maxAllowedInsertions(int inTree) {
226	// max. species allowed to insert (in one pass) into a tree with 'inTree' leafs
227	return inTree/2;
228	}
229	inline int calcInsertNow(int toInsert, int inTree) {
230	// calculate number of species added in next pass
231	return std::min(toInsert, maxAllowedInsertions(inTree));
232	}
233
234	static long calc_steps(int toInsert, int inTree) {
235	ap_assert((toInsert+inTree) >= 2);
236
237	if (!toInsert) return 0;
238	if (!inTree) return 1 + calc_steps(toInsert-2, 2);
239
240	int edges = leafs_2_edges(inTree, UNROOTED);
241	int insertNow = calcInsertNow(toInsert, inTree);
242
243	return (edges+1)*insertNow + calc_steps(toInsert-insertNow, inTree+insertNow); // +1 for final step (=actual insertion of species)
244	}
245
246	class AP_subtree { // defines a subtree
247	AP_tree_nlen *subNode;
248	AP_tree_nlen *upNode;
249
250	bool valid() const { return subNode && upNode; }
251
252	public:
253	AP_subtree() : subNode(NULp), upNode(NULp) {}
254	AP_subtree(AP_tree_edge e, AP_tree_nlen sub_node) :
255	subNode(sub_node),
256	upNode(e->otherNode(subNode))
257	{}
258
259	AP_tree_edge *edgeToSubtree() const { ap_assert(valid()); return upNode->edgeTo(subNode); }
260	AP_tree_nlen *subtreeRoot() const { return subNode; }
261
262	void setSubtreeRoot(AP_tree_nlen *new_subtree) {
263	ap_assert(upNode->edgeTo(new_subtree));
264	subNode = new_subtree;
265	}
266	};
267
268	struct EdgeBetween : private AP_subtree {
269	// semantically same as AP_tree_edge, but survives tree-modifications which modify edges (like insert+moveNextTo/moveTo)
270
271	EdgeBetween() {}
272	EdgeBetween(AP_tree_edge *e) : AP_subtree(e, e->sonNode()) {}
273	AP_tree_edge *find() const { return edgeToSubtree(); }
274	};
275
276	struct BestEdge {
277	Mutations pars;
278	EdgeBetween between; // need to store pair of AP_tree_nlen here
279	// (using AP_tree_edge is not stable; may move elsewhere by calls insert() or moveNextTo()!)
280
281	BestEdge() : pars(-1) {}
282	BestEdge(const EdgeBetween& betw, Mutations p) : pars(p), between(betw) {}
283
284	AP_tree_edge *edge() const { return between.find(); }
285	};
286
287	struct NodeInsertOrder {
288	bool operator() (AP_tree_nlen i, AP_tree_nlen j) { return strcmp(i->name, j->name)<0; }
289	};
290
291	typedef InsertedSpecies::const_iterator InsertSpeciesIterator;
292
293	static void insert_species_into_tree(const InsertSpeciesIterator begin, const InsertSpeciesIterator end, arb_progress& progress) {
294	typedef map<AP_tree_nlen*, BestEdge> BestEdge4Node;
295	BestEdge4Node bestpos;
296
297	ap_assert(begin != end);
298
299	{
300	ap_main->remember();
301
302	EdgeChain chain(rootEdge(), ANY_EDGE, false);
303	ap_assert(chain.size()>0);
304
305	bool speciesInserted = false;
306
307	while (chain) {
308	AP_tree_edge edge = chain; ++chain;
309	edge->set_root();
310
311	EdgeBetween betweenNodes(edge);
312
313	InsertSpeciesIterator curr = begin;
314	AP_tree_nlen species = curr++;
315
316	if (speciesInserted) {
317	species->moveTo(edge);
318	}
319	else {
320	species->insert(edge->sonNode()); // edge is root-edge -> son does not matter
321	speciesInserted = true;
322	}
323
324	species->set_root(); // => only needs one combine when exchanging species
325
326	Mutations pars = rootNode()->costs();
327	BestEdge4Node::iterator found = bestpos.find(species);
328	if (found == bestpos.end() \|\| pars<found->second.pars) {
329	bestpos[species] = BestEdge(betweenNodes, pars);
330	}
331	++progress;
332
333	AP_tree_nlen *rot_node = rootNode()->get_leftson(); // rot=rest of tree
334	if (rot_node == species) {
335	rot_node = rot_node->get_brother();
336	}
337	ap_assert(rot_node->get_brother() == species);
338
339	AP_combinableSeq *rot_seq = rot_node->get_seq();
340	Mutations rot_costs = rot_node->stored_costs();
341
342	ap_assert(species->stored_costs() == 0); // leaf has no mutations
343
344	while (1) {
345	if (curr == end) break;
346
347	AP_tree_nlen nextSpec = curr++;
348	AP_combinableSeq *nextSeq = nextSpec->get_seq();
349
350	pars = nextSeq->mutations_if_combined_with(rot_seq) + rot_costs;
351	found = bestpos.find(nextSpec);
352	if (found == bestpos.end() \|\| pars<found->second.pars) {
353	bestpos[nextSpec] = BestEdge(betweenNodes, pars);
354	}
355	++progress;
356	}
357	}
358
359	ap_main->revert();
360	}
361
362	// create insert lists for each used insert position:
363	typedef list<AP_tree_nlen*> NodeList;
364	typedef map<AP_tree_edge*, NodeList> NodesAtEdge;
365
366	NodesAtEdge atEdge;
367	for (InsertSpeciesIterator s = begin; s != end; ++s) {
368	const BestEdge& best = bestpos[*s];
369	AP_tree_edge *edge = best.edge();
370
371	ap_assert(edge != NULp);
372
373	NodesAtEdge::iterator at = atEdge.find(edge);
374	if (at == atEdge.end()) {
375	atEdge[edge] = NodeList(1, *s);
376	}
377	else {
378	at->second.push_back(*s);
379	}
380	}
381
382	#if defined(DEVEL_RALF)
383	// testcode: test whether all found edges are members of the tree
384	// (got some problem with insert/REMOVE while root is next to inserted/removed node)
385
386	set<AP_tree_edge*> edgeInTree;
387	{
388	EdgeChain chain(rootEdge(), ANY_EDGE, false);
389	while (chain) {
390	AP_tree_edge edge = chain; ++chain;
391	edgeInTree.insert(edge);
392	}
393
394	for (BestEdge4Node::iterator b = bestpos.begin(); b != bestpos.end(); ++b) {
395	AP_tree_edge *e = b->second.edge();
396
397	if (edgeInTree.find(e) == edgeInTree.end()) {
398	GBK_terminate("remembered edge has been removed from tree");
399	}
400	}
401	}
402	#endif
403
404	// build list of edges where insert takes place (value=iterator into 'atEdge')
405	// => insert in determined order
406	typedef list<NodesAtEdge::iterator> InsertOrder;
407	InsertOrder insertOrder;
408	{
409	EdgeChain chain(rootEdge(), ANY_EDGE, false);
410	while (chain) {
411	AP_tree_edge edge = chain; ++chain;
412
413	NodesAtEdge::iterator at = atEdge.find(edge);
414	if (at != atEdge.end()) {
415	insertOrder.push_back(at);
416	}
417	}
418	}
419
420	typedef list<AP_subtree> OptiList;
421	OptiList optiPos;
422
423	// insert species to tree according to insert-lists:
424	for (InsertOrder::iterator o = insertOrder.begin(); o != insertOrder.end(); ++o) {
425	NodesAtEdge::iterator e = *o;
426	AP_tree_edge *edge = e->first;
427	NodeList& nodes = e->second;
428
429	edge->set_root();
430
431	AP_tree_nlen *brother = edge->sonNode();
432	size_t nodes_size = nodes.size();
433
434	#if defined(ASSERTION_USED)
435	ap_assert(bestpos[nodes.front()].edge() == edge);
436	#endif
437
438	if (nodes_size == 1) {
439	nodes.front()->insert(brother);
440	ASSERT_VALID_TREE(rootNode());
441	}
442	else {
443	bool atLeaf = brother->is_leaf();
444	if (!atLeaf && edge->is_leaf_edge()) { // at leaf edge -> make sure brother points to leaf node
445	brother = edge->notSonNode();
446	ap_assert(brother->is_leaf());
447	atLeaf = true;
448	}
449
450	#if defined(UNIT_TESTS)
451	if (RUNNING_TEST()) {
452	// use a determined order to insert multiple species at one position.
453	// Does not produce "better" topologies, just makes result independent from insert order.
454	typedef vector<AP_tree_nlen*> NodeVector;
455
456	NodeVector toSort(nodes.begin(), nodes.end());
457	sort(toSort.begin(), toSort.end(), NodeInsertOrder());
458	nodes = NodeList(toSort.begin(), toSort.end());
459	}
460	#endif
461
462	AP_tree_nlen *at = brother;
463	for (NodeList::iterator n = nodes.begin(); n != nodes.end(); ++n) {
464	(*n)->insert(at);
465	at = *n; // only insert 1st node at 'brother', insert following nodes next to previously added nodes
466	}
467
468	ASSERT_VALID_TREE(rootNode());
469
470	AP_tree_nlen *ourFather = brother->get_father();
471	AP_tree_nlen *addedSubtree = brother->get_brother(); // contains all added species
472	ap_assert(addedSubtree->is_ancestor_of(at));
473
474	if (atLeaf) {
475	// if inserted at leaf edge -> perform NNI at parent edge (i.e. including the leaf)
476	AP_tree_edge *toRest = ourFather->nextEdge();
477	for (int i = 0; i<2; ++i) {
478	AP_tree_nlen *rest = toRest->otherNode(ourFather);
479	if (rest != brother && rest != addedSubtree) {
480	break;
481	}
482	toRest = ourFather->nextEdge(toRest);
483	}
484
485	optiPos.push_back(AP_subtree(toRest, ourFather));
486	ap_assert(optiPos.back().subtreeRoot() == ourFather);
487	}
488	else {
489	if (nodes_size>2) { // if inserted at internal edge && only 2 species inserted -> NNI makes no sense
490	// Store (directed) edge to brother (for later optimization of subtree):
491	AP_tree_edge *subEdge = ourFather->edgeTo(addedSubtree);
492	optiPos.push_back(AP_subtree(subEdge, addedSubtree));
493	ap_assert(optiPos.back().subtreeRoot() == addedSubtree);
494	}
495	}
496	}
497	progress.inc_by(nodes_size);
498	}
499
500	// Optimize all inserts of multiple species at one position:
501	{
502	arb_suppress_progress suppress_child; // suppress implicit progress count caused by nni_rec
503
504	AP_FLOAT curr_pars = rootNode()->costs();
505	AP_FLOAT prev_pars = curr_pars;
506
507	int loop = 0;
508
509	do {
510	++loop;
511	prev_pars = curr_pars;
512	for (OptiList::iterator op = optiPos.begin(); op != optiPos.end(); ++op) {
513	AP_tree_edge *subtreeEdge = op->edgeToSubtree();
514	AP_tree_nlen *subtreeRoot = op->subtreeRoot();
515
516	subtreeEdge->set_root();
517	ap_assert(subtreeEdge->isConnectedTo(subtreeRoot));
518	AP_tree_nlen *father = subtreeEdge->otherNode(subtreeRoot);
519
520	AP_FLOAT this_pars;
521	while (1) {
522	ap_assert(subtreeEdge->isConnectedTo(father)); // otherwise block fails
523	this_pars = subtreeEdge->nni_rec(SKIP_LEAF_EDGES, AP_BL_NNI_ONLY, father, false);
524	if (!(this_pars<curr_pars)) {
525	ap_assert(!(this_pars>curr_pars));
526	break;
527	}
528	curr_pars = this_pars;
529	}
530
531	ap_assert(subtreeEdge->isConnectedTo(father)); // otherwise next command fails
532	AP_tree_nlen *newSubtreeRoot = subtreeEdge->otherNode(father);
533	if (newSubtreeRoot != subtreeRoot) {
534	op->setSubtreeRoot(newSubtreeRoot);
535	}
536	}
537	}
538	while (curr_pars<prev_pars);
539	}
540	}
541
542	static void insert_all_species_into_tree(GB_HASH*& hash) {
543	// inserts all species (from hash) into tree
544
545	AP_tree_nlen *tree = rootNode();
546
547	int inTree = tree ? tree->count_leafs() : 0;
548	int toInsert = GBS_hash_elements(hash);
549
550	ap_assert(toInsert);
551
552	long steps = calc_steps(toInsert, inTree);
553	arb_progress progress(steps);
554
555	// move species to insert to a stack
556	InsertedSpecies speciesToInsert;
557	speciesToInsert.reserve(toInsert);
558
559	if (maxAllowedInsertions(inTree)<toInsert) {
560	// insert longest sequences first
561	GBS_hash_do_sorted_loop(hash, toInserted, sort_sequences_by_length, &speciesToInsert);
562	}
563	else {
564	// insert all sequences (order should not matter)
565	GBS_hash_do_loop(hash, toInserted, &speciesToInsert);
566	}
567	GBS_free_hash(hash);
568	hash = NULp;
569
570	ap_assert(toInsert != 2); // @@@ need to test this case
571
572	InsertSpeciesIterator curr = speciesToInsert.begin();
573	InsertSpeciesIterator end = speciesToInsert.end();
574
575	AP_tree_edge *oldRootEdge = NULp;
576	if (!tree) { // create initial tree
577	AP_pars_root *troot = ap_main->get_tree_root();
578
579	AP_tree_nlen s1 = curr++;
580	AP_tree_nlen s2 = curr++;
581
582	s1->initial_insert(s2, troot);
583
584	inTree = 2;
585	toInsert -= 2;
586
587	++progress;
588	}
589	else {
590	oldRootEdge = rootEdge();
591	}
592
593	ASSERT_VALID_TREE(rootNode());
594
595	while (1) {
596	int insertNow = calcInsertNow(toInsert, inTree);
597	ap_assert(insertNow<=toInsert);
598	if (insertNow == toInsert) break;
599
600	{
601	InsertSpeciesIterator partEnd = curr;
602	advance(partEnd, insertNow);
603
604	insert_species_into_tree(curr, partEnd, progress);
605	curr = partEnd;
606	}
607
608	toInsert -= insertNow;
609	inTree += insertNow;
610	}
611
612	insert_species_into_tree(curr, end, progress);
613
614	if (oldRootEdge) oldRootEdge->set_root(); // set root back to old position
615	}
616
617	enum AddWhat {
618	NT_ADD_MARKED,
619	NT_ADD_SELECTED,
620	};
621
622	static void nt_add(AWT_graphic_parsimony *agt, AddWhat what, bool quick) {
623	GB_ERROR error = NULp;
624
625	AP_tree *oldrootleft = NULp;
626	AP_tree *oldrootright = NULp;
627	{
628	AP_tree_nlen *root = rootNode();
629	if (root) {
630	root->reset_subtree_layout();
631	oldrootleft = root->get_leftson();
632	oldrootright = root->get_rightson();
633	}
634	}
635
636	GB_HASH *hash = NULp;
637	GBDATA *gb_main = agt->get_gbmain();
638	{
639	GB_transaction ta(gb_main);
640	switch (what) {
641	case NT_ADD_SELECTED: {
642	char *name = GBT_readOrCreate_string(gb_main, AWAR_SPECIES_NAME, "");
643	if (name && strlen(name)) {
644	GBDATA *gb_species = GBT_find_species(gb_main, name);
645	if (gb_species) {
646	hash = GBS_create_hash(1, GB_MIND_CASE);
647	GBS_write_hash(hash, name, (long)gb_species);
648	}
649	else error = GBS_global_string("Selected Species (%s) not found", name);
650	}
651	else error = "Please select a species";
652	free(name);
653	break;
654	}
655	case NT_ADD_MARKED: {
656	hash = GBT_create_marked_species_hash(gb_main);
657	break;
658	}
659	}
660	}
661
662	if (!error) {
663	ap_assert(hash);
664
665	arb_progress progress(quick ? "Quick add" : "Add + NNI");
666
667	NT_remove_species_in_tree_from_hash(rootNode(), hash);
668
669	size_t species_count = GBS_hash_elements(hash);
670	InsertPerfMeter insertPerf("(quick-)add", species_count);
671
672	{
673	GB_transaction ta(gb_main);
674	GBS_hash_do_loop(hash, transform_gbd_to_leaf, NULp);
675	}
676	{
677	size_t skipped = species_count - GBS_hash_elements(hash);
678	if (skipped) {
679	GBT_message(gb_main, GBS_global_string("Skipped %zu species (no data?)", skipped));
680	}
681	}
682	if (GBS_hash_elements(hash)) {
683	insert_all_species_into_tree(hash);
684	}
685	else {
686	GBT_message(gb_main, "No species (left) to insert");
687	}
688
689	if (rootNode()) {
690	if (oldrootleft) {
691	if (oldrootleft->father == oldrootright) oldrootleft->set_root();
692	else oldrootright->set_root();
693	}
694	else {
695	ARB_edge innermost = rootNode()->get_tree_root()->find_innermost_edge();
696	innermost.set_root();
697	}
698
699	if (!quick) {
700	arb_suppress_progress quiet; // @@@ use weighted progress (#789) for loop below
701
702	Mutations pars_prev = rootNode()->costs();
703	rootNode()->compute_tree(); // see AP_tree_edge.cxx@flags_broken_by_moveNextTo
704	progress.subtitle("local optimize (repeated NNI)");
705	while (1) {
706	rootEdge()->nni_rec(EdgeSpec(SKIP_UNMARKED_EDGES\|SKIP_LEAF_EDGES), AP_BL_NNI_ONLY, NULp, true);
707	Mutations pars_curr = rootNode()->costs();
708	if (pars_curr == pars_prev) break;
709	ap_assert(pars_curr<pars_prev);
710	pars_prev = pars_curr;
711	}
712	}
713
714	{
715	arb_suppress_progress ignore;
716	rootEdge()->calc_branchlengths();
717	}
718
719	ASSERT_VALID_TREE(rootNode());
720	rootNode()->compute_tree();
721	}
722	else {
723	error = "Tree lost (no leafs left)";
724	}
725
726	insertPerf.dump(stdout);
727	}
728
729	if (hash) GBS_free_hash(hash);
730	if (error) aw_message(error);
731
732	// @@@ quick-add w/o NNI should sort according to original tree
733	agt->reorderTree(BIG_BRANCHES_TO_TOP);
734	}
735
736	// ------------------------------------------
737	// Adding partial sequences to tree
738
739	class PartialSequence {
740	GBDATA *gb_species;
741	mutable AP_tree_nlen *self; // self converted to leaf (ready for insertion)
742	const AP_tree_nlen *best_full_match; // full sequence position which matched best
743	long overlap; // size of overlapping region
744	long penalty; // weighted mismatches
745	bool released;
746	bool multi_match;
747	string multi_list; // list of equal-rated insertion-points (not containing self)
748
749	AP_tree_nlen *get_self() const {
750	if (!self) {
751	ap_assert(!released); // request not possible, because leaf has already been released!
752
753	self = (AP_tree_nlen*)transform_gbd_to_leaf(GBT_read_name(gb_species), (long)gb_species, NULp);
754	ap_assert(self);
755	}
756	return self;
757	}
758
759	public:
760	PartialSequence(GBDATA *gb_species_) :
761	gb_species(gb_species_),
762	self(NULp),
763	best_full_match(NULp),
764	overlap(0),
765	penalty(LONG_MAX),
766	released(false),
767	multi_match(false)
768	{}
769	PartialSequence(const PartialSequence& other)
770	: gb_species(other.gb_species),
771	self(other.self),
772	best_full_match(other.best_full_match),
773	overlap(other.overlap),
774	penalty(other.penalty),
775	released(other.released),
776	multi_match(other.multi_match),
777	multi_list(other.multi_list)
778	{
779	ap_assert(!self); // copying self not implemented
780	}
781	DECLARE_ASSIGNMENT_OPERATOR(PartialSequence);
782	~PartialSequence() { ap_assert(!self); }
783
784	GBDATA *get_species() const { return gb_species; }
785	const AP_tree_nlen *get_best_match() const { return best_full_match; }
786	AP_FLOAT get_branchlength() const { return AP_FLOAT(penalty)/overlap; }
787	void test_match(const AP_tree_nlen *leaf_full);
788	bool is_multi_match() const { return multi_match; }
789
790	const char *get_name() const {
791	const char *name = get_self()->name;
792	ap_assert(name);
793	return name;
794	}
795
796	string get_multilist() const {
797	ap_assert(is_multi_match());
798	return string(best_full_match->name)+multi_list;
799	}
800
801	AP_tree_nlen *release() {
802	AP_tree_nlen *s = self;
803	self = NULp;
804	released = true;
805	return s;
806	}
807
808	void dump(const char *whichMatch) const {
809	ap_assert(best_full_match);
810	printf("%s match for '%s' is '%s' (overlap=%li penalty=%li)\n",
811	whichMatch, get_name(), best_full_match->name,
812	overlap, penalty);
813	}
814
815	};
816
817	void PartialSequence::test_match(const AP_tree_nlen *leaf_full) {
818	long curr_overlap;
819	long curr_penalty;
820
821	leaf_full->get_seq()->partial_match(get_self()->get_seq(), &curr_overlap, &curr_penalty);
822
823	bool better = false;
824
825	if (curr_overlap > overlap) {
826	better = true;
827	}
828	else if (curr_overlap == overlap) {
829	if (curr_penalty<penalty) {
830	better = true;
831	}
832	else if (curr_penalty == penalty) {
833	// found two equal-rated insertion points -> store data for warning
834	#if defined(DEBUG)
835	if (!multi_match) dump("better");
836	printf("Another equal match is against '%s' (overlap=%li penalty=%li)\n", leaf_full->name, curr_overlap, curr_penalty);
837	#endif // DEBUG
838
839	multi_match = true;
840	multi_list.append(1, '/');
841	multi_list.append(leaf_full->name);
842	}
843	}
844
845	if (better) {
846	overlap = curr_overlap;
847	penalty = curr_penalty;
848	best_full_match = leaf_full;
849	multi_match = false;
850	multi_list = "";
851
852	#if defined(DEBUG)
853	dump("better");
854	#endif
855	}
856	#if defined(DEBUG)
857	else if (!multi_match) {
858	printf("Worse match against '%s' (overlap=%li penalty=%li)\n", leaf_full->name, curr_overlap, curr_penalty);
859	}
860	#endif
861	}
862
863	static GB_ERROR nt_best_partial_match_rec(list<PartialSequence>& partial, const AP_tree_nlen *tree) {
864	GB_ERROR error = NULp;
865
866	if (tree) {
867	if (tree->is_leaf() && tree->name) {
868	if (tree->gb_node) {
869	int is_partial = GBT_is_partial(tree->gb_node, 0, true); // marks undef as 'full sequence'
870	if (is_partial == 0) { // do not consider other partial sequences
871	list<PartialSequence>::iterator i = partial.begin();
872	list<PartialSequence>::iterator e = partial.end();
873	for (; i != e; ++i) {
874	i->test_match(tree);
875	}
876	}
877	else if (is_partial == -1) {
878	error = GB_await_error();
879	}
880	}
881	}
882	else {
883	error = nt_best_partial_match_rec(partial, tree->get_leftson());
884	if (!error) error = nt_best_partial_match_rec(partial, tree->get_rightson());
885	}
886	}
887	return error;
888	}
889
890	static void count_partial_and_full(const AP_tree_nlen at, int partial, int full, int zombies, int default_value, bool define_if_undef) {
891	if (at->is_leaf()) {
892	if (at->gb_node) {
893	int is_partial = GBT_is_partial(at->gb_node, default_value, define_if_undef);
894	if (is_partial) ++(*partial);
895	else ++(*full);
896	}
897	else {
898	++(*zombies);
899	}
900	}
901	else {
902	count_partial_and_full(at->get_leftson(), partial, full, zombies, default_value, define_if_undef);
903	count_partial_and_full(at->get_rightson(), partial, full, zombies, default_value, define_if_undef);
904	}
905	}
906
907	static const AP_tree_nlen find_least_deep_leaf(const AP_tree_nlen at, int depth, int *min_depth) {
908	if (depth >= *min_depth) {
909	return NULp; // already found better or equal
910	}
911
912	if (at->is_leaf()) {
913	if (at->gb_node) {
914	*min_depth = depth;
915	return at;
916	}
917	return NULp;
918	}
919
920	const AP_tree_nlen *left = find_least_deep_leaf(at->get_leftson(), depth+1, min_depth);
921	const AP_tree_nlen *right = find_least_deep_leaf(at->get_rightson(), depth+1, min_depth);
922
923	return right ? right : left;
924	}
925	inline AP_tree_nlen find_least_deep_leaf(AP_tree_nlen at, int depth, int *min_depth) {
926	return const_cast<AP_tree_nlen>(find_least_deep_leaf(const_cast<const AP_tree_nlen>(at), depth, min_depth));
927	}
928
929	static void push_partial(const char , long val, void cd_partial) {
930	list<PartialSequence> partial = reinterpret_cast<list<PartialSequence> >(cd_partial);
931	partial->push_back(PartialSequence((GBDATA*)val));
932	}
933
934	// -------------------------------
935	// Add Partial sequences
936
937	static void nt_add_partial(AWT_graphic_parsimony *agt) {
938	GB_ERROR error = NULp;
939	GBDATA *gb_main = agt->get_gbmain();
940
941	GB_begin_transaction(gb_main);
942
943	int full_marked_sequences = 0;
944
945	arb_progress part_add_progress("Adding partial sequences");
946
947	{
948	list<PartialSequence> partial;
949	{
950	GB_HASH *partial_hash = GBS_create_hash(GBT_get_species_count(gb_main), GB_MIND_CASE);
951
952	int marked_found = 0;
953	int partial_marked_sequences = 0;
954	int no_data = 0; // no data in alignment
955
956	for (GBDATA *gb_marked = GBT_first_marked_species(gb_main);
957	!error && gb_marked;
958	gb_marked = GBT_next_marked_species(gb_marked))
959	{
960	++marked_found;
961
962	if (GBT_find_sequence(gb_marked, ap_main->get_aliname())) { // species has sequence in alignment
963	const char *name = GBT_read_name(gb_marked);
964
965	switch (GBT_is_partial(gb_marked, 1, true)) { // marks undef as 'partial sequence'
966	case 0: { // full sequences
967	GBT_message(gb_main, GBS_global_string("'%s' is a full sequence (cannot add partial)", name));
968	++full_marked_sequences;
969	break;
970	}
971	case 1: // partial sequences
972	++partial_marked_sequences;
973	GBS_write_hash(partial_hash, name, (long)gb_marked);
974	break;
975	case -1: // error
976	error = GB_await_error();
977	break;
978	default:
979	ap_assert(0);
980	break;
981	}
982	}
983	else {
984	no_data++;
985	}
986	}
987
988	if (!error && !marked_found) error = "There are no marked species";
989
990	if (!error) {
991	NT_remove_species_in_tree_from_hash(rootNode(), partial_hash); // skip all species which are in tree
992	GBS_hash_do_const_loop(partial_hash, push_partial, &partial); // build partial list from hash
993
994	int partials_already_in_tree = partial_marked_sequences - partial.size();
995
996	if (no_data>0) GBT_message(gb_main, GBS_global_string("%i marked species have no data in '%s'", no_data, ap_main->get_aliname()));
997	if (full_marked_sequences>0) GBT_message(gb_main, GBS_global_string("%i marked species are declared full sequences", full_marked_sequences));
998	if (partials_already_in_tree>0) GBT_message(gb_main, GBS_global_string("%i marked species are already in tree", partials_already_in_tree));
999
1000	if (partial.empty()) error = "No species left to add";
1001	}
1002
1003	GBS_free_hash(partial_hash);
1004	}
1005
1006	if (!error) error = GBT_add_new_changekey(gb_main, "ARB_partial", GB_INT);
1007
1008	if (!error) {
1009	rootNode()->reset_subtree_layout();
1010
1011	// find best matching full sequence for each partial sequence
1012	error = nt_best_partial_match_rec(partial, rootNode());
1013
1014	list<PartialSequence>::iterator i = partial.begin();
1015	list<PartialSequence>::iterator e = partial.end();
1016
1017	arb_progress part_insert_progress(partial.size());
1018
1019	#if defined(DEBUG)
1020	// show results :
1021	for (; i != e; ++i) i->dump("best");
1022	i = partial.begin();
1023	#endif // DEBUG
1024
1025	for (; i != e && !error; ++i) {
1026	const char *name = i->get_name();
1027
1028	if (i->is_multi_match()) {
1029	GBT_message(gb_main, GBS_global_string("Insertion of '%s' is ambiguous.\n"
1030	"(took first of equal scored insertion points: %s)",
1031	name, i->get_multilist().c_str()));
1032	}
1033
1034	AP_tree_nlen *part_leaf = i->release();
1035	AP_tree_nlen full_seq = const_cast<AP_tree_nlen>(i->get_best_match());
1036	AP_tree_nlen *brother = full_seq->get_brother();
1037	int is_partial = 0;
1038	AP_tree_nlen *target = NULp;
1039
1040	if (brother->is_leaf()) {
1041	if (brother->gb_node) {
1042	is_partial = GBT_is_partial(brother->gb_node, 0, true);
1043
1044	if (is_partial) { // brother is partial sequence
1045	target = brother; // insert as brother of brother
1046	}
1047	else {
1048	target = full_seq; // insert as brother of full_seq
1049	}
1050	}
1051	else {
1052	error = "There are zombies in your tree - please remove them";
1053	}
1054	}
1055	else {
1056	int partial_count = 0;
1057	int full_count = 0;
1058	int zombie_count = 0;
1059
1060	count_partial_and_full(brother, &partial_count, &full_count, &zombie_count, 0, true);
1061
1062	if (zombie_count) {
1063	error = "There are zombies in your tree - please remove them";
1064	}
1065	else if (full_count) {
1066	// brother is a subtree containing full sequences
1067	// -> add new brother to full_seq found above
1068	target = full_seq;
1069	}
1070	else { // brother subtree only contains partial sequences
1071	// find one of the least-deep leafs
1072	int depth = INT_MAX;
1073	target = find_least_deep_leaf(brother, 0, &depth);
1074	is_partial = 1;
1075	}
1076	}
1077
1078
1079	if (!error) {
1080	#if defined(DEBUG)
1081	printf("inserting '%s'\n", name);
1082	#endif // DEBUG
1083	part_leaf->insert(target);
1084
1085	// we need to create the sequence of the father node!
1086	AP_tree_nlen *father = part_leaf->get_father();
1087	father->costs();
1088
1089	// ensure full-sequence is always on top
1090	if (father->rightson == target) {
1091	father->swap_sons();
1092	}
1093
1094	if (!error) { // now correct the branch lengths modified by insert()
1095	// calc the original branchlen (of target leaf branch)
1096	GBT_LEN orglen = father->get_branchlength()+target->get_branchlength();
1097
1098	if (is_partial) { // we have a subtree of partial sequences
1099	target->set_branchlength(orglen); // restore original branchlength
1100	father->set_branchlength(0); // all father branches are zero length
1101	}
1102	else { // we have a subtree of one full+one partial sequence
1103	ap_assert(full_seq->get_father() == father);
1104
1105	father->set_branchlength(orglen); // father branch represents original length (w/o partial seq)
1106	full_seq->set_branchlength(0); // full seq has no sub-branch length
1107	}
1108	part_leaf->set_branchlength(i->get_branchlength());
1109	printf("Adding with branchlength=%f\n", i->get_branchlength());
1110	}
1111	}
1112	else {
1113	destroy(part_leaf);
1114	}
1115
1116	part_insert_progress.inc_and_check_user_abort(error);
1117	}
1118	}
1119	}
1120
1121	if (full_marked_sequences) {
1122	GBT_message(gb_main, GBS_global_string("%i marked full sequences were not added", full_marked_sequences));
1123	}
1124
1125	if (error) {
1126	GBT_message(gb_main, error);
1127	GB_abort_transaction(gb_main);
1128	}
1129	else {
1130	GB_commit_transaction(gb_main);
1131	agt->exports.request_save();
1132	}
1133	}
1134
1135	static void NT_add_partial_and_update(UNFIXED, TREE_canvas *ntw) {
1136	AWT_auto_refresh allowed_on(ntw);
1137	nt_add_partial(AWT_TREE_PARS(ntw));
1138	}
1139
1140	// -------------------------------
1141	// add marked / selected
1142
1143	static void nt_add_and_update(AWT_canvas *ntw, AddWhat what, bool quick) {
1144	AWT_auto_refresh allowed_on(ntw);
1145	nt_add(AWT_TREE_PARS(ntw), what, quick);
1146	}
1147
1148	static void NT_add_and_NNI(UNFIXED, TREE_canvas *ntw, AddWhat what) { nt_add_and_update(ntw, what, false); }
1149	static void NT_add_quick (UNFIXED, TREE_canvas *ntw, AddWhat what) { nt_add_and_update(ntw, what, true); }
1150
1151	// ------------------------------------------
1152	// remove and add marked / selected
1153
1154	static void nt_reAdd(AWT_graphic_parsimony *agt, AddWhat what, bool quick) {
1155	if (agt->get_root_node()) {
1156	ap_assert(what == NT_ADD_MARKED); // code below will misbehave for NT_ADD_SELECTED
1157	agt->get_tree_root()->remove_leafs(AWT_REMOVE_MARKED);
1158	nt_add(agt, what, quick);
1159	}
1160	}
1161
1162	static void nt_reAdd_and_update(AWT_canvas *ntw, AddWhat what, bool quick) {
1163	AWT_auto_refresh allowed_on(ntw);
1164	nt_reAdd(AWT_TREE_PARS(ntw), what, quick);
1165	}
1166
1167	static void NT_reAdd_and_NNI(UNFIXED, TREE_canvas *ntw, AddWhat what) { nt_reAdd_and_update(ntw, what, false); }
1168	static void NT_reAdd_quick (UNFIXED, TREE_canvas *ntw, AddWhat what) { nt_reAdd_and_update(ntw, what, true); }
1169
1170	// --------------------------------------------------------------------------------
1171
1172	static void calc_branchlengths_and_reorder(AWT_graphic_parsimony *agt) {
1173	arb_progress progress("Calculating branchlengths");
1174	rootEdge()->calc_branchlengths();
1175	agt->reorderTree(BIG_BRANCHES_TO_TOP);
1176	}
1177
1178	static void NT_calc_branchlengths_reorder_and_update(AW_window , TREE_canvas ntw) {
1179	AWT_auto_refresh allowed_on(ntw);
1180	calc_branchlengths_and_reorder(AWT_TREE_PARS(ntw));
1181	}
1182
1183	static void NT_bootstrap(AW_window , TREE_canvas ntw, bool limit_only) {
1184	arb_progress progress("Calculating bootstrap limit");
1185	AWT_auto_refresh allowed_on(ntw);
1186	AP_BL_MODE mode = AP_BL_MODE((limit_only ? AP_BL_BOOTSTRAP_LIMIT : AP_BL_BOOTSTRAP_ESTIMATE)\|AP_BL_BL_ONLY);
1187
1188	rootEdge()->nni_rec(ANY_EDGE, mode, NULp, true);
1189	AWT_graphic_tree *agt = AWT_TREE(ntw);
1190	agt->reorderTree(BIG_BRANCHES_TO_TOP);
1191	agt->set_logical_root_to(agt->get_root_node());
1192	}
1193
1194	static void optimizeTree(AWT_graphic_parsimony *agt, const KL_Settings& settings) {
1195	arb_progress progress("Optimizing tree");
1196	agt->get_parsimony().optimize_tree(rootNode(), settings, progress);
1197	ASSERT_VALID_TREE(rootNode());
1198	calc_branchlengths_and_reorder(agt);
1199	}
1200	static void NT_optimize(AW_window , TREE_canvas ntw) {
1201	AWT_auto_refresh allowed_on(ntw);
1202	optimizeTree(AWT_TREE_PARS(ntw), KL_Settings(ntw->awr));
1203	}
1204
1205	static void recursiveNNI(AWT_graphic_parsimony *agt, EdgeSpec whichEdges) {
1206	arb_progress progress("Recursive NNI");
1207	Mutations orgPars = rootNode()->costs();
1208	Mutations prevPars = orgPars;
1209	progress.subtitle(GBS_global_string("best=%li", orgPars));
1210
1211	{
1212	arb_suppress_progress quiet; // @@@ use weighted progress (#789) for loop below
1213
1214	while (!progress.aborted()) {
1215	Mutations currPars = rootEdge()->nni_rec(whichEdges, AP_BL_NNI_ONLY, NULp, true);
1216	if (currPars == prevPars) break; // no improvement -> abort
1217	progress.subtitle(GBS_global_string("best=%li (gain=%li)", currPars, orgPars-currPars));
1218	prevPars = currPars;
1219	}
1220	calc_branchlengths_and_reorder(agt);
1221	}
1222	}
1223
1224	static void NT_recursiveNNI(AW_window , TREE_canvas ntw) {
1225	AWT_auto_refresh allowed_on(ntw);
1226	EdgeSpec whichEdges = KL_Settings(ntw->awr).whichEdges;
1227	recursiveNNI(AWT_TREE_PARS(ntw), whichEdges);
1228	}
1229
1230	static int calculate_default_random_repeat(long leafs) {
1231	double balanced_depth = log10(leafs) / log10(2);
1232	int repeat = int(balanced_depth*2.0 + .5);
1233	if (repeat<1) repeat = 1;
1234	return repeat;
1235	}
1236
1237	static void update_random_repeat(AW_root awr, AWT_graphic_parsimony agt) {
1238	long leafs = agt->get_root_node()->count_leafs();
1239	int repeat = calculate_default_random_repeat(leafs);
1240	awr->awar(AWAR_RAND_REPEAT)->write_int(repeat);
1241	}
1242
1243	static void mixtree_and_calclengths(AWT_graphic_parsimony *agt, int repeat, int percent, EdgeSpec whichEdges) {
1244	arb_progress progress("Randomizing tree", 2L);
1245
1246	// @@@ possible candidate for weighted progress (#789)
1247	progress.subtitle("mixing");
1248	rootEdge()->mixTree(repeat, percent, whichEdges);
1249	++progress;
1250
1251	progress.subtitle("calculating branchlengths");
1252	rootEdge()->calc_branchlengths();
1253	++progress;
1254
1255	agt->exports.request_save();
1256	}
1257
1258	static void randomMixTree(AW_window aww, TREE_canvas ntw) {
1259	AWT_auto_refresh allowed_on(ntw);
1260	AW_root *awr = aww->get_root();
1261
1262	mixtree_and_calclengths(AWT_TREE_PARS(ntw), awr->awar(AWAR_RAND_REPEAT)->read_int(), awr->awar(AWAR_RAND_PERCENT)->read_int(), KL_Settings(awr).whichEdges);
1263	{
1264	ARB_edge newRootEdge = rootNode()->get_tree_root()->find_innermost_edge();
1265	newRootEdge.son()->set_root();
1266	}
1267	AWT_TREE_PARS(ntw)->reorderTree(BIG_BRANCHES_TO_TOP);
1268	}
1269
1270
1271	static AWT_config_mapping_def optimizer_config_mapping[] = {
1272	{ AWAR_OPTI_MARKED_ONLY, "marked_only" },
1273	{ AWAR_OPTI_SKIP_FOLDED, "skip_folded" },
1274
1275	// { AWAR_RAND_REPEAT, "rand_repeat" }, // do not store (use treesize-dependent default)
1276	{ AWAR_RAND_PERCENT, "rand_percent" },
1277
1278	{ AWAR_KL_MAXDEPTH, "maxdepth" },
1279	{ AWAR_KL_INCDEPTH, "incdepth" },
1280
1281	{ AWAR_KL_STATIC_ENABLED, "static" },
1282	{ AWAR_KL_STATIC_DEPTH1, "s_depth1" },
1283	{ AWAR_KL_STATIC_DEPTH2, "s_depth2" },
1284	{ AWAR_KL_STATIC_DEPTH3, "s_depth3" },
1285	{ AWAR_KL_STATIC_DEPTH4, "s_depth4" },
1286	{ AWAR_KL_STATIC_DEPTH5, "s_depth5" },
1287
1288	{ AWAR_KL_DYNAMIC_ENABLED, "dynamic" },
1289	{ AWAR_KL_DYNAMIC_START, "start" },
1290	{ AWAR_KL_DYNAMIC_MAXX, "maxx" },
1291	{ AWAR_KL_DYNAMIC_MAXY, "maxy" },
1292
1293	{ NULp, NULp }
1294	};
1295
1296	static AWT_predefined_config optimizer_predefined_configs[] = {
1297	{
1298	"*minimum_static_reduction",
1299	"Sets paths allowed by static reduction to maximum\n(causing the minimal reduction)",
1300	"s_depth1='8';s_depth2='6';s_depth3='6';s_depth4='6';s_depth5='6';static='1'" // only defines/affects settings related to static path reduction
1301	},
1302	{
1303	"*whole_tree_level8",
1304	"Level-8-optimization of whole tree\n(no path reduction)",
1305	"dynamic='0';incdepth='0';marked_only='0';maxdepth='8';skip_folded='0';static='0'"
1306	},
1307	{ NULp, NULp, NULp }
1308	};
1309
1310	static AW_window createOptimizeWindow(AW_root aw_root, TREE_canvas *ntw) {
1311	AW_window_simple *aws = new AW_window_simple;
1312	aws->init(aw_root, "TREE_OPTIMIZE", "Tree optimization");
1313	aws->load_xfig("pars/tree_opti.fig");
1314
1315	aws->at("close");
1316	aws->callback(AW_POPDOWN);
1317	aws->create_button("CLOSE", "CLOSE", "C");
1318
1319	aws->at("help");
1320	aws->callback(makeHelpCallback("pa_optimizer.hlp"));
1321	aws->create_button("HELP", "HELP", "H");
1322
1323	aws->at("marked");
1324	aws->label("Only subtrees containing marked species");
1325	aws->create_toggle(AWAR_OPTI_MARKED_ONLY);
1326
1327	aws->at("folded");
1328	aws->label("Do not modify folded subtrees");
1329	aws->create_toggle(AWAR_OPTI_SKIP_FOLDED);
1330
1331	aws->button_length(18);
1332
1333	aws->at("rec_nni");
1334	aws->callback(makeWindowCallback(NT_recursiveNNI, ntw));
1335	aws->create_button("REC_NNI", "Recursive NNI", "N");
1336
1337	aws->at("heuristic");
1338	aws->callback(makeWindowCallback(NT_optimize, ntw));
1339	aws->create_button("HEURISTIC", "Heuristic\noptimizer", "H");
1340
1341	aws->at("config");
1342	AWT_insert_config_manager(aws, AW_ROOT_DEFAULT, "treeopti", optimizer_config_mapping, NULp, optimizer_predefined_configs);
1343
1344	aws->at("settings");
1345	aws->callback(makeCreateWindowCallback(create_kernighan_properties_window));
1346	aws->create_button("SETTINGS", "Settings", "S");
1347
1348	aws->at("randomize");
1349	aws->callback(makeWindowCallback(randomMixTree, ntw));
1350	aws->create_button("RANDOMIZE", "Randomize tree", "R");
1351
1352	aws->button_length(5);
1353
1354	aws->at("repeat"); aws->create_input_field(AWAR_RAND_REPEAT);
1355	aws->at("percent"); aws->create_input_field(AWAR_RAND_PERCENT);
1356
1357	return aws;
1358	}
1359
1360	// -----------------------
1361	// test functions
1362
1363	#if defined(TESTMENU)
1364	static void refreshTree(AWT_canvas *ntw) {
1365	GB_transaction ta(ntw->gb_main);
1366	AWT_auto_refresh allowed_on(ntw);
1367	ntw->request_save_and_zoom_reset();
1368	}
1369
1370	static void setBranchlens(AP_tree_nlen *node, double newLen) {
1371	node->setBranchlen(newLen, newLen);
1372
1373	if (!node->is_leaf()) {
1374	setBranchlens(node->get_leftson(), newLen);
1375	setBranchlens(node->get_rightson(), newLen);
1376	}
1377	}
1378
1379	static void TESTMENU_setBranchlen(AW_window , AWT_canvas ntw) {
1380	AP_tree_nlen *root = rootNode();
1381
1382	setBranchlens(root, 1.0);
1383	refreshTree(ntw);
1384	}
1385
1386	static void TESTMENU_treeStats(AW_window *) {
1387	ARB_tree_info tinfo;
1388	AP_tree_nlen *root = rootNode();
1389
1390	if (root) {
1391	{
1392	GB_transaction ta(root->get_tree_root()->get_gb_main());
1393	root->calcTreeInfo(tinfo);
1394	}
1395
1396	puts("Tree stats:");
1397
1398	printf("nodes =%6zu\n", tinfo.nodes());
1399	printf(" inner =%6zu\n", tinfo.innerNodes);
1400	printf(" groups =%6zu\n", tinfo.groups);
1401	printf(" leafs =%6zu\n", tinfo.leafs);
1402	printf(" unlinked =%6zu (zombies?)\n", tinfo.unlinked);
1403	printf(" linked =%6zu\n", tinfo.linked());
1404	printf(" marked =%6zu\n", tinfo.marked);
1405	}
1406	else {
1407	puts("No tree");
1408	}
1409	}
1410
1411	static void TESTMENU_sortTreeByName(AW_window , AWT_canvas ntw) {
1412	AP_tree_nlen *root = rootNode();
1413
1414	root->sortByName();
1415	refreshTree(ntw);
1416	}
1417
1418	static void init_TEST_menu(AW_window_menu_modes awm, AWT_canvas ntw) {
1419	awm->create_menu("Test[debug]", "g", AWM_ALL);
1420
1421	awm->insert_menu_topic("treestat", "Tree statistics", "s", "", AWM_ALL, TESTMENU_treeStats);
1422	awm->insert_menu_topic("setlens", "Set branchlens", "b", "", AWM_ALL, makeWindowCallback(TESTMENU_setBranchlen, ntw));
1423	awm->insert_menu_topic("sorttreebyname", "Sort tree by name", "o", "", AWM_ALL, makeWindowCallback(TESTMENU_sortTreeByName, ntw));
1424	}
1425	#endif // TESTMENU
1426
1427	static GB_ERROR pars_check_size(AW_root awr, GB_ERROR& warning, const adfiltercbstruct filterDef) {
1428	GB_ERROR error = NULp;
1429	warning = NULp;
1430
1431	char *tree_name = awr->awar(AWAR_TREE)->read_string();
1432	char *filter = awr->awar(filterDef->def_filter)->read_string();
1433	long ali_len = 0;
1434
1435	if (strlen(filter)) {
1436	int i;
1437	for (i=0; filter[i]; i++) {
1438	if (filter[i] != '0') ali_len++;
1439	}
1440	}
1441	else {
1442	char *ali_name = awr->awar(AWAR_ALIGNMENT)->read_string();
1443	ali_len = GBT_get_alignment_len(ap_main->get_gb_main(), ali_name);
1444	if (ali_len<=0) {
1445	error = "Please select a valid alignment";
1446	GB_clear_error();
1447	}
1448	free(ali_name);
1449	}
1450
1451	if (!error) {
1452	long tree_size = GBT_size_of_tree(ap_main->get_gb_main(), tree_name);
1453	if (tree_size == -1) {
1454	error = "Please select an existing tree";
1455	}
1456	else {
1457	size_t expected_memuse = (ali_len * tree_size * 4 / 1024);
1458	if (expected_memuse > GB_get_usable_memory()) {
1459	warning = GBS_global_string("Estimated memory usage (%s) exceeds physical memory (will swap)\n"
1460	"(did you specify a filter?)",
1461	GBS_readable_size(expected_memuse, "b"));
1462	}
1463	}
1464	}
1465
1466	free(filter);
1467	free(tree_name);
1468
1469	ap_assert(!GB_have_error());
1470	return error;
1471	}
1472
1473	static void pars_reset_optimal_parsimony(AW_window *aww) {
1474	AW_root *awr = aww->get_root();
1475	awr->awar(AWAR_BEST_PARSIMONY)->write_int(awr->awar(AWAR_PARSIMONY)->read_int());
1476	}
1477
1478	class LowDataCheck {
1479	int leafs; // counts leafs with insufficiant data
1480	int inner; // same for inner nodes
1481
1482	public:
1483	LowDataCheck() : leafs(0), inner(0) {}
1484
1485	void count(AP_tree_nlen *node);
1486
1487	int get_leafs() const { return leafs; }
1488	int get_inner() const { return inner; }
1489	};
1490
1491	void LowDataCheck::count(AP_tree_nlen *node) {
1492	const AP_combinableSeq *seq = node->get_seq();
1493	AP_FLOAT bases = seq->weighted_base_count();
1494
1495	if (node->is_leaf()) {
1496	if (bases<MIN_SEQUENCE_LENGTH) ++leafs;
1497	}
1498	else {
1499	if (bases<MIN_SEQUENCE_LENGTH) ++inner;
1500
1501	count(node->get_leftson());
1502	count(node->get_rightson());
1503	}
1504	}
1505
1506	static void PARS_infomode_cb(UNFIXED, TREE_canvas *canvas, AWT_COMMAND_MODE mode) {
1507	AWT_trigger_remote_action(NULp, canvas->gb_main, "ARB_NT:species_info");
1508	nt_mode_event(NULp, canvas, mode);
1509	}
1510
1511	static void pars_start_cb(AW_window aw_parent, WeightedFilter wfilt, const PARS_commands *cmds) {
1512	AW_root *awr = aw_parent->get_root();
1513	GBDATA *gb_main = ap_main->get_gb_main();
1514	GB_begin_transaction(gb_main);
1515	{
1516	GB_ERROR warning;
1517	GB_ERROR error = pars_check_size(awr, warning, wfilt->get_adfiltercbstruct());
1518
1519	if (warning && !error) {
1520	char *question = GBS_global_string_copy("%s\nDo you want to continue?", warning);
1521	bool cont = aw_ask_sure("swap_warning", question);
1522	free(question);
1523
1524	if (!cont) error = "User abort";
1525
1526	}
1527
1528	if (error) {
1529	aw_message(error);
1530	GB_commit_transaction(gb_main);
1531	return;
1532	}
1533	}
1534
1535
1536	AW_window_menu_modes *awm = new AW_window_menu_modes;
1537	awm->init(awr, "ARB_PARSIMONY", "ARB_PARSIMONY", 400, 200);
1538
1539	GLOBAL_PARS->generate_tree(wfilt);
1540
1541	TREE_canvas *ntw;
1542	{
1543	AP_tree_display_style prev_style = global_tree()->get_tree_style();
1544	global_tree()->set_tree_style(AP_LIST_SIMPLE, NULp); // avoid NDS warnings during startup
1545	ntw = new TREE_canvas(gb_main, awm, awm->get_window_id(), global_tree(), awr->awar(AWAR_TREE));
1546	global_tree()->set_tree_style(prev_style, ntw);
1547	}
1548
1549	{
1550	GB_ERROR error = NULp;
1551	arb_progress progress("loading tree");
1552	NT_reload_tree_event(awr, ntw, false); // load tree (but do not expose - first zombies need to be removed)
1553	if (!global_tree()->get_root_node()) {
1554	error = "Failed to load the selected tree";
1555	}
1556	else {
1557	AP_tree_edge::initialize(rootNode()); // builds edges
1558	long removed = global_tree_root()->remove_leafs(AWT_REMOVE_ZOMBIES);
1559
1560	PARS_tree_init(global_tree());
1561	removed += global_tree_root()->remove_leafs(AWT_RemoveType(AWT_REMOVE_ZOMBIES \| AWT_REMOVE_NO_SEQUENCE));
1562
1563	if (!global_tree()->get_root_node()) {
1564	const char *aliname = global_tree_root()->get_aliview()->get_aliname();
1565	error = GBS_global_string("Less than 2 species contain data in '%s'\n"
1566	"Tree vanished", aliname);
1567	}
1568	else if (removed) {
1569	aw_message(GBS_global_string("Removed %li leafs (zombies or species w/o data in alignment)", removed));
1570	}
1571
1572	error = GB_end_transaction(ntw->gb_main, error);
1573	if (!error) {
1574	progress.subtitle("Calculating inner nodes");
1575	GLOBAL_PARS->get_root_node()->costs();
1576
1577	progress.subtitle("Checking amount of data");
1578	LowDataCheck lowData;
1579	lowData.count(GLOBAL_PARS->get_root_node());
1580
1581	bool warned = false;
1582	if (lowData.get_inner()>0) {
1583	aw_message(GBS_global_string("Inner nodes with insufficient data: %i", lowData.get_inner()));
1584	warned = true;
1585	}
1586	if (lowData.get_leafs()>0) {
1587	aw_message(GBS_global_string("Species with insufficient data: %i", lowData.get_leafs()));
1588	warned = true;
1589	}
1590	#define _STRGIZE(x) #x
1591	#define INT2STR(i) _STRGIZE(i)
1592	if (warned) {
1593	aw_message("Warning: low sequence data (<" INT2STR(MIN_SEQUENCE_LENGTH) " bp) detected! (filter too restrictive?)");
1594	}
1595	}
1596	}
1597	if (error) aw_popup_exit(error);
1598	}
1599
1600	if (cmds->add_marked) NT_add_quick(NULp, ntw, NT_ADD_MARKED);
1601	if (cmds->add_selected) NT_add_quick(NULp, ntw, NT_ADD_SELECTED);
1602	if (cmds->calc_branch_lengths) NT_calc_branchlengths_reorder_and_update(awm, ntw);
1603	if (cmds->calc_bootstrap) NT_bootstrap(awm, ntw, 0);
1604	if (cmds->quit) pars_exit(awm);
1605
1606	GB_transaction ta(ntw->gb_main);
1607
1608	#if defined(DEBUG)
1609	AWT_create_debug_menu(awm);
1610	#endif // DEBUG
1611
1612	awm->create_menu("File", "F", AWM_ALL);
1613	{
1614	insert_macro_menu_entry(awm, false);
1615	awm->insert_menu_topic("print_tree", "Print Tree ...", "P", "tree2prt.hlp", AWM_ALL, makeWindowCallback(AWT_popup_print_window, static_cast<AWT_canvas*>(ntw)));
1616	awm->insert_menu_topic("quit", "Quit", "Q", "quit.hlp", AWM_ALL, pars_exit);
1617	}
1618
1619	awm->create_menu("Species", "S", AWM_ALL);
1620	{
1621	NT_insert_mark_submenus(awm, ntw, 0);
1622
1623	}
1624	awm->create_menu("Tree", "T", AWM_ALL);
1625	{
1626
1627	awm->insert_menu_topic("nds", "NDS (Node Display Setup) ...", "N", "props_nds.hlp", AWM_ALL, makeCreateWindowCallback(AWT_create_nds_window, ntw->gb_main));
1628
1629	awm->sep______________();
1630	awm->insert_menu_topic("tree_print", "Print tree ...", "P", "tree2prt.hlp", AWM_ALL, makeWindowCallback(AWT_popup_print_window, static_cast<AWT_canvas*>(ntw)));
1631	awm->insert_menu_topic("tree_2_xfig", "Export tree to XFIG ...", "F", "tree2file.hlp", AWM_ALL, makeWindowCallback(AWT_popup_tree_export_window, static_cast<AWT_canvas*>(ntw)));
1632	awm->sep______________();
1633	NT_insert_collapse_submenu(awm, ntw);
1634	awm->sep______________();
1635	awm->insert_sub_menu("Remove Species from Tree", "R");
1636	{
1637	awm->insert_menu_topic("tree_remove_deleted", "Remove Zombies", "Z", "trm_del.hlp", AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_REMOVE_ZOMBIES));
1638	awm->insert_menu_topic("tree_remove_marked", "Remove Marked", "M", "trm_mrkd.hlp", AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_REMOVE_MARKED));
1639	awm->insert_menu_topic("tree_keep_marked", "Keep Marked", "K", "tkeep_mrkd.hlp", AWM_ALL, makeWindowCallback(NT_remove_leafs, ntw, AWT_KEEP_MARKED));
1640	}
1641	awm->close_sub_menu();
1642	awm->insert_sub_menu("Add Species to Tree", "A");
1643	{
1644	awm->insert_menu_topic("add_marked", "Add Marked Species", "M", "pa_quick.hlp", AWM_ALL, makeWindowCallback(NT_add_quick, ntw, NT_ADD_MARKED));
1645	awm->insert_menu_topic("add_marked_nni", "Add Marked Species + Local Optimization (NNI)", "N", "pa_add.hlp", AWM_ALL, makeWindowCallback(NT_add_and_NNI, ntw, NT_ADD_MARKED));
1646	awm->insert_menu_topic("rm_add_marked", "Remove & Add Marked Species", "R", "pa_quick.hlp", AWM_ALL, makeWindowCallback(NT_reAdd_quick, ntw, NT_ADD_MARKED));
1647	awm->insert_menu_topic("rm_add_marked_nni\|", "Remove & Add Marked + Local Optimization (NNI)", "L", "pa_add.hlp", AWM_ALL, makeWindowCallback(NT_reAdd_and_NNI, ntw, NT_ADD_MARKED));
1648	awm->sep______________();
1649	awm->insert_menu_topic("add_marked_partial", "Add Marked Partial Species", "P", "pa_partial.hlp", AWM_ALL, makeWindowCallback(NT_add_partial_and_update, ntw));
1650	awm->sep______________();
1651	awm->insert_menu_topic("add_selected", "Add Selected Species", "S", "pa_quick.hlp", AWM_ALL, makeWindowCallback(NT_add_quick, ntw, NT_ADD_SELECTED));
1652	awm->insert_menu_topic("add_selected_nni", "Add Selected Species + Local Optimization (NNI)", "O", "pa_add.hlp", AWM_ALL, makeWindowCallback(NT_add_and_NNI, ntw, NT_ADD_SELECTED));
1653	}
1654	awm->close_sub_menu();
1655	awm->sep______________();
1656	awm->insert_menu_topic("optimize", "Tree Optimization ...", "O", "pa_optimizer.hlp", AWM_ALL, makeCreateWindowCallback(createOptimizeWindow, ntw));
1657	awm->insert_menu_topic("reset", "Reset optimal parsimony", "s", "pa_reset.hlp", AWM_ALL, pars_reset_optimal_parsimony);
1658	awm->sep______________();
1659	awm->insert_menu_topic("beautify_tree", "Beautify Tree", "B", "resorttree.hlp", AWM_ALL, makeWindowCallback(NT_resort_tree_cb, ntw, BIG_BRANCHES_TO_TOP));
1660	awm->insert_menu_topic("calc_branch_lengths", "Calculate Branch Lengths", "L", "pa_branchlengths.hlp", AWM_ALL, makeWindowCallback(NT_calc_branchlengths_reorder_and_update, ntw));
1661	awm->sep______________();
1662	awm->insert_menu_topic("calc_upper_bootstrap_indep", "Calculate Upper Bootstrap Limit (dependent NNI)", "U", "pa_bootstrap.hlp", AWM_ALL, makeWindowCallback(NT_bootstrap, ntw, false));
1663	awm->insert_menu_topic("calc_upper_bootstrap_dep", "Calculate Upper Bootstrap Limit (independent NNI)", "i", "pa_bootstrap.hlp", AWM_ALL, makeWindowCallback(NT_bootstrap, ntw, true));
1664	awm->insert_menu_topic("tree_remove_remark", "Remove bootstrap values", "v", "trm_boot.hlp", AWM_ALL, makeWindowCallback(NT_remove_bootstrap, ntw));
1665	}
1666
1667	#if defined(TESTMENU)
1668	init_TEST_menu(awm, ntw);
1669	#endif // TESTMENU
1670
1671	awm->create_menu("Reset", "R", AWM_ALL);
1672	{
1673	awm->insert_menu_topic("reset_logical_zoom", "Logical Zoom", "L", "rst_log_zoom.hlp", AWM_ALL, makeWindowCallback(NT_reset_lzoom_cb, ntw));
1674	awm->insert_menu_topic("reset_physical_zoom", "Physical Zoom", "P", "rst_phys_zoom.hlp", AWM_ALL, makeWindowCallback(NT_reset_pzoom_cb, ntw));
1675	}
1676
1677	awm->create_menu("Properties", "P", AWM_ALL);
1678	{
1679	awm->insert_menu_topic("props_menu", "Frame settings ...", "F", "props_frame.hlp", AWM_ALL, AW_preset_window);
1680	awm->insert_menu_topic("props_tree2", "Tree options", "o", "nt_tree_settings.hlp", AWM_ALL, TREE_create_settings_window);
1681	awm->insert_menu_topic("props_tree", "Tree colors & fonts", "c", "color_props.hlp", AWM_ALL, makeCreateWindowCallback(AW_create_gc_window, ntw->gc_manager));
1682	awm->insert_menu_topic("props_kl", "Optimizer settings (KL)", "K", "kernlin.hlp", AWM_ALL, makeCreateWindowCallback(create_kernighan_properties_window));
1683	awm->sep______________();
1684	AW_insert_common_property_menu_entries(awm);
1685	awm->sep______________();
1686	awm->insert_menu_topic("save_props", "Save Defaults (pars.arb)", "D", "savedef.hlp", AWM_ALL, AW_save_properties);
1687	}
1688	awm->button_length(5);
1689
1690	awm->insert_help_topic("ARB_PARSIMONY help", "P", "arb_pars.hlp", AWM_ALL, makeHelpCallback("arb_pars.hlp"));
1691
1692	// ----------------------
1693	// mode buttons
1694	//
1695	// keep them synchronized as far as possible with those in ARB_PARSIMONY
1696	// see ../NTREE/NT_extern.cxx@keepModesSynchronized
1697
1698	awm->create_mode("mode_select.xpm", "mode_select.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SELECT));
1699	awm->create_mode("mode_mark.xpm", "mode_mark.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_MARK));
1700	awm->create_mode("mode_group.xpm", "mode_group.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_GROUP));
1701	awm->create_mode("mode_zoom.xpm", "mode_pzoom.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_ZOOM));
1702	awm->create_mode("mode_lzoom.xpm", "mode_lzoom.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_LZOOM));
1703
1704	awm->create_mode("mode_info.xpm", "mode_info.hlp", AWM_ALL, makeWindowCallback(PARS_infomode_cb, ntw, AWT_MODE_INFO));
1705	// reserve mode-locations (to put the modes below at the same position as in ARB_NT)
1706	awm->create_mode("mode_empty.xpm", "mode.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_EMPTY));
1707
1708	// topology-modification-modes
1709	awm->create_mode("mode_setroot.xpm", "mode_setroot.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SETROOT));
1710	awm->create_mode("mode_swap.xpm", "mode_swap.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_SWAP));
1711	awm->create_mode("mode_move.xpm", "mode_move.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_MOVE));
1712
1713	awm->create_mode("mode_nni.xpm", "mode_nni.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_NNI));
1714	awm->create_mode("mode_kernlin.xpm", "mode_kernlin.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_KERNINGHAN));
1715	awm->create_mode("mode_optimize.xpm", "mode_optimize.hlp", AWM_ALL, makeWindowCallback(nt_mode_event, ntw, AWT_MODE_OPTIMIZE));
1716
1717	awm->at(5, 2);
1718	awm->auto_space(0, -2);
1719	awm->shadow_width(1);
1720
1721
1722	int db_treex, db_treey;
1723	awm->get_at_position(&db_treex, &db_treey);
1724	awm->callback(makeHelpCallback("nt_tree_select.hlp"));
1725	awm->button_length(16);
1726	awm->help_text("nt_tree_select.hlp");
1727	awm->create_button(NULp, AWAR_TREE);
1728
1729
1730	int db_stackx, db_stacky;
1731	awm->label_length(8);
1732	awm->label("Stored");
1733	awm->get_at_position(&db_stackx, &db_stacky);
1734	awm->button_length(6);
1735	awm->callback(makeHelpCallback("ap_stack.hlp"));
1736	awm->help_text("ap_stack.hlp");
1737	awm->create_button(NULp, AWAR_STACKPOINTER);
1738
1739	int db_parsx, db_parsy;
1740	awm->label_length(14);
1741	awm->label("Current Pars:");
1742	awm->get_at_position(&db_parsx, &db_parsy);
1743
1744	awm->button_length(10);
1745	awm->create_button(NULp, AWAR_PARSIMONY, NULp, "+");
1746
1747	awm->button_length(0);
1748
1749	awm->callback(makeWindowCallback(NT_jump_cb, ntw, AP_JUMP_BY_BUTTON));
1750	awm->help_text("tr_jump.hlp");
1751	awm->create_button("JUMP", "Jump");
1752
1753	awm->callback(makeHelpCallback("arb_pars.hlp"));
1754	awm->help_text("help.hlp");
1755	awm->create_button("HELP", "HELP", "H");
1756
1757	awm->at_newline();
1758
1759	awm->button_length(8);
1760	awm->at_x(db_stackx);
1761	awm->callback(makeWindowCallback(AP_user_pop_cb, ntw));
1762	awm->help_text("ap_stack.hlp");
1763	awm->create_button("POP", "RESTORE");
1764
1765	awm->button_length(6);
1766	awm->callback(AP_user_push_cb);
1767	awm->help_text("ap_stack.hlp");
1768	awm->create_button("PUSH", "STORE");
1769
1770	awm->at_x(db_parsx);
1771	awm->label_length(14);
1772	awm->label("Optimal Pars:");
1773
1774	awm->button_length(10);
1775	awm->create_button(NULp, AWAR_BEST_PARSIMONY, NULp, "+");
1776
1777	awm->button_length(0);
1778	awm->auto_space(0, -2);
1779
1780	awm->at_x(db_treex);
1781	awm->callback(makeWindowCallback(NT_set_tree_style, ntw, AP_TREE_RADIAL));
1782	awm->help_text("tr_type_radial.hlp");
1783	awm->create_button("RADIAL_TREE", "#radial.xpm");
1784
1785	awm->callback(makeWindowCallback(NT_set_tree_style, ntw, AP_TREE_NORMAL));
1786	awm->help_text("tr_type_list.hlp");
1787	awm->create_button("LIST_TREE", "#dendro.xpm");
1788
1789	awm->at_newline();
1790	awm->at(db_treex, awm->get_at_yposition());
1791
1792	{
1793	SmartPtr<AW_at_storage> maxSize(AW_at_storage::make(awm, AW_AT_MAXSIZE));
1794
1795	awm->button_length(AWAR_FOOTER_MAX_LEN);
1796	awm->create_button(NULp, AWAR_FOOTER);
1797	awm->at_newline();
1798	awm->restore_at_from(*maxSize);
1799	}
1800
1801	awm->get_at_position(&db_treex, &db_treey);
1802	awm->set_info_area_height(db_treey);
1803
1804	awm->set_bottom_area_height(0);
1805
1806	aw_parent->hide(); // hide parent
1807	awm->show();
1808
1809	TREE_install_update_callbacks(ntw);
1810
1811	update_random_repeat(awr, AWT_TREE_PARS(ntw));
1812	AP_user_push_cb(aw_parent); // push initial tree
1813	set_keep_ghostnodes(); // make sure no stacked nodes get deleted
1814	}
1815
1816	static AW_window create_pars_init_window(AW_root awr, const PARS_commands *cmds) {
1817	AW_window_simple *aws = new AW_window_simple;
1818	aws->init(awr, "PARS_PROPS", "SET PARSIMONY OPTIONS");
1819	aws->load_xfig("pars/init.fig");
1820
1821	aws->button_length(10);
1822	aws->label_length(10);
1823
1824	aws->callback(pars_exit);
1825	aws->at("close");
1826	aws->create_button("ABORT", "ABORT", "A");
1827
1828	aws->callback(makeHelpCallback("arb_pars_init.hlp"));
1829	aws->at("help");
1830	aws->create_button("HELP", "HELP", "H");
1831
1832	GBDATA *gb_main = ap_main->get_gb_main();
1833	WeightedFilter *weighted_filter = // do NOT free (bound to callbacks)
1834	new WeightedFilter(gb_main, aws->get_root(), AWAR_FILTER_NAME, AWAR_COLUMNSTAT_NAME, aws->get_root()->awar_string(AWAR_ALIGNMENT));
1835
1836	aws->at("filter");
1837	aws->callback(makeCreateWindowCallback(awt_create_select_filter_win, weighted_filter->get_adfiltercbstruct()));
1838	aws->create_button("SELECT_FILTER", AWAR_FILTER_NAME);
1839
1840	aws->at("weights");
1841	aws->callback(makeCreateWindowCallback(COLSTAT_create_selection_window, weighted_filter->get_column_stat()));
1842	aws->sens_mask(AWM_EXP);
1843	aws->create_button("SELECT_CSP", AWAR_COLUMNSTAT_NAME);
1844	aws->sens_mask(AWM_ALL);
1845
1846	aws->at("alignment");
1847	awt_create_ALI_selection_list(gb_main, aws, AWAR_ALIGNMENT, "*=");
1848
1849	aws->at("tree");
1850	awt_create_TREE_selection_list(gb_main, aws, AWAR_TREE, true);
1851
1852	aws->callback(makeWindowCallback(pars_start_cb, weighted_filter, cmds));
1853	aws->at("go");
1854	aws->create_button("GO", "GO", "G");
1855
1856	return aws;
1857	}
1858
1859	KL_Settings::KL_Settings(AW_root *aw_root) {
1860	maxdepth = aw_root->awar(AWAR_KL_MAXDEPTH)->read_int();
1861	incdepth = aw_root->awar(AWAR_KL_INCDEPTH)->read_int();
1862
1863	Static.enabled = aw_root->awar(AWAR_KL_STATIC_ENABLED)->read_int();
1864	Static.depth[0] = 2; // always test both possibilities at starting edge
1865	Static.depth[1] = aw_root->awar(AWAR_KL_STATIC_DEPTH1)->read_int();
1866	Static.depth[2] = aw_root->awar(AWAR_KL_STATIC_DEPTH2)->read_int();
1867	Static.depth[3] = aw_root->awar(AWAR_KL_STATIC_DEPTH3)->read_int();
1868	Static.depth[4] = aw_root->awar(AWAR_KL_STATIC_DEPTH4)->read_int();
1869	Static.depth[5] = aw_root->awar(AWAR_KL_STATIC_DEPTH5)->read_int();
1870
1871	Dynamic.enabled = aw_root->awar(AWAR_KL_DYNAMIC_ENABLED)->read_int();
1872	Dynamic.start = aw_root->awar(AWAR_KL_DYNAMIC_START)->read_int();
1873	Dynamic.maxx = aw_root->awar(AWAR_KL_DYNAMIC_MAXX)->read_int();
1874	Dynamic.maxy = aw_root->awar(AWAR_KL_DYNAMIC_MAXY)->read_int();
1875	Dynamic.type = (KL_DYNAMIC_THRESHOLD_TYPE)aw_root->awar(AWAR_KL_FUNCTION_TYPE)->read_int();
1876
1877	whichEdges = ANY_EDGE;
1878	if (aw_root->awar(AWAR_OPTI_MARKED_ONLY)->read_int()) whichEdges = EdgeSpec(whichEdges\|SKIP_UNMARKED_EDGES);
1879	if (aw_root->awar(AWAR_OPTI_SKIP_FOLDED)->read_int()) whichEdges = EdgeSpec(whichEdges\|SKIP_FOLDED_EDGES);
1880	}
1881	#if defined(UNIT_TESTS)
1882	KL_Settings::KL_Settings() {
1883	// set default values
1884	maxdepth = 15;
1885
1886	Static.enabled = true;
1887	Static.depth[0] = 2; // always test both possibilities at starting edge
1888	Static.depth[1] = 8;
1889	Static.depth[2] = 6;
1890	Static.depth[3] = 6;
1891	Static.depth[4] = 6;
1892	Static.depth[5] = 6;
1893
1894	Dynamic.enabled = true;
1895	Dynamic.start = 100;
1896	Dynamic.maxy = 150;
1897	Dynamic.maxx = 6;
1898
1899	// these values do not seem to have any effect (i.e. are not covered by unit-tests):
1900	incdepth = 4;
1901
1902	// const setting (not configurable)
1903	Dynamic.type = AP_QUADRAT_START;
1904	whichEdges = EdgeSpec(SKIP_UNMARKED_EDGES\|SKIP_FOLDED_EDGES);
1905	}
1906	#endif
1907
1908	static void create_optimize_vars(AW_root *aw_root, AW_default props) {
1909	// kernighan
1910
1911	aw_root->awar_int(AWAR_OPTI_MARKED_ONLY, 1, props);
1912	aw_root->awar_int(AWAR_OPTI_SKIP_FOLDED, 1, props);
1913
1914	aw_root->awar_int(AWAR_KL_MAXDEPTH, 15, props);
1915	aw_root->awar_int(AWAR_KL_INCDEPTH, 4, props);
1916
1917	aw_root->awar_int(AWAR_KL_STATIC_ENABLED, 1, props);
1918	aw_root->awar_int(AWAR_KL_STATIC_DEPTH1, 5, props)->set_minmax(1, 8);
1919	aw_root->awar_int(AWAR_KL_STATIC_DEPTH2, 3, props)->set_minmax(1, 6);
1920	aw_root->awar_int(AWAR_KL_STATIC_DEPTH3, 2, props)->set_minmax(1, 6);
1921	aw_root->awar_int(AWAR_KL_STATIC_DEPTH4, 2, props)->set_minmax(1, 6);
1922	aw_root->awar_int(AWAR_KL_STATIC_DEPTH5, 1, props)->set_minmax(1, 6);
1923
1924	aw_root->awar_int(AWAR_KL_DYNAMIC_ENABLED, 1, props);
1925	aw_root->awar_int(AWAR_KL_DYNAMIC_START, 100, props);
1926	aw_root->awar_int(AWAR_KL_DYNAMIC_MAXX, 6, props);
1927	aw_root->awar_int(AWAR_KL_DYNAMIC_MAXY, 150, props);
1928
1929	aw_root->awar_int(AWAR_KL_FUNCTION_TYPE, AP_QUADRAT_START, props);
1930	}
1931
1932	static void pars_create_all_awars(AW_root awr, AW_default aw_def, GBDATA gb_main) {
1933	awr->awar_string(AWAR_SPECIES_NAME, "", gb_main);
1934	awr->awar_string(AWAR_FOOTER, "", aw_def);
1935
1936	// copy currently selected alignment to awar:
1937	{
1938	GB_transaction ta(gb_main);
1939	char *dali = GBT_get_default_alignment(gb_main);
1940	awr->awar_string(AWAR_ALIGNMENT, dali, gb_main)->write_string(dali);
1941	free(dali);
1942	}
1943
1944	awt_create_filter_awars(awr, aw_def, AWAR_FILTER_NAME, AWAR_ALIGNMENT);
1945	awt_set_awar_to_valid_filter_good_for_tree_methods(gb_main, awr, AWAR_FILTER_NAME);
1946
1947	awr->awar_int(AWAR_PARS_TYPE, PARS_WAGNER, gb_main);
1948
1949	{
1950	GB_transaction ta(gb_main);
1951	GBDATA *gb_tree_name = GB_search(gb_main, AWAR_TREE, GB_STRING);
1952	char *tree_name = GB_read_string(gb_tree_name);
1953
1954	awr->awar_string(AWAR_TREE, "", aw_def)->write_string(tree_name);
1955	free(tree_name);
1956	}
1957
1958	awr->awar_int(AWAR_PARSIMONY, 0, aw_def);
1959	awr->awar_int(AWAR_BEST_PARSIMONY, 0, aw_def);
1960	awr->awar_int(AWAR_STACKPOINTER, 0, aw_def);
1961
1962	awr->awar_int(AWAR_RAND_REPEAT, 1, aw_def)->set_minmax(1, 1000000); // default value is overwritten by update_random_repeat()
1963	awr->awar_int(AWAR_RAND_PERCENT, 50, aw_def)->set_minmax(1, 100);
1964
1965	create_optimize_vars(awr, aw_def);
1966	create_nds_vars(awr, aw_def, gb_main, false);
1967
1968	TREE_create_awars(awr, gb_main);
1969
1970	#if defined(DEBUG)
1971	AWT_create_db_browser_awars(awr, aw_def);
1972	#endif // DEBUG
1973
1974	GB_ERROR error = ARB_init_global_awars(awr, aw_def, gb_main);
1975	if (error) aw_message(error);
1976	}
1977
1978	static AW_root *AD_map_viewer_aw_root = NULp;
1979
1980	void PARS_map_viewer(GBDATA *gb_species, AD_MAP_VIEWER_TYPE vtype) {
1981	// Note: sync with ../NTREE/ad_spec.cxx@launch_MapViewer_cb
1982
1983	if (AD_map_viewer_aw_root &&
1984	gb_species &&
1985	(vtype == ADMVT_SELECT \|\| vtype == ADMVT_INFO))
1986	{
1987	AD_map_viewer_aw_root->awar(AWAR_SPECIES_NAME)->write_string(GBT_read_name(gb_species));
1988	}
1989	}
1990
1991	int ARB_main(int argc, char *argv[]) {
1992	aw_initstatus();
1993
1994	GB_shell shell;
1995	AW_root *aw_root = AWT_create_root("pars.arb", "ARB_PARS", need_macro_ability(), &argc, &argv);
1996	AD_map_viewer_aw_root = aw_root;
1997
1998	ap_main = new AP_main;
1999	GLOBAL_PARS = new ArbParsimony();
2000
2001	const char *db_server = ":";
2002
2003	PARS_commands cmds;
2004
2005	while (argc>=2 && argv[1][0] == '-') {
2006	argc--;
2007	argv++;
2008	if (!strcmp(argv[0], "-quit")) cmds.quit = 1;
2009	else if (!strcmp(argv[0], "-add_marked")) cmds.add_marked = 1;
2010	else if (!strcmp(argv[0], "-add_selected")) cmds.add_selected = 1;
2011	else if (!strcmp(argv[0], "-calc_branchlengths")) cmds.calc_branch_lengths = 1;
2012	else if (!strcmp(argv[0], "-calc_bootstrap")) cmds.calc_bootstrap = 1;
2013	else {
2014	fprintf(stderr, "Unknown option '%s'\n", argv[0]);
2015
2016	printf(" Options: Meaning:\n"
2017	"\n"
2018	" -add_marked add marked species (without changing topology)\n"
2019	" -add_selected add selected species (without changing topology)\n"
2020	" -calc_branchlengths calculate branch lengths only\n"
2021	" -calc_bootstrap estimate bootstrap values\n"
2022	" -quit quit after performing operations\n"
2023	);
2024
2025	exit(EXIT_FAILURE);
2026	}
2027	}
2028
2029
2030	if (argc==2) db_server = argv[1];
2031
2032	GB_ERROR error = ap_main->open(db_server);
2033	if (!error) {
2034	GBDATA *gb_main = ap_main->get_gb_main();
2035	error = configure_macro_recording(aw_root, "ARB_PARS", gb_main);
2036
2037	if (!error) {
2038	#if defined(DEBUG)
2039	AWT_announce_db_to_browser(gb_main, GBS_global_string("ARB-database (%s)", db_server));
2040	#endif // DEBUG
2041
2042	pars_create_all_awars(aw_root, AW_ROOT_DEFAULT, gb_main);
2043
2044	AW_window *aww = create_pars_init_window(aw_root, &cmds);
2045	aww->show();
2046
2047	AWT_install_cb_guards();
2048	aw_root->main_loop();
2049	}
2050	}
2051
2052	if (error) aw_popup_exit(error);
2053	return EXIT_SUCCESS;
2054	}
2055
2056
2057	// --------------------------------------------------------------------------------
2058
2059	#ifdef UNIT_TESTS
2060	#include <arb_file.h>
2061	#include <arb_diff.h>
2062	#include <test_unit.h>
2063	#include <AP_seq_dna.hxx>
2064	#include <AP_seq_protein.hxx>
2065	#include "test_env.h"
2066
2067	// #define AUTO_UPDATE_IF_CHANGED // uncomment to auto update expected results
2068
2069	static arb_test::match_expectation topologyEquals(AP_tree_nlen root_node, const char treefile_base) {
2070	using namespace arb_test;
2071	expectation_group fulfilled;
2072
2073	char *outfile = GBS_global_string_copy("pars/%s.tree", treefile_base);
2074	char *expected = GBS_global_string_copy("%s.expected", outfile);
2075	bool update = false;
2076
2077	{
2078	FILE *out = fopen(outfile, "wt");
2079	fulfilled.add(that(out).does_differ_from_NULL());
2080	char *newick = GBT_tree_2_newick(root_node, NewickFormat(nLENGTH\|nWRAP), false);
2081	fputs(newick, out);
2082	free(newick);
2083	fclose(out);
2084	}
2085
2086	if (GB_is_regularfile(expected)) {
2087	bool match_exp_topo = textfiles_have_difflines(outfile,expected,0);
2088	#if defined(AUTO_UPDATE_IF_CHANGED)
2089	if (!match_exp_topo) update = true;
2090	#endif
2091	if (!update) fulfilled.add(that(match_exp_topo).is_equal_to(true));
2092	}
2093	else {
2094	update = true;
2095	}
2096
2097	if (update) TEST_COPY_FILE(outfile, expected);
2098	TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(outfile));
2099
2100	free(expected);
2101	free(outfile);
2102
2103	return all().ofgroup(fulfilled);
2104	}
2105
2106	template<class ENV>
2107	arb_test::match_expectation calcCostsCausesCombines(ENV& env, AP_FLOAT exp_pars, long exp_combines) {
2108	using namespace arb_test;
2109	expectation_group fulfilled;
2110
2111	long combines_b4_costCalc = env.combines_performed();
2112	fulfilled.add(that(combines_b4_costCalc).is_equal_to(0));
2113
2114	AP_FLOAT new_pars = env.root_node()->costs();
2115	long combines_by_costCalc = env.combines_performed();
2116
2117	fulfilled.add(that(new_pars).fulfills(epsilon_similar(0.001), exp_pars));
2118	fulfilled.add(that(combines_by_costCalc).is_equal_to(exp_combines));
2119
2120	return all().ofgroup(fulfilled);
2121	}
2122
2123	#define TEST_EXPECT_SAVED_TOPOLOGY(env,exp_topo) TEST_EXPECTATION(topologyEquals(env.root_node(), exp_topo))
2124	#define TEST_EXPECT_SAVED_TOPOLOGY__BROKEN(env,exp_topo,got_topo) TEST_EXPECTATION__BROKEN(topologyEquals(env.root_node(), exp_topo), topologyEquals(env.root_node(), got_topo))
2125
2126	#define TEST_EXPECT_PARSVAL(env,exp_pars) TEST_EXPECT_EQUAL(env.root_node()->costs(), exp_pars);
2127	#define TEST_EXPECT_ONLY_PARSVAL_COMBINES(env,exp_pars,exp_combines) TEST_EXPECTATION(calcCostsCausesCombines(env, exp_pars, exp_combines))
2128	// use TEST_EXPECT_ONLY_PARSVAL_COMBINES when
2129	// - no combines occurred (or combines were just tested using TEST_EXPECT_COMBINES_PERFORMED) and
2130	// - topology was modified, so that calculation of costs causes new combines.
2131	#define TEST_EXPECT_KNOWN_PARSVAL(env,exp_pars) TEST_EXPECT_ONLY_PARSVAL_COMBINES(env,exp_pars,0)
2132
2133	enum TopoMod {
2134	MOD_REMOVE_MARKED,
2135
2136	MOD_QUICK_READD,
2137	MOD_QUICK_ADD,
2138	MOD_READD_NNI,
2139
2140	MOD_ADD_PARTIAL,
2141
2142	MOD_CALC_LENS,
2143	MOD_OPTI_NNI,
2144	MOD_OPTI_GLOBAL,
2145
2146	MOD_MIX_TREE,
2147	};
2148
2149	template <typename SEQ>
2150	static void modifyTopology(PARSIMONY_testenv<SEQ>& env, TopoMod mod) {
2151	switch (mod) {
2152	case MOD_REMOVE_MARKED:
2153	env.graphic_tree()->get_tree_root()->remove_leafs(AWT_REMOVE_MARKED);
2154	break;
2155
2156	case MOD_QUICK_READD:
2157	nt_reAdd(env.graphic_tree(), NT_ADD_MARKED, true);
2158	break;
2159
2160	case MOD_QUICK_ADD:
2161	nt_add(env.graphic_tree(), NT_ADD_MARKED, true);
2162	break;
2163
2164	case MOD_READD_NNI:
2165	nt_reAdd(env.graphic_tree(), NT_ADD_MARKED, false);
2166	break;
2167
2168	case MOD_ADD_PARTIAL:
2169	nt_add_partial(env.graphic_tree());
2170	break;
2171
2172	case MOD_CALC_LENS:
2173	calc_branchlengths_and_reorder(env.graphic_tree());
2174	break;
2175
2176	case MOD_OPTI_NNI: // only marked/unfolded
2177	recursiveNNI(env.graphic_tree(), env.get_KL_settings().whichEdges);
2178	break;
2179
2180	case MOD_OPTI_GLOBAL:
2181	optimizeTree(env.graphic_tree(), env.get_KL_settings());
2182	break;
2183
2184	case MOD_MIX_TREE: {
2185	long leafs = rootNode()->count_leafs();
2186	mixtree_and_calclengths(env.graphic_tree(), calculate_default_random_repeat(leafs), 100, ANY_EDGE);
2187	break;
2188	}
2189	}
2190	}
2191
2192	template <typename SEQ>
2193	static arb_test::match_expectation modifyingTopoResultsIn(TopoMod mod, const char *topo, long pars_expected, PARSIMONY_testenv<SEQ>& env, bool restore) {
2194	using namespace arb_test;
2195	expectation_group fulfilled;
2196
2197	TEST_EXPECT_VALID_TREE(env.root_node());
2198
2199	Level upc = env.get_user_push_counter();
2200	Level fl = env.get_frame_level();
2201
2202	if (restore) {
2203	env.push();
2204	TEST_EXPECT_VALID_TREE(env.root_node());
2205	}
2206
2207	AWT_graphic_exports& exports = env.graphic_tree()->exports;
2208	exports.clear_save_request();
2209	modifyTopology(env, mod);
2210	if (topo) {
2211	fulfilled.add(topologyEquals(env.root_node(), topo));
2212	if (mod != MOD_REMOVE_MARKED) { // remove_leafs doesn't request save
2213	fulfilled.add(that(exports.needs_save()).is_equal_to(true));
2214	}
2215	}
2216
2217	fulfilled.add(that(allBranchlengthsAreDefined(env.root_node())).is_equal_to(true));
2218
2219	if (pars_expected != -1) {
2220	fulfilled.add(that(env.root_node()->costs()).is_equal_to(pars_expected));
2221	}
2222
2223	if (restore) {
2224	TEST_EXPECT_VALID_TREE(env.root_node());
2225	TEST_VALIDITY(env.pop_will_produce_valid_tree());
2226	env.pop();
2227	bool blen_def_after_pop = allBranchlengthsAreDefined(env.root_node());
2228	fulfilled.add(that(blen_def_after_pop).is_equal_to(true));
2229	}
2230
2231	TEST_EXPECT_EQUAL(fl, env.get_frame_level());
2232	TEST_EXPECT_EQUAL(upc, env.get_user_push_counter());
2233
2234	TEST_EXPECT_VALID_TREE(env.root_node());
2235
2236	return all().ofgroup(fulfilled);
2237	}
2238
2239	static arb_test::match_expectation movingRootDoesntAffectCosts(long pars_expected) {
2240	using namespace arb_test;
2241	expectation_group fulfilled;
2242
2243	long pars_min = LONG_MAX;
2244	long pars_max = LONG_MIN;
2245
2246	for (int depth_first = 0; depth_first<=1; ++depth_first) {
2247	for (int push_local = 0; push_local<=1; ++push_local) {
2248	EdgeChain chain(rootEdge(), ANY_EDGE, depth_first);
2249
2250	if (!push_local) ap_main->remember();
2251	while (chain) {
2252	AP_tree_edge edge = chain; ++chain;
2253
2254	if (push_local) ap_main->remember();
2255	edge->set_root();
2256	long pars = rootNode()->costs();
2257	pars_min = std::min(pars, pars_min);
2258	pars_max = std::max(pars, pars_max);
2259	if (push_local) ap_main->revert();
2260	}
2261	if (!push_local) ap_main->revert();
2262	}
2263	}
2264
2265	fulfilled.add(that(pars_min).is_equal_to(pars_expected));
2266	fulfilled.add(that(pars_max).is_equal_to(pars_expected));
2267
2268	return all().ofgroup(fulfilled);
2269	}
2270
2271	static GBDATA copy_to(GBDATA gb_species, const char *newShortname) {
2272	GBDATA *gb_species_data = GB_get_father(gb_species);
2273	GBDATA *gb_new_species = GB_create_container(gb_species_data, "species");
2274
2275	GB_ERROR error = NULp;
2276	if (!gb_new_species) {
2277	error = GB_await_error();
2278	}
2279	else {
2280	error = GB_copy_dropProtectMarksAndTempstate(gb_new_species, gb_species);
2281	if (!error) error = GBT_write_string(gb_new_species, "name", newShortname);
2282	}
2283
2284	ap_assert(contradicted(gb_new_species, error));
2285	return gb_new_species;
2286	}
2287
2288	inline void mark_only(GBDATA *gb_species) {
2289	GBDATA *gb_main = GB_get_root(gb_species);
2290	GB_transaction ta(gb_main);
2291	GBT_mark_all(gb_main, 0);
2292	GB_write_flag(gb_species, 1);
2293	}
2294	inline void mark(GBDATA *gb_species) {
2295	GBDATA *gb_main = GB_get_root(gb_species);
2296	GB_transaction ta(gb_main);
2297	GB_write_flag(gb_species, 1);
2298	}
2299	inline void mark_all(GBDATA *gb_main) {
2300	GB_transaction ta(gb_main);
2301	GBT_mark_all(gb_main, 1);
2302	}
2303
2304	inline int is_partial(GBDATA *gb_species) {
2305	GB_transaction ta(gb_species);
2306	return GBT_is_partial(gb_species, -1, false);
2307	}
2308
2309	template <typename SEQ>
2310	static arb_test::match_expectation addedAsBrotherOf(const char name, const char allowedBrothers, PARSIMONY_testenv<SEQ>& env) {
2311	using namespace arb_test;
2312	expectation_group fulfilled;
2313
2314	AP_tree_nlen *node_in_tree = env.root_node()->findLeafNamed(name);
2315	ap_assert(node_in_tree);
2316	fulfilled.add(that(node_in_tree).does_differ_from_NULL());
2317
2318	const char *brother = node_in_tree->get_brother()->name;
2319	ap_assert(brother);
2320	fulfilled.add(that(allowedBrothers).does_contain(brother));
2321
2322	return all().ofgroup(fulfilled);
2323	}
2324
2325	template <typename SEQ>
2326	static arb_test::match_expectation addingPartialResultsIn(GBDATA gb_added_species, const char allowedBrothers, const char *topo, int pars_expected, PARSIMONY_testenv<SEQ>& env) {
2327	using namespace arb_test;
2328	expectation_group fulfilled;
2329
2330	mark_only(gb_added_species);
2331	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2332	fulfilled.add(modifyingTopoResultsIn(MOD_ADD_PARTIAL, topo, pars_expected, env, false));
2333	fulfilled.add(that(is_partial(gb_added_species)).is_equal_to(1));
2334
2335	const char *name = GBT_read_name(gb_added_species);
2336	fulfilled.add(addedAsBrotherOf(name, allowedBrothers, env));
2337
2338	return all().ofgroup(fulfilled);
2339	}
2340
2341	static int seqDiff(GBDATA gb_main, const char aliname, const char species1, const char species2, int startPos, int endPos) {
2342	GB_transaction ta(gb_main);
2343
2344	GBDATA *gb_species1 = GBT_expect_species(gb_main, species1);
2345	GBDATA *gb_species2 = GBT_expect_species(gb_main, species2);
2346	int diffs = -1;
2347
2348	if (gb_species1 && gb_species2) {
2349	GBDATA *gb_seq1 = GBT_find_sequence(gb_species1, aliname);
2350	GBDATA *gb_seq2 = GBT_find_sequence(gb_species2, aliname);
2351
2352	if (gb_seq1 && gb_seq2) {
2353	char *seq1 = GB_read_string(gb_seq1);
2354	char *seq2 = GB_read_string(gb_seq2);
2355
2356	int maxPos1 = strlen(seq1)-1;
2357	#if defined(ASSERTION_USED)
2358	int maxPos2 = strlen(seq2)-1;
2359	#endif
2360	ap_assert(maxPos1 == maxPos2);
2361
2362	if (endPos>maxPos1) endPos = maxPos1;
2363
2364	diffs = 0;
2365	for (int p = startPos; p<=endPos; ++p) { // LOOP_ VECTORIZED[!<9.1] // @@@ fails in RELEASE code! // IRRELEVANT_LOOP
2366	diffs += seq1[p] != seq2[p];
2367	}
2368
2369	free(seq2);
2370	free(seq1);
2371	}
2372	}
2373
2374	return diffs;
2375	}
2376
2377	static GBDATA createPartialSeqFrom(GBDATA gb_main, const char aliname, const char dest_species, const char *source_species, int startPos, int endPos) {
2378	GB_transaction ta(gb_main);
2379
2380	GBDATA *gb_result = NULp;
2381	GBDATA *gb_source_species = GBT_expect_species(gb_main, source_species);
2382
2383	if (gb_source_species) {
2384	GBDATA *gb_dest_species = copy_to(gb_source_species, dest_species);
2385	GBDATA *gb_dest_seq = GBT_find_sequence(gb_dest_species, aliname); // =same as source seq
2386	char *seq = GB_read_string(gb_dest_seq);
2387
2388	if (seq) {
2389	int maxPos = strlen(seq)-1;
2390
2391	startPos = std::min(startPos, maxPos);
2392	endPos = std::min(endPos, maxPos);
2393
2394	if (startPos>0) memset(seq, '.', startPos);
2395	if (endPos<maxPos) memset(seq+endPos+1, '.', maxPos-endPos);
2396
2397	GB_ERROR error = GB_write_string(gb_dest_seq, seq);
2398	if (error) GB_export_error(error);
2399	else {
2400	gb_result = gb_dest_species; // success
2401	#if defined(DEBUG)
2402	fprintf(stderr, "created partial '%s' from '%s' (seq='%s')\n", dest_species, source_species, seq);
2403	#endif
2404	}
2405
2406	free(seq);
2407	}
2408	}
2409
2410	return gb_result;
2411	}
2412
2413	static GB_ERROR modifyOneBase(GBDATA gb_species, const char aliname, char cOld, char cNew) {
2414	GB_transaction ta(gb_species);
2415	GB_ERROR error = "failed to modifyOneBase";
2416
2417	GBDATA *gb_seq = GBT_find_sequence(gb_species, aliname);
2418	if (gb_seq) {
2419	char *seq = GB_read_string(gb_seq);
2420	if (seq) {
2421	char *B = strchr(seq, cOld);
2422	if (!B) {
2423	error = "does not contain base in modifyOneBase";
2424	}
2425	else {
2426	B[0] = cNew;
2427	error = GB_write_string(gb_seq, seq);
2428	}
2429	free(seq);
2430	}
2431	}
2432
2433	return error;
2434	}
2435
2436	static long unmark_unwanted(const char , long cd_gbd, void) {
2437	GBDATA gbd = (GBDATA)cd_gbd;
2438	GB_write_flag(gbd, 0);
2439	return 0;
2440	}
2441
2442	void TEST_nucl_tree_modifications() {
2443	const char *aliname = "ali_5s";
2444
2445	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
2446	TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
2447	TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
2448
2449	const int PARSIMONY_ORG = 302;
2450	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 14);
2451
2452	// [NUCOPTI] opposed to protein tests below the initial tree here is NOT optimized! compare .@PROTOPTI
2453	// -> removing and adding species produces a better tree (for add+NNI)
2454	//
2455	// diff initial->removed : http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-removed.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
2456	// diff initial->add-quick: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-add-quick.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
2457	// diff initial->add-NNI: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/nucl-add-NNI.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fnucl-initial.tree.expected
2458
2459	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, "nucl-removed", PARSIMONY_ORG-94, env, true)); // test remove-marked only (same code as part of nt_reAdd)
2460	TEST_EXPECT_COMBINES_PERFORMED(env, 3);
2461
2462	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "nucl-add-quick", PARSIMONY_ORG-18, env, true)); // test quick-add
2463	TEST_EXPECT_COMBINES_PERFORMED(env, 400);
2464
2465	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_READD_NNI, "nucl-add-NNI", PARSIMONY_ORG-20, env, true)); // test add + NNI
2466	TEST_EXPECT_COMBINES_PERFORMED(env, 503);
2467
2468	// test partial-add
2469	{
2470	GBDATA *gb_main = env.gbmain();
2471
2472	// create 2 non-overlapping partial species
2473	const int SPLIT = 55;
2474	GBDATA *CorGlutP = createPartialSeqFrom(gb_main, aliname, "CorGlutP", "CorGluta", 0, SPLIT);
2475	GBDATA *CloButyP = createPartialSeqFrom(gb_main, aliname, "CloButyP", "CloButyr", SPLIT+1, INT_MAX);
2476	GBDATA *CloButyM = createPartialSeqFrom(gb_main, aliname, "CloButyM", "CloButyr", SPLIT+1, INT_MAX);
2477	TEST_EXPECT_NO_ERROR(modifyOneBase(CloButyM, aliname, 'G', 'C')); // change first 'G' into 'C'
2478
2479	TEST_VALIDITY(env.all_available_pops_will_produce_valid_trees()); // no push yet (does nothing)
2480
2481	// test partials differ from full:
2482	TEST_REJECT_ZERO(seqDiff(gb_main, aliname, "CorGlutP", "CorGluta", 0, INT_MAX));
2483	TEST_REJECT_ZERO(seqDiff(gb_main, aliname, "CloButyP", "CloButyr", 0, INT_MAX));
2484	TEST_REJECT_ZERO(seqDiff(gb_main, aliname, "CloButyM", "CloButyr", 0, INT_MAX));
2485	// test partials created from CloButyr differ in partial range:
2486	TEST_REJECT_ZERO(seqDiff(gb_main, aliname, "CloButyM", "CloButyP", SPLIT+1, INT_MAX));
2487
2488	// test condition that "CloButyr and CloButy2 do NOT differ in seq-range of partial" (otherwise test below makes no sense!)
2489	TEST_EXPECT_ZERO(seqDiff(gb_main, aliname, "CloButyr", "CloButy2", SPLIT+1, INT_MAX));
2490
2491	// test that "CloButyr and CloButy2 DO differ in whole seq-range" (otherwise inserting into tree is non-deterministic)
2492	TEST_REJECT_ZERO(seqDiff(gb_main, aliname, "CloButyr", "CloButy2", 0, INT_MAX));
2493
2494	// add CloButyP (and undo)
2495	{
2496	env.push();
2497
2498	// CloButyr and CloButy2 do not differ in seq-range of partial -> any of both may be chosen as brother.
2499	// behavior should be changed with #605
2500	TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", "nucl-addPart-CloButyP", PARSIMONY_ORG, env));
2501	TEST_EXPECT_COMBINES_PERFORMED(env, 6);
2502	env.pop();
2503	}
2504
2505	{
2506	env.push();
2507	TEST_EXPECTATION(addingPartialResultsIn(CorGlutP, "CorGluta", "nucl-addPart-CorGlutP", PARSIMONY_ORG, env)); // add CorGlutP
2508	TEST_EXPECT_COMBINES_PERFORMED(env, 5); // @@@ partial-add should not perform combines at all (maybe caused by cost-recalc?)
2509	TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", "nucl-addPart-CorGlutP-CloButyP", PARSIMONY_ORG, env)); // also add CloButyP
2510	TEST_EXPECT_COMBINES_PERFORMED(env, 6);
2511	env.pop();
2512	}
2513
2514	// now add CorGlutP as full, then CloButyP and CloButyM as partials
2515	{
2516	env.push();
2517
2518	mark_only(CorGlutP);
2519	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2520	{
2521	GB_transaction ta(gb_main);
2522	TEST_EXPECT_NO_ERROR(GBT_write_int(CorGlutP, "ARB_partial", 0)); // revert species to "full"
2523	}
2524
2525	const int PARSIMONY_ADDED = PARSIMONY_ORG; // value after adding CorGlutP (as full-length sequence)
2526
2527	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "nucl-addPartialAsFull-CorGlutP", PARSIMONY_ADDED, env, false));
2528	TEST_EXPECT_COMBINES_PERFORMED(env, 230);
2529	TEST_EXPECT_EQUAL(is_partial(CorGlutP), 0); // check CorGlutP was added as full sequence
2530	TEST_EXPECTATION(addedAsBrotherOf("CorGlutP", "CorGluta", env)); // partial created from CorGluta gets inserted next to CorGluta
2531
2532	// add CloButyP as partial.
2533	// as expected it is placed next to matching full sequences (does not differ in partial range)
2534	TEST_EXPECTATION(addingPartialResultsIn(CloButyP, "CloButyr;CloButy2", NULp, PARSIMONY_ADDED, env));
2535	TEST_EXPECT_COMBINES_PERFORMED(env, 6);
2536
2537	// CloButyM differs slightly in overlap with CloButyr/CloButy2, but has no overlap with CorGlutP
2538	// shows bug described in #609 is fixed:
2539	TEST_EXPECTATION(addingPartialResultsIn(CloButyM, "CloButyP", "nucl-addPart-bug609",
2540	PARSIMONY_ADDED+1, // @@@ known bug - partial should not affect parsimony value; possibly related to ../HELP_SOURCE/oldhelp/pa_partial.hlp@WARNINGS
2541	env));
2542	TEST_EXPECT_COMBINES_PERFORMED(env, 7);
2543	env.pop();
2544	}
2545	}
2546
2547	TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
2548
2549	const int PARSIMONY_NNI_MARKED = PARSIMONY_ORG-18;
2550	const int PARSIMONY_NNI_ALL = PARSIMONY_ORG-18;
2551	const int PARSIMONY_OPTI_MARKED = PARSIMONY_ORG-25;
2552	const int PARSIMONY_OPTI_VISIBLE = PARSIMONY_ORG-26;
2553	const int PARSIMONY_OPTI_ALL = PARSIMONY_ORG-36;
2554
2555	{
2556	env.push();
2557	TEST_EXPECTATION(movingRootDoesntAffectCosts(PARSIMONY_ORG));
2558	TEST_EXPECT_COMBINES_PERFORMED(env, 342);
2559	env.pop();
2560	}
2561
2562	// ------------------------------
2563	// test optimize (some)
2564
2565	// mark initially marked species
2566	{
2567	GB_transaction ta(env.gbmain());
2568	GBT_restore_marked_species(env.gbmain(), "CorAquat;CorGluta;CurCitre;CloButyr;CloButy2;CytAquat");
2569	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2570	}
2571
2572	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_ORG);
2573
2574	// test branchlength calculation
2575	// (optimizations below implicitely recalculates branchlengths)
2576	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "nucl-calclength", PARSIMONY_ORG, env, false));
2577	TEST_EXPECT_COMBINES_PERFORMED(env, 120);
2578
2579	// test whether branchlength calculation depends on root-position
2580	{
2581	AP_tree_edge *orgRootEdge = rootEdge();
2582
2583	env.push();
2584
2585	const char *tested_roots[] = {
2586	"CloButyr",
2587	"CloTyro4",
2588	"CloTyrob",
2589	"CloInnoc",
2590	};
2591
2592	for (size_t r = 0; r<ARRAY_ELEMS(tested_roots); ++r) {
2593	const char *leafName = tested_roots[r];
2594	env.root_node()->findLeafNamed(leafName)->set_root();
2595	calc_branchlengths_and_reorder(env.graphic_tree());
2596	orgRootEdge->set_root();
2597	env.graphic_tree()->reorderTree(BIG_BRANCHES_TO_TOP);
2598
2599	TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-calclength");
2600	}
2601	TEST_EXPECT_COMBINES_PERFORMED(env, 517);
2602
2603	env.pop();
2604	}
2605
2606	// test optimize (some)
2607	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "nucl-opti-NNI", PARSIMONY_NNI_MARKED, env, true)); // test recursive NNI
2608	TEST_EXPECT_COMBINES_PERFORMED(env, 208);
2609
2610	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-marked-global", PARSIMONY_OPTI_MARKED, env, true)); // test recursive NNI+KL
2611	TEST_EXPECT_COMBINES_PERFORMED(env, 18518);
2612
2613	{
2614	KL_Settings& KL = env.get_KL_settings();
2615	LocallyModify<EdgeSpec> target(KL.whichEdges, EdgeSpec(KL.whichEdges&~SKIP_UNMARKED_EDGES)); // ignore marks; skip folded
2616
2617	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-visible-global", PARSIMONY_OPTI_VISIBLE, env, true)); // same result as if all species marked (see below)
2618	TEST_EXPECT_COMBINES_PERFORMED(env, 34925);
2619
2620	KL.whichEdges = EdgeSpec(KL.whichEdges&~SKIP_FOLDED_EDGES); // ignore marks and folding
2621
2622	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-global", PARSIMONY_OPTI_ALL, env, true)); // same result as if all species marked and all groups unfolded (see below)
2623	TEST_EXPECT_COMBINES_PERFORMED(env, 124811);
2624	}
2625
2626	// -----------------------------
2627	// test optimize (all)
2628
2629	// mark all species
2630	mark_all(env.gbmain());
2631	// unmark species not in tree
2632	{
2633	GB_transaction ta(env.gbmain());
2634	GB_HASH *markedNotInTree = GBT_create_marked_species_hash(env.gbmain());
2635	NT_remove_species_in_tree_from_hash(env.root_node(), markedNotInTree);
2636	GBS_hash_do_loop(markedNotInTree, unmark_unwanted, NULp);
2637	GBS_free_hash(markedNotInTree);
2638	}
2639	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2640	TEST_EXPECT_EQUAL(GBT_count_marked_species(env.gbmain()), 15);
2641
2642	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_ORG);
2643
2644	// test branchlength calculation
2645	// (optimizations below implicitely recalculates branchlengths)
2646	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "nucl-calclength", PARSIMONY_ORG, env, false));
2647	TEST_EXPECT_COMBINES_PERFORMED(env, 120);
2648
2649	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "nucl-opti-all-NNI", PARSIMONY_NNI_ALL, env, true)); // test recursive NNI
2650	TEST_EXPECT_COMBINES_PERFORMED(env, 242);
2651
2652	{
2653	env.push();
2654	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-visible-global", PARSIMONY_OPTI_VISIBLE, env, false)); // test recursive NNI+KL
2655	TEST_EXPECT_COMBINES_PERFORMED(env, 34925);
2656
2657	TEST_EXPECTATION(movingRootDoesntAffectCosts(PARSIMONY_OPTI_VISIBLE));
2658	TEST_EXPECT_COMBINES_PERFORMED(env, 336);
2659	env.pop();
2660	}
2661
2662	// unfold groups
2663	{
2664	AP_tree_nlen *CloTyrob = env.root_node()->findLeafNamed("CloTyrob");
2665	AP_tree_nlen *group = CloTyrob->get_father();
2666	ap_assert(group->gr.grouped);
2667	group->gr.grouped = false; // unfold the only folded group
2668
2669	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "nucl-opti-global", PARSIMONY_OPTI_ALL, env, true)); // test recursive NNI+KL
2670	TEST_EXPECT_COMBINES_PERFORMED(env, 124811);
2671	}
2672
2673	// test re-add all (i.e. test "create tree from scratch")
2674	// Note: trees generated below are NO LONGER better than optimized trees! (see also r13651)
2675
2676	// quick add:
2677	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "nucl-readdall-quick", PARSIMONY_ORG-7, env, true));
2678	TEST_EXPECT_COMBINES_PERFORMED(env, 439);
2679
2680	// quick add + NNI:
2681	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_READD_NNI, "nucl-readdall-NNI", PARSIMONY_ORG-8, env, true));
2682	TEST_EXPECT_COMBINES_PERFORMED(env, 613);
2683
2684	// test adding a too short sequence
2685	// (has to be last test, because it modifies seq data) << ------------ !!!!!
2686	{
2687	env.push();
2688
2689	AP_tree_nlen *CloTyrob = env.root_node()->findLeafNamed("CloTyrob");
2690	mark_only(CloTyrob->gb_node);
2691	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2692
2693	// modify sequence of CloTyrob (keep only some bases)
2694	{
2695	GB_transaction ta(env.gbmain());
2696	GBDATA *gb_seq = GBT_find_sequence(CloTyrob->gb_node, aliname);
2697
2698	char *seq = GB_read_string(gb_seq);
2699	int keep_bases = MIN_SEQUENCE_LENGTH-1;
2700
2701	for (int i = 0; seq[i]; ++i) {
2702	if (!GAP::is_std_gap(seq[i])) {
2703	if (keep_bases) --keep_bases;
2704	else seq[i] = '.';
2705	}
2706	}
2707
2708	GB_topSecurityLevel unsecured(gb_seq);
2709	TEST_EXPECT_NO_ERROR(GB_write_string(gb_seq, seq));
2710	free(seq);
2711	}
2712
2713	// remove CloTyrob
2714	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULp, PARSIMONY_ORG-1, env, false));
2715	TEST_EXPECT_COMBINES_PERFORMED(env, 4);
2716	TEST_EXPECT_EQUAL(env.root_node()->count_leafs(), 14);
2717
2718	// attempt to add CloTyrob (should fail because sequence too short) and CorGluta (should stay, because already in tree)
2719	TEST_REJECT_NULL(env.root_node()->findLeafNamed("CorGluta")); // has to be in tree
2720	{
2721	GB_transaction ta(env.gbmain());
2722	GBT_restore_marked_species(env.gbmain(), "CloTyrob;CorGluta");
2723	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2724	}
2725
2726	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_ADD, NULp, PARSIMONY_ORG-1, env, false));
2727	TEST_EXPECT_COMBINES_PERFORMED(env, 110); // @@@ why does this perform combines at all?
2728	TEST_EXPECT_EQUAL(env.root_node()->count_leafs(), 14); // ok, CloTyrob was not added
2729	TEST_REJECT_NULL(env.root_node()->findLeafNamed("CorGluta")); // has to be in tree
2730
2731	env.pop();
2732	}
2733	}
2734
2735	void TEST_prot_tree_modifications() {
2736	const char *aliname = "ali_tuf_pro";
2737
2738	PARSIMONY_testenv<AP_sequence_protein> env("TEST_prot.arb", aliname);
2739	TEST_EXPECT_NO_ERROR(env.load_tree("tree_prot_opti"));
2740	TEST_EXPECT_SAVED_TOPOLOGY(env, "prot-initial");
2741
2742	const int PARSIMONY_ORG = 1081;
2743	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 10);
2744
2745	// [PROTOPTI] opposed to nucleid tests above the initial tree here is already optimized! compare .@NUCOPTI
2746	// -> adding species approximately reproduces initial topology
2747	//
2748	// diff initial->add-quick: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/prot-add-quick.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fprot-initial.tree.expected
2749	// diff initial->add-NNI: http://bugs.arb-home.de/changeset/HEAD/branches/pars/UNIT_TESTER/run/pars/prot-add-NNI.tree.expected?old=HEAD&old_path=branches%2Fpars%2FUNIT_TESTER%2Frun%2Fpars%2Fprot-initial.tree.expected
2750	//
2751	// Note: comparing these two diffs also demonstrates why quick-adding w/o NNI suffers
2752
2753	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, "prot-removed", PARSIMONY_ORG-146, env, true)); // test remove-marked only (same code as part of nt_reAdd)
2754	TEST_EXPECT_COMBINES_PERFORMED(env, 5);
2755
2756	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "prot-add-quick", PARSIMONY_ORG, env, true)); // test quick-add
2757	TEST_EXPECT_COMBINES_PERFORMED(env, 213);
2758
2759	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_READD_NNI, "prot-add-NNI", PARSIMONY_ORG, env, true)); // test add + NNI
2760	TEST_EXPECT_COMBINES_PERFORMED(env, 262);
2761
2762	// test partial-add
2763	{
2764	GBDATA *gb_main = env.gbmain();
2765
2766	// create 2 non-overlapping partial species
2767	GBDATA *MucRaceP = createPartialSeqFrom(gb_main, aliname, "MucRaceP", "MucRacem", 0, 60+4); // (+4 = dots inserted into DB at left end)
2768	GBDATA *StrCoelP = createPartialSeqFrom(gb_main, aliname, "StrCoelP", "StrCoel9", 66-1+4, 184-1+4);
2769	GBDATA *StrCoelM = createPartialSeqFrom(gb_main, aliname, "StrCoelM", "StrCoel9", 66-1+4, 184-1+4);
2770	TEST_EXPECT_NO_ERROR(modifyOneBase(StrCoelM, aliname, 'Y', 'H')); // change first 'Y' into 'H'
2771
2772	// add StrCoelP (and undo)
2773	{
2774	env.push();
2775	// StrCoel9 and StrRamo3 do not differ in seq-range of partial -> any of both may be chosen as brother.
2776	// behavior should be changed with #605
2777	TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", "prot-addPart-StrCoelP", PARSIMONY_ORG, env));
2778	TEST_EXPECT_COMBINES_PERFORMED(env, 4);
2779	env.pop();
2780	}
2781
2782	{
2783	env.push();
2784	TEST_EXPECTATION(addingPartialResultsIn(MucRaceP, "MucRacem", "prot-addPart-MucRaceP", PARSIMONY_ORG, env)); // add MucRaceP
2785	TEST_EXPECT_COMBINES_PERFORMED(env, 6);
2786	TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", "prot-addPart-MucRaceP-StrCoelP", PARSIMONY_ORG, env)); // also add StrCoelP
2787	TEST_EXPECT_COMBINES_PERFORMED(env, 4);
2788	env.pop();
2789	}
2790
2791	// now add MucRaceP as full, then StrCoelP and StrCoelM as partials
2792	{
2793	env.push();
2794
2795	mark_only(MucRaceP);
2796	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2797	{
2798	GB_transaction ta(gb_main);
2799	TEST_EXPECT_NO_ERROR(GBT_write_int(MucRaceP, "ARB_partial", 0)); // revert species to "full"
2800	}
2801
2802	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "prot-addPartialAsFull-MucRaceP", PARSIMONY_ORG, env, false));
2803	TEST_EXPECT_COMBINES_PERFORMED(env, 156);
2804	TEST_EXPECT_EQUAL(is_partial(MucRaceP), 0); // check MucRaceP was added as full sequence
2805	TEST_EXPECTATION(addedAsBrotherOf("MucRaceP", "Eukarya EF-Tu", env)); // partial created from MucRacem gets inserted next to this group
2806	// Note: looks ok. group contains MucRacem, AbdGlauc and 4 other species
2807
2808	// add StrCoelP as partial.
2809	// as expected it is placed next to matching full sequences (does not differ in partial range)
2810	TEST_EXPECTATION(addingPartialResultsIn(StrCoelP, "StrCoel9;StrRamo3", NULp, PARSIMONY_ORG, env));
2811	TEST_EXPECT_COMBINES_PERFORMED(env, 4);
2812
2813	// StrCoelM differs slightly in overlap with StrCoel9/StrRamo3, but has no overlap with MucRaceP
2814	// shows bug described in #609 is fixed:
2815	TEST_EXPECTATION(addingPartialResultsIn(StrCoelM, "StrCoelP", "prot-addPart-bug609",
2816	PARSIMONY_ORG+1, // @@@ known bug - partial should not affect parsimony value; possibly related to ../HELP_SOURCE/oldhelp/pa_partial.hlp@WARNINGS
2817	env));
2818	TEST_EXPECT_COMBINES_PERFORMED(env, 5);
2819	env.pop();
2820	}
2821	}
2822
2823	TEST_EXPECT_SAVED_TOPOLOGY(env, "prot-initial");
2824
2825	const unsigned mixseed = 8164724;
2826
2827	const long PARSIMONY_MIXED = PARSIMONY_ORG + 1519;
2828	const long PARSIMONY_NNI_MARKED = PARSIMONY_ORG + 1053;
2829	const long PARSIMONY_NNI_ALL = PARSIMONY_ORG;
2830	const long PARSIMONY_OPTI_MARKED = PARSIMONY_ORG;
2831	const long PARSIMONY_OPTI_ALL = PARSIMONY_ORG; // no gain (initial tree already is optimized)
2832
2833	// ------------------------------------------------------
2834	// mix tree (original tree already is optimized)
2835
2836	GB_random_seed(mixseed);
2837	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_MIX_TREE, "prot-mixed", PARSIMONY_MIXED, env, false));
2838	TEST_EXPECT_COMBINES_PERFORMED(env, 89);
2839
2840	{
2841	env.push();
2842	TEST_EXPECTATION(movingRootDoesntAffectCosts(PARSIMONY_MIXED));
2843	TEST_EXPECT_COMBINES_PERFORMED(env, 234);
2844	env.pop();
2845	}
2846
2847	// ------------------------------
2848	// test optimize (some)
2849
2850	// mark initially marked species
2851	{
2852	GB_transaction ta(env.gbmain());
2853
2854	GBT_restore_marked_species(env.gbmain(), "CytLyti6;StrRamo3;MucRace2;SacCere5");
2855	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2856	}
2857
2858	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_MIXED);
2859
2860	// test branchlength calculation
2861	// (optimizations below implicitely recalculates branchlengths)
2862	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "prot-calclength", PARSIMONY_MIXED, env, false));
2863	TEST_EXPECT_COMBINES_PERFORMED(env, 79);
2864
2865	// test whether branchlength calculation depends on root-position
2866	{
2867	AP_tree_edge *orgRootEdge = rootEdge();
2868
2869	env.push();
2870
2871	const char *tested_roots[] = {
2872	// "CytLyti6", // no effect on branchlengths
2873	"TaxOcell",
2874	"MucRace3",
2875	"StrCoel9",
2876	};
2877
2878	for (size_t r = 0; r<ARRAY_ELEMS(tested_roots); ++r) {
2879	TEST_ANNOTATE(tested_roots[r]);
2880	const char *leafName = tested_roots[r];
2881	env.root_node()->findLeafNamed(leafName)->set_root();
2882	calc_branchlengths_and_reorder(env.graphic_tree());
2883	orgRootEdge->set_root();
2884	env.graphic_tree()->reorderTree(BIG_BRANCHES_TO_TOP);
2885
2886	TEST_EXPECT_SAVED_TOPOLOGY(env, "prot-calclength");
2887	}
2888	TEST_EXPECT_COMBINES_PERFORMED(env, 265);
2889
2890	env.pop();
2891	}
2892
2893	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "prot-opti-NNI", PARSIMONY_NNI_MARKED, env, true)); // test recursive NNI
2894	TEST_EXPECT_COMBINES_PERFORMED(env, 246);
2895
2896	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "prot-opti-marked-global", PARSIMONY_OPTI_MARKED, env, true)); // test recursive NNI+KL
2897	TEST_EXPECT_COMBINES_PERFORMED(env, 2810);
2898
2899	// -----------------------------
2900	// test optimize (all)
2901
2902	// mark all species
2903	mark_all(env.gbmain());
2904	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
2905	TEST_EXPECT_EQUAL(GBT_count_marked_species(env.gbmain()), 14);
2906
2907	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_MIXED);
2908
2909	// test branchlength calculation
2910	// (optimizations below implicitely recalculates branchlengths)
2911	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_CALC_LENS, "prot-calclength", PARSIMONY_MIXED, env, false));
2912	TEST_EXPECT_COMBINES_PERFORMED(env, 79);
2913
2914	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_NNI, "prot-opti-all-NNI", PARSIMONY_NNI_ALL, env, true)); // test recursive NNI
2915	TEST_EXPECT_COMBINES_PERFORMED(env, 359);
2916
2917	{
2918	env.push();
2919	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_OPTI_GLOBAL, "prot-opti-global", PARSIMONY_OPTI_ALL, env, false)); // test recursive NNI+KL
2920	TEST_EXPECT_COMBINES_PERFORMED(env, 1690);
2921
2922	TEST_EXPECTATION(movingRootDoesntAffectCosts(PARSIMONY_OPTI_ALL));
2923	TEST_EXPECT_COMBINES_PERFORMED(env, 215);
2924	env.pop();
2925	}
2926	}
2927
2928	void TEST_node_stack() {
2929	// test was used to fix #620
2930
2931	const char *aliname = "ali_5s";
2932	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
2933	TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
2934	TEST_EXPECT_SAVED_TOPOLOGY(env, "nucl-initial");
2935
2936	const int PARSIMONY_ORG = 302;
2937	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 14);
2938
2939	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2940	TEST_EXPECT(env.root_node()->has_valid_root_remarks());
2941
2942	// test set root to CytAquat + pop (works)
2943	{
2944	env.push();
2945	env.root_node()->findLeafNamed("CytAquat")->set_root();
2946	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2947	env.pop();
2948	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2949	}
2950
2951	TEST_EXPECT(env.root_node()->has_valid_root_remarks());
2952
2953	// test set root to CloButyr + pop (works)
2954	{
2955	env.push();
2956	env.root_node()->findLeafNamed("CloButyr")->set_root();
2957	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2958	env.pop();
2959	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2960	}
2961
2962	TEST_EXPECT(env.root_node()->has_valid_root_remarks());
2963
2964	// test set root to CloBifer + set root to CloTyrob + pop (works)
2965	// Note: both species are in same subtree (of root)
2966	{
2967	env.push();
2968
2969	env.root_node()->findLeafNamed("CloBifer")->set_root();
2970	env.root_node()->findLeafNamed("CloTyrob")->set_root();
2971
2972	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2973	env.pop();
2974	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2975	}
2976
2977	TEST_EXPECT(env.root_node()->has_valid_root_remarks());
2978
2979	// test set root to CytAquat + set root to CloButyr + pop (failed, fixed by [13138])
2980	TEST_EXPECT_COMBINES_PERFORMED(env, 0);
2981	for (int calcCostsBetween = 0; calcCostsBetween<2; ++calcCostsBetween) {
2982	TEST_ANNOTATE(GBS_global_string("calcCostsBetween=%i", calcCostsBetween));
2983
2984	TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
2985
2986	env.push();
2987
2988	env.root_node()->findLeafNamed("CytAquat")->set_root();
2989
2990	if (calcCostsBetween) {
2991	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 2);
2992	}
2993
2994	env.root_node()->findLeafNamed("CloButyr")->set_root();
2995
2996	TEST_VALIDITY(env.root_node()->sequence_state_valid());
2997	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 6);
2998
2999	env.pop();
3000
3001	TEST_VALIDITY(env.root_node()->sequence_state_valid());
3002	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_ORG);
3003	TEST_EXPECT_VALID_TREE(env.root_node());
3004	}
3005
3006	{
3007	env.push();
3008	{
3009	env.push();
3010
3011	env.root_node()->findLeafNamed("CloInnoc")->moveNextTo(env.root_node()->findLeafNamed("CytAquat"), 0.5);
3012	TEST_EXPECT_VALID_TREE(env.root_node());
3013	env.root_node()->findLeafNamed("CloInnoc")->set_root();
3014	TEST_EXPECT_VALID_TREE(env.root_node());
3015	env.root_node()->findLeafNamed("CytAquat")->moveNextTo(env.root_node()->findLeafNamed("CloPaste"), 0.5);
3016	TEST_EXPECT_VALID_TREE(env.root_node());
3017	env.root_node()->findLeafNamed("CloPaste")->set_root();
3018	TEST_EXPECT_VALID_TREE(env.root_node());
3019	env.root_node()->findLeafNamed("CloPaste")->moveNextTo(env.root_node()->findLeafNamed("CloInnoc"), 0.5);
3020	TEST_EXPECT_VALID_TREE(env.root_node());
3021
3022	{
3023	AP_tree_nlen *son_of_brother;
3024	AP_tree_nlen *brother_of_father;
3025
3026	// COVER1: son of root -> grandson of root
3027	{
3028	AP_tree_nlen *son_of_root = env.root_node()->get_leftson();
3029	ap_assert(son_of_root);
3030
3031	son_of_brother = son_of_root->get_brother()->get_leftson();
3032	son_of_root->moveNextTo(son_of_brother, 0.5);
3033	TEST_EXPECT_VALID_TREE(env.root_node());
3034	}
3035
3036	// COVER2: grandson of root -> son of brother
3037	{
3038	AP_tree_nlen *son_of_root = env.root_node()->get_leftson();
3039	AP_tree_nlen *grandson_of_root = son_of_root->get_brother()->get_rightson();
3040	ap_assert(grandson_of_root);
3041
3042	son_of_brother = grandson_of_root->get_brother()->get_leftson();
3043	grandson_of_root->moveNextTo(son_of_brother, 0.5);
3044	TEST_EXPECT_VALID_TREE(env.root_node());
3045	}
3046
3047	AP_tree_nlen *some_leaf = env.root_node()->findLeafNamed("CloBifer");
3048	ap_assert(some_leaf);
3049
3050	// COVER3: some leaf -> son of brother
3051	son_of_brother = some_leaf->get_brother()->get_leftson();
3052	some_leaf->moveNextTo(son_of_brother, 0.5);
3053	TEST_EXPECT_VALID_TREE(env.root_node());
3054
3055	// COVER4: some leaf -> son of brother
3056	brother_of_father = some_leaf->get_father()->get_brother();
3057	some_leaf->moveNextTo(brother_of_father, 0.5);
3058	TEST_EXPECT_VALID_TREE(env.root_node());
3059
3060	// test forbidden moves:
3061	TEST_EXPECT_ERROR_CONTAINS(some_leaf->cantMoveNextTo(some_leaf->get_father()), "Already there");
3062	TEST_EXPECT_ERROR_CONTAINS(some_leaf->cantMoveNextTo(some_leaf->get_brother()), "Already there");
3063	}
3064
3065	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG+5, 6);
3066
3067	env.pop();
3068	TEST_EXPECT_VALID_TREE(env.root_node());
3069	}
3070	env.pop();
3071
3072	TEST_EXPECT_KNOWN_PARSVAL(env, PARSIMONY_ORG);
3073	TEST_EXPECT_VALID_TREE(env.root_node());
3074	}
3075
3076	// remove + quick add marked + pop() both works
3077	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, "nucl-add-quick", PARSIMONY_ORG-18, env, true)); // test quick-add
3078	TEST_EXPECT_COMBINES_PERFORMED(env, 400);
3079
3080	TEST_EXPECT(env.root_node()->has_valid_root_remarks());
3081
3082	// remove + quick-add marked + pop() quick-add -> corrupts tree
3083	// (root-edge is lost)
3084	{
3085	env.push();
3086	TEST_EXPECT_VALID_TREE(env.root_node());
3087	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULp, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
3088	TEST_EXPECT_COMBINES_PERFORMED(env, 0);
3089	TEST_EXPECT_VALID_TREE(env.root_node());
3090
3091	TEST_VALIDITY(env.all_available_pops_will_produce_valid_trees());
3092
3093	env.push();
3094	TEST_EXPECT_VALID_TREE(env.root_node());
3095	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, NULp, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
3096	TEST_EXPECT_COMBINES_PERFORMED(env, 400);
3097	TEST_EXPECT_VALID_TREE(env.root_node());
3098	TEST_VALIDITY(env.all_available_pops_will_produce_valid_trees());
3099	env.pop();
3100
3101	TEST_VALIDITY(env.root_node()->has_valid_edges());
3102
3103	env.pop();
3104	TEST_EXPECT_VALID_TREE(env.root_node());
3105	TEST_EXPECT_COMBINES_PERFORMED(env, 0);
3106	}
3107
3108	// same as above, but with only 1 species marked
3109	const char *testSingle[] = {
3110	"CytAquat", // CytAquat is the only grandson of root (CytAquat located in lower subtree)
3111	"CloBifer", // two father nodes between CloBifer and root (CloBifer located in upper subtree)
3112	"CloPaste", // two father nodes between CloPaste and root (CloPaste located in upper subtree)
3113	"CorGluta", // three father nodes between CorGluta and root (CorGluta located in lower subtree)
3114	"CelBiazo", // two father nodes between CelBiazo and root
3115	NULp
3116	};
3117
3118	for (int i = 0; testSingle[i]; ++i) {
3119	for (int swapped = 0; swapped<2; ++swapped) {
3120	TEST_ANNOTATE(GBS_global_string("single=%s swapped=%i", testSingle[i], swapped));
3121
3122	env.push();
3123	TEST_EXPECT_VALID_TREE(env.root_node());
3124	{
3125	AP_tree_nlen *old_rightson = env.root_node()->get_rightson();
3126	env.root_node()->get_leftson()->get_rightson()->set_root();
3127	old_rightson->get_leftson()->set_root();
3128	old_rightson->set_root();
3129
3130	ap_assert(env.root_node()->get_rightson() == old_rightson);
3131	}
3132	TEST_EXPECT_VALID_TREE(env.root_node());
3133
3134	mark_only(env.root_node()->findLeafNamed(testSingle[i])->gb_node);
3135	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
3136
3137	env.push();
3138	if (swapped) {
3139	env.root_node()->swap_sons();
3140	}
3141
3142	TEST_EXPECT_VALID_TREE(env.root_node());
3143	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULp, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
3144	TEST_EXPECT_VALID_TREE(env.root_node());
3145
3146	env.push();
3147	TEST_EXPECT_VALID_TREE(env.root_node());
3148	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, NULp, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
3149	TEST_EXPECT_VALID_TREE(env.root_node());
3150	env.pop();
3151
3152	TEST_EXPECT_VALID_TREE(env.root_node());
3153	env.pop();
3154	TEST_EXPECT_VALID_TREE(env.root_node());
3155	env.pop();
3156	TEST_EXPECT_VALID_TREE(env.root_node());
3157	}
3158	}
3159	TEST_EXPECT_COMBINES_PERFORMED(env, 2120);
3160
3161	// similar to above (remove+add a grandson of root; here grandson is a subtree with 4 species)
3162
3163	for (int remove_from_lower_subtree = 0; remove_from_lower_subtree<2; ++remove_from_lower_subtree) {
3164	TEST_ANNOTATE(GBS_global_string("remove_from_lower_subtree=%i", remove_from_lower_subtree));
3165
3166	// mark a complete subtree (which - as a whole - forms a grandson of root). subtree is located in upper part of the tree
3167	mark_only(env.root_node()->findLeafNamed("CloButy2")->gb_node);
3168	mark(env.root_node()->findLeafNamed("CloButyr")->gb_node);
3169	mark(env.root_node()->findLeafNamed("CloCarni")->gb_node);
3170	mark(env.root_node()->findLeafNamed("CloPaste")->gb_node);
3171	env.compute_tree(); // species marks affect order of node-chain (used in nni_rec)
3172
3173	env.push();
3174	if (remove_from_lower_subtree) {
3175	env.root_node()->swap_sons();
3176	}
3177	TEST_EXPECT_VALID_TREE(env.root_node());
3178	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_REMOVE_MARKED, NULp, -1, env, false)); // test remove-marked only (same code as part of nt_reAdd)
3179	TEST_EXPECT_VALID_TREE(env.root_node());
3180
3181	env.push();
3182	TEST_EXPECT_VALID_TREE(env.root_node());
3183	TEST_EXPECTATION(modifyingTopoResultsIn(MOD_QUICK_READD, NULp, -1, env, false)); // test quick-add (same code as part of nt_reAdd)
3184	TEST_EXPECT_VALID_TREE(env.root_node());
3185	env.pop();
3186
3187	TEST_VALIDITY(env.root_node()->has_valid_edges()); // now always valid
3188
3189	env.pop();
3190	TEST_EXPECT_VALID_TREE(env.root_node());
3191	}
3192	TEST_EXPECT_COMBINES_PERFORMED(env, 584);
3193	}
3194
3195	void TEST_node_edge_resources() {
3196	const char *aliname = "ali_5s";
3197
3198	// document memsize of nodes and edges:
3199
3200	#define STATE_STACK_SIZE sizeof(StateStack) // 8 (Cxx11) or 16 (older C++); maybe 4/8 in 32bit
3201
3202	#if defined(ARB_64)
3203	TEST_EXPECT_EQUAL(sizeof(AP_tree_nlen), 184 + STATE_STACK_SIZE);
3204	TEST_EXPECT_EQUAL(sizeof(AP_tree), 136);
3205	TEST_EXPECT_EQUAL(sizeof(ARB_seqtree), 88);
3206	TEST_EXPECT_EQUAL(sizeof(TreeNode), 80);
3207	#else // !defined(ARB_64)
3208	TEST_EXPECT_EQUAL(sizeof(AP_tree_nlen), 112 + STATE_STACK_SIZE);
3209	TEST_EXPECT_EQUAL(sizeof(AP_tree), 84);
3210	TEST_EXPECT_EQUAL(sizeof(ARB_seqtree), 48);
3211	TEST_EXPECT_EQUAL(sizeof(TreeNode), 44);
3212	#endif
3213
3214
3215	#if defined(ARB_64)
3216	TEST_EXPECT_EQUAL(sizeof(AP_tree_edge), 64);
3217	#else // !defined(ARB_64)
3218	TEST_EXPECT_EQUAL(sizeof(AP_tree_edge), 32);
3219	#endif
3220
3221	PARSIMONY_testenv<AP_sequence_parsimony> env("TEST_trees.arb", aliname);
3222	TEST_EXPECT_NO_ERROR(env.load_tree("tree_test"));
3223
3224	const int PARSIMONY_ORG = 302;
3225	TEST_EXPECT_ONLY_PARSVAL_COMBINES(env, PARSIMONY_ORG, 14);
3226
3227	AP_tree_nlen *CloButyr = env.root_node()->findLeafNamed("CloButyr");
3228	AP_tree_nlen *CloButy2 = env.root_node()->findLeafNamed("CloButy2");
3229	TEST_EXPECT_EQUAL(CloButyr->get_brother()->name, CloButy2->name); // test they are brothers
3230
3231	AP_tree_nlen *CorAquat = env.root_node()->findLeafNamed("CorAquat");
3232	AP_tree_nlen *CurCitre = env.root_node()->findLeafNamed("CurCitre");
3233	TEST_EXPECT_EQUAL(CorAquat->get_brother()->name, CurCitre->name); // test they are brothers
3234
3235	CorAquat->REMOVE();
3236
3237	for (int test = 1; test<=8; ++test) {
3238	// test == 1 -> provokes common nodes+edges in revert+accept
3239	// test == 2 -> provokes common nodes+edges in revert+accept
3240	// test == 3 -> provokes common nodes+edges in revert+accept
3241	// tests 4-7 do not provoke common nodes or edges
3242
3243	for (int mode = 0; mode<=3; ++mode) {
3244	bool accept_outer = mode&2;
3245	bool accept_inner = mode&1;
3246
3247	TEST_ANNOTATE(GBS_global_string("accept_outer=%i accept_inner=%i (mode=%i, test=%i)", accept_outer, accept_inner, mode, test));
3248
3249	TEST_EXPECT_NULL(CorAquat->get_father());
3250	TEST_EXPECT(CloButyr->get_brother() == CloButy2);
3251	TEST_EXPECT_VALID_TREE(env.root_node());
3252
3253	env.push();
3254
3255	switch (test) {
3256	case 1: CorAquat->insert(CurCitre); break;
3257	case 2: CorAquat->insert(CurCitre); break;
3258	case 3: break;
3259	case 4: CloButyr->REMOVE(); break;
3260	case 5: CloButyr->REMOVE(); break;
3261	case 6: break;
3262	case 7: CloButyr->moveNextTo(CurCitre, 0.5); break;
3263	case 8: break;
3264	default: ap_assert(0); break;
3265	}
3266	TEST_EXPECT_VALID_TREE(env.root_node());
3267
3268	{
3269	env.push();
3270
3271	switch (test) {
3272	case 1: CorAquat->REMOVE(); break;
3273	case 2: break;
3274	case 3: CorAquat->insert(CurCitre); break;
3275	case 4: CloButyr->insert(CloButy2); break;
3276	case 5: break;
3277	case 6: CloButyr->REMOVE(); break;
3278	case 7: CloButyr->moveNextTo(CloButy2, 0.5); break;
3279	case 8: CorAquat->insert(CurCitre); CorAquat->REMOVE(); break;
3280	default: ap_assert(0); break;
3281	}
3282	TEST_EXPECT_VALID_TREE(env.root_node());
3283
3284	env.accept_if(accept_inner);
3285	}
3286
3287	switch (test) {
3288	case 1: break;
3289	case 2: CorAquat->REMOVE(); break;
3290	case 3: if (CorAquat->father) CorAquat->REMOVE(); break;
3291	case 4: break;
3292	case 5: CloButyr->insert(CloButy2); break;
3293	case 6: if (!CloButyr->father) CloButyr->insert(CloButy2); break;
3294	case 7: CloButyr->REMOVE(); break;
3295	case 8: break;
3296	default: ap_assert(0); break;
3297	}
3298	TEST_EXPECT_VALID_TREE(env.root_node());
3299
3300	env.accept_if(accept_outer);
3301
3302	// manually revert changes (outside any stack frame)
3303	if (CorAquat->father) CorAquat->REMOVE();
3304	if (!CloButyr->father) CloButyr->insert(CloButy2);
3305	}
3306	}
3307
3308	CorAquat->insert(CurCitre);
3309
3310	TEST_EXPECT_PARSVAL(env, PARSIMONY_ORG);
3311	env.combines_performed(); // accept any no of combines
3312	}
3313
3314	#endif // UNIT_TESTS
3315
3316	// --------------------------------------------------------------------------------
3317

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

Download in other formats: