source: branches/nameserver/SL/AP_TREE/AP_Tree.cxx

Last change on this file was 18125, checked in by westram, 5 years ago
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 52.1 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : AP_Tree.cxx                                       //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "AP_Tree.hxx"
12#include "AP_TreeShader.hxx"
13
14#include <AP_filter.hxx>
15#include <aw_msg.hxx>
16#include <arb_progress.h>
17#include <arb_str.h>
18#include <ad_cb.h>
19
20#include <math.h>
21#include <map>
22#include <climits>
23#include <algorithm>
24
25using namespace std;
26
27/*!***************************************************************************************
28 ************           Rates                       **********
29 *****************************************************************************************/
30void AP_rates::print() {
31    AP_FLOAT max;
32    int i;
33
34    max = 0.0;
35    for (i=0; i<rate_len; i++) {
36        if (rates[i] > max) max = rates[i];
37    }
38    printf("rates:");
39    for (i=0; i<rate_len; i++) {
40        putchar('0' + (int)(rates[i]/max*9.9));
41    }
42    printf("\n");
43}
44
45AP_rates::AP_rates() {
46    memset ((char *)this, 0, sizeof(AP_rates));
47}
48
49char *AP_rates::init(AP_filter *fil) {
50    if (fil->get_timestamp() <= this->update) return NULp;
51
52    rate_len = fil->get_filtered_length();
53    delete [] rates;
54    rates = new AP_FLOAT[rate_len];
55    for (int i=0; i<rate_len; i++) { // LOOP_VECTORIZED // tested down to gcc 5.5.0 (may fail on older gcc versions)
56        rates[i] = 1.0;
57    }
58    this->update = fil->get_timestamp();
59    return NULp;
60}
61
62char *AP_rates::init(AP_FLOAT * ra, AP_filter *fil) {
63    if (fil->get_timestamp() <= this->update) return NULp;
64
65    rate_len = fil->get_filtered_length();
66    delete [] rates;
67    rates = new AP_FLOAT[rate_len];
68    int i, j;
69    for (j=i=0; i<rate_len; j++) {
70        if (fil->use_position(j)) {
71            rates[i++] = ra[j];
72        }
73    }
74    this->update = fil->get_timestamp();
75    return NULp;
76}
77
78AP_rates::~AP_rates() { delete [] rates; }
79
80
81/*!***************************************************************************************
82************           AP_tree_root                    **********
83*****************************************************************************************/
84
85AP_tree_root::AP_tree_root(AliView *aliView, AP_sequence *seq_proto, bool add_delete_callbacks, const group_scaling *scaling)
86    : ARB_seqtree_root(aliView, seq_proto, add_delete_callbacks),
87      root_changed_cb(NULp), root_changed_cd(NULp),
88      node_deleted_cb(NULp), node_deleted_cd(NULp),
89      gScale(scaling),
90      gb_tree_gone(NULp),
91      gone_tree_name(NULp),
92      tree_timer(0),
93      species_timer(0),
94      rates(NULp)
95{
96    GBDATA         *gb_main = get_gb_main();
97    GB_transaction  ta(gb_main);
98
99    gb_species_data = GBT_get_species_data(gb_main);
100}
101
102AP_tree_root::~AP_tree_root() {
103    predelete();
104    free(gone_tree_name);
105    ap_assert(!get_root_node());
106}
107
108bool AP_tree_root::is_tree_updated() {
109    GBDATA *gbtree = get_gb_tree();
110    if (gbtree) {
111        GB_transaction ta(gbtree);
112        return GB_read_clock(gbtree)>tree_timer;
113    }
114    return true;
115}
116
117bool AP_tree_root::is_species_updated() {
118    if (gb_species_data) {
119        GB_transaction ta(gb_species_data);
120        return GB_read_clock(gb_species_data)>species_timer;
121    }
122    return true;
123}
124
125void AP_tree_root::update_timers() {
126    if (gb_species_data) {
127        GB_transaction  ta(GB_get_root(gb_species_data));
128        GBDATA         *gbtree = get_gb_tree();
129        if (gbtree) tree_timer = GB_read_clock(gbtree);
130        species_timer          = GB_read_clock(gb_species_data);
131    }
132}
133
134/*!***************************************************************************************
135************           AP_tree                     **********
136*****************************************************************************************/
137
138static void ap_tree_node_deleted(GBDATA *, AP_tree *tree) {
139    tree->gb_node = NULp;
140}
141
142AP_tree::~AP_tree() {
143    if (gr.callback_exists && gb_node) {
144        GB_remove_callback(gb_node, GB_CB_DELETE, makeDatabaseCallback(ap_tree_node_deleted, this));
145    }
146
147    AP_tree_root *root = get_tree_root();
148    if (root) root->inform_about_delete(this);
149}
150
151void AP_tree::initial_insert(AP_tree *new_brother, AP_tree_root *troot) {
152    ap_assert(troot);
153    ap_assert(is_leaf());
154    ap_assert(new_brother->is_leaf());
155    ap_assert(!troot->get_root_node());
156
157    ASSERT_VALID_TREE(this);
158    ASSERT_VALID_TREE(new_brother);
159
160    AP_tree *new_root = DOWNCAST(AP_tree*, troot->makeNode());
161
162    new_root->leftson  = this;
163    new_root->rightson = new_brother;
164    new_root->father   = NULp;
165
166    father              = new_root;
167    new_brother->father = new_root;
168
169    new_root->leftlen  = 0.5;
170    new_root->rightlen = 0.5;
171
172    troot->change_root(NULp, new_root);
173
174    set_tree_root(troot);
175    new_brother->set_tree_root(troot);
176}
177
178void AP_tree::insert(AP_tree *new_brother) {
179    ASSERT_VALID_TREE(this);
180    ASSERT_VALID_TREE(new_brother);
181    ap_assert(new_brother->get_tree_root()->get_root_node()->has_valid_root_remarks());
182
183    AP_tree  *new_tree = DOWNCAST(AP_tree*, new_brother->get_tree_root()->makeNode());
184    AP_FLOAT  laenge;
185
186    if (new_brother->is_son_of_root()) {
187        new_brother->get_father()->remove_root_remark();
188    }
189
190    new_tree->leftson  = this;
191    new_tree->rightson = new_brother;
192    new_tree->father   = new_brother->father;
193    father             = new_tree;
194
195    if (new_brother->father) {
196        if (new_brother->father->leftson == new_brother) {
197            laenge = new_brother->father->leftlen = new_brother->father->leftlen*.5;
198            new_brother->father->leftson = new_tree;
199        }
200        else {
201            laenge = new_brother->father->rightlen = new_brother->father->rightlen*.5;
202            new_brother->father->rightson = new_tree;
203        }
204    }
205    else {
206        laenge = 0.5;
207    }
208    new_tree->leftlen   = laenge;
209    new_tree->rightlen  = laenge;
210    new_brother->father = new_tree;
211
212    AP_tree_root *troot = new_brother->get_tree_root();
213    ap_assert(troot); // Note: initial_insert() has to be used to build initial tree
214
215    if (!new_tree->father) troot->change_root(new_brother, new_tree);
216    new_tree->set_tree_root(troot);
217    set_tree_root(troot);
218
219    ASSERT_VALID_TREE(troot->get_root_node());
220    ap_assert(get_tree_root()->get_root_node()->has_valid_root_remarks());
221}
222
223void AP_tree_root::change_root(TreeNode *oldroot, TreeNode *newroot) {
224    if (root_changed_cb) { // @@@ better call after calling base::change_root?
225        root_changed_cb(root_changed_cd, DOWNCAST(AP_tree*, oldroot), DOWNCAST(AP_tree*, newroot));
226    }
227    if (!oldroot) {
228        ap_assert(newroot);
229        if (gb_tree_gone) {                         // when tree was temporarily deleted (e.g. by 'Remove & add all')
230            set_gb_tree(gb_tree_gone);              // re-use previous DB entry
231            gb_tree_gone = NULp;
232        }
233    }
234    if (!newroot) {                                 // tree empty
235        GBDATA *gbtree = get_gb_tree();
236        if (gbtree) {
237            ap_assert(!gb_tree_gone);               // no tree should be remembered yet
238            gb_tree_gone = gbtree;                  // remember for deletion (done in AP_tree::save)
239        }
240    }
241    ARB_seqtree_root::change_root(oldroot, newroot);
242}
243
244void AP_tree_root::inform_about_delete(AP_tree *del) {
245    if (node_deleted_cb) node_deleted_cb(node_deleted_cd, del);
246}
247
248void AP_tree_root::set_root_changed_callback(AP_rootChangedCb cb, void *cd) {
249    root_changed_cb = cb;
250    root_changed_cd = cd;
251}
252
253void AP_tree_root::set_node_deleted_callback(AP_nodeDelCb cb, void *cd) {
254    node_deleted_cb = cb;
255    node_deleted_cd = cd;
256}
257
258
259AP_tree *AP_tree::REMOVE() {
260    // Remove this + father from tree. Father node will be destroyed.
261    // Caller has to destroy the removed node (if intended).
262    //
263    // Warning: when removing the 2nd to last node from the tree,
264    // the whole tree will be removed.
265    // In that case both leaf nodes remain undestroyed.
266
267    ASSERT_VALID_TREE(this);
268    if (!father) {
269        get_tree_root()->change_root(this, NULp); // tell AP_tree_root that the root node has been removed
270        forget_origin(); // removed nodes are rootless
271    }
272    else {
273        AP_tree      *brother     = get_brother();  // brother remains in tree
274        GBT_LEN       brothersLen = brother->get_branchlength();
275        AP_tree      *fath        = get_father();   // fath of this is removed as well
276        ARB_seqtree  *grandfather = fath->get_father();
277        AP_tree_root *troot       = get_tree_root();
278
279        if (fath->gb_node) {                      // move inner information to remaining subtree
280            if (!brother->gb_node && !brother->is_leaf()) {
281                brother->gb_node = fath->gb_node;
282                fath->gb_node    = NULp;
283            }
284        }
285
286        remove_remarks_from_this_and_parent(); // remove remarks of this + father
287
288        if (grandfather) {
289            brother->unlink_from_father();
290
291            bool wasLeftSon = fath->is_leftson();
292            fath->unlink_from_father();
293
294            if (wasLeftSon) {
295                ap_assert(!grandfather->leftson);
296                grandfather->leftlen += brothersLen;
297                grandfather->leftson  = brother;
298            }
299            else {
300                ap_assert(!grandfather->rightson);
301                grandfather->rightlen += brothersLen;
302                grandfather->rightson  = brother;
303            }
304            brother->father = grandfather;
305            if (!grandfather->father) {
306                ap_assert(brother->is_son_of_root());
307                if (!brother->is_leaf()) brother->remove_remark();
308            }
309        }
310        else {                                      // father is root, make brother the new root
311            if (brother->is_leaf()) {
312                troot->change_root(fath, NULp);     // erase tree from root
313                brother->unlink_from_father();      // do not automatically delete brother
314            }
315            else {
316                brother->unlink_from_father();
317                troot->change_root(fath, brother);
318            }
319        }
320
321        ap_assert(fath == father);
322
323        ASSERT_VALID_TREE_OR_NULL(troot->get_root_node());
324
325        troot->inform_about_delete(fath);
326        troot->inform_about_delete(this);
327
328        fath->forget_origin();
329        ASSERT_VALID_TREE(fath);
330
331        unlink_from_father();
332        destroy(fath, troot);
333        ASSERT_VALID_TREE(this);
334    }
335    return this;
336}
337
338GB_ERROR AP_tree::cantMoveNextTo(AP_tree *new_brother) {
339    GB_ERROR error = NULp;
340
341    if (!father)                                error = "Can't move the root of the tree";
342    else if (!new_brother->father)              error = "Can't move to the root of the tree";
343    else if (new_brother->father == father)     error = "Already there";
344    else if (new_brother == father)             error = "Already there";
345    else if (!father->father)                   error = "Can't move son of root";
346    else if (new_brother->is_inside(this))      error = "Can't move a subtree into itself";
347
348    return error;
349}
350
351void AP_tree::moveNextTo(AP_tree *new_brother, AP_FLOAT rel_pos) {
352    // rel_pos: 0.0 -> branch at father; 1.0 -> branch at brother; 0.5 -> branch at half distance between father and brother
353
354    // @@@ "move subtree" needs better correction for groups (esp. if moving from root-of-group or when keeled groups are involved; see #785)
355
356    ap_assert(father);
357    ap_assert(new_brother);
358    ap_assert(new_brother->father);
359    ap_assert(new_brother->father != father);       // already there
360    ap_assert(new_brother != father);               // already there
361
362    ap_assert(!new_brother->is_inside(this));       // can't move tree into itself
363    ap_assert(get_tree_root()->get_root_node()->has_valid_root_remarks());
364
365    remove_remarks_from_this_and_parent();
366    get_brother()->smart_remove_remark();
367    new_brother->remove_remarks_from_this_and_parent();
368
369    if (father->leftson != this) get_father()->swap_sons();
370
371    AP_tree *new_root = NULp;
372    if (!father->father) { // move son of root
373        ap_assert(!father->has_group_info());
374
375        get_father()->remove_root_remark();
376
377        new_root         = get_brother();
378        new_root->father = NULp;
379
380        ap_assert(!new_root->is_leaf()); // a leaf is no valid tree (should be impossible, because new_brother!=brother)
381    }
382    else {
383        AP_tree *grandfather = get_father()->get_father();
384        if (!grandfather->father) { // move grandchild of root
385            grandfather->remove_root_remark();
386        }
387        if (grandfather->leftson == father) {
388            grandfather->leftlen += father->rightlen;
389            grandfather->leftson  = father->rightson;
390        }
391        else {
392            grandfather->rightlen += father->rightlen;
393            grandfather->rightson  = father->rightson;
394        }
395        father->rightson->father = grandfather;
396    }
397
398    AP_tree  *new_tree       = get_father();
399    AP_tree  *brother_father = new_brother->get_father();
400    AP_FLOAT  laenge;
401
402    if (brother_father->leftson == new_brother) {
403        laenge  = brother_father->leftlen;
404        laenge -= brother_father->leftlen = laenge * rel_pos;
405        brother_father->leftson = new_tree;
406    }
407    else {
408        laenge  = brother_father->rightlen;
409        laenge -= brother_father->rightlen = laenge * rel_pos;
410        brother_father->rightson = new_tree;
411    }
412
413    new_tree->rightlen  = laenge;
414    new_brother->father = new_tree;
415    new_tree->rightson  = new_brother;
416    new_tree->father    = brother_father;
417
418    if (new_root) {
419        new_tree->get_tree_root()->change_root(new_tree, new_root);
420        new_root->remove_root_remark();
421    }
422
423    ap_assert(new_tree->get_tree_root()->get_root_node()->has_valid_root_remarks());
424}
425
426inline int tree_read_byte(GBDATA *tree, const char *key, int init) {
427    if (tree) {
428        GBDATA *gbd = GB_entry(tree, key);
429        if (gbd) return GB_read_byte(gbd);
430    }
431    return init;
432}
433
434inline float tree_read_float(GBDATA *tree, const char *key, float init) {
435    if (tree) {
436        GBDATA *gbd = GB_entry(tree, key);
437        if (gbd) return GB_read_float(gbd);
438    }
439    return init;
440}
441
442
443
444//! moves all node/leaf information from struct TreeNode to AP_tree
445void AP_tree::load_node_info() {
446    gr.spread          = tree_read_float(gb_node, "spread",          1.0);
447    gr.left_angle      = tree_read_float(gb_node, "left_angle",      0.0);
448    gr.right_angle     = tree_read_float(gb_node, "right_angle",     0.0);
449    gr.left_linewidth  = tree_read_byte (gb_node, "left_linewidth",  0);
450    gr.right_linewidth = tree_read_byte (gb_node, "right_linewidth", 0);
451    gr.grouped         = tree_read_byte (gb_node, "grouped",         0);
452}
453
454void AP_tree::load_subtree_info() {
455    load_node_info();
456    if (!is_leaf()) {
457        get_leftson()->load_subtree_info();
458        get_rightson()->load_subtree_info();
459    }
460}
461
462#if defined(DEBUG)
463#if defined(DEVEL_RALF)
464#define DEBUG_tree_write_byte
465#endif // DEVEL_RALF
466#endif // DEBUG
467
468
469static GB_ERROR tree_write_byte(GBDATA *gb_tree, AP_tree *node, short i, const char *key, int init) {
470    GBDATA   *gbd;
471    GB_ERROR  error = NULp;
472    if (i==init) {
473        if (node->gb_node) {
474            gbd = GB_entry(node->gb_node, key);
475            if (gbd) {
476#if defined(DEBUG_tree_write_byte)
477                printf("[tree_write_byte] deleting db entry %p\n", gbd);
478#endif // DEBUG_tree_write_byte
479                GB_delete(gbd);
480            }
481        }
482    }
483    else {
484        if (!node->gb_node) {
485            node->gb_node = GB_create_container(gb_tree, "node");
486#if defined(DEBUG_tree_write_byte)
487            printf("[tree_write_byte] created node-container %p\n", node->gb_node);
488#endif // DEBUG_tree_write_byte
489        }
490        gbd = GB_entry(node->gb_node, key);
491        if (!gbd) {
492            gbd = GB_create(node->gb_node, key, GB_BYTE);
493#if defined(DEBUG_tree_write_byte)
494            printf("[tree_write_byte] created db entry %p\n", gbd);
495#endif // DEBUG_tree_write_byte
496        }
497        error = GB_write_byte(gbd, i);
498    }
499    return error;
500}
501
502static GB_ERROR tree_write_float(GBDATA *gb_tree, AP_tree *node, float f, const char *key, float init) {
503    GB_ERROR error = NULp;
504    if (f==init) {
505        if (node->gb_node) {
506            GBDATA *gbd    = GB_entry(node->gb_node, key);
507            if (gbd) error = GB_delete(gbd);
508        }
509    }
510    else {
511        if (!node->gb_node) {
512            node->gb_node             = GB_create_container(gb_tree, "node");
513            if (!node->gb_node) error = GB_await_error();
514        }
515        if (!error) error = GBT_write_float(node->gb_node, key, f);
516    }
517    return error;
518}
519
520GB_ERROR AP_tree::tree_write_tree_rek(GBDATA *gb_tree) {
521    GB_ERROR error = NULp;
522    if (!is_leaf()) {
523        error             = get_leftson()->tree_write_tree_rek(gb_tree);
524        if (!error) error = get_rightson()->tree_write_tree_rek(gb_tree);
525
526        if (!error) error = tree_write_float(gb_tree, this, gr.spread,          "spread",          1.0);
527        if (!error) error = tree_write_float(gb_tree, this, gr.left_angle,      "left_angle",      0.0);
528        if (!error) error = tree_write_float(gb_tree, this, gr.right_angle,     "right_angle",     0.0);
529        if (!error) error = tree_write_byte (gb_tree, this, gr.left_linewidth,  "left_linewidth",  0);
530        if (!error) error = tree_write_byte (gb_tree, this, gr.right_linewidth, "right_linewidth", 0);
531        if (!error) error = tree_write_byte (gb_tree, this, gr.grouped,         "grouped",         0);
532    }
533    return error;
534}
535
536GB_ERROR AP_tree_root::saveToDB() {
537    GB_ERROR  error  = GB_push_transaction(get_gb_main());
538    if (get_gb_tree()) {
539        error = get_root_node()->tree_write_tree_rek(get_gb_tree());
540    }
541    else {
542        ap_assert(!gb_tree_gone);      // should have been handled by caller (e.g. in AWT_graphic_tree::save)
543    }
544    if (!error) {
545        if (!get_gb_tree() && gone_tree_name) { // tree was deleted before
546            GBDATA *gb_tree_exists = GBT_find_tree(get_gb_main(), gone_tree_name);
547            if (gb_tree_exists) {
548                error = "tree already exists";
549            }
550            else {
551                error = GBT_write_tree(get_gb_main(), gone_tree_name, get_root_node());
552                if (!error) {
553                    gb_tree_exists = GBT_find_tree(get_gb_main(), gone_tree_name);
554                    ap_assert(gb_tree_exists);
555                    if (gb_tree_exists) {
556                        set_gb_tree_and_name(GBT_find_tree(get_gb_main(), gone_tree_name), gone_tree_name);
557                        aw_message(GBS_global_string("Recreated previously deleted '%s'", gone_tree_name));
558                        freenull(gone_tree_name);
559                    }
560                }
561            }
562
563            if (error) aw_message(GBS_global_string("Failed to recreate '%s' (Reason: %s)", gone_tree_name, error));
564        }
565
566        if (!error) error = ARB_seqtree_root::saveToDB();
567    }
568    if (!error) update_timers();
569
570    return GB_end_transaction(get_gb_main(), error);
571}
572
573inline GBDATA *find_group_name_entry(TreeNode *node) { return node->has_group_info() ? GB_entry(node->gb_node, "group_name") : NULp; }
574
575inline void TreeNode::swap_node_info(TreeNode *other, bool ofKeeledGroups) {
576    if (ofKeeledGroups) {
577        // member 'inverseLeft' cannot be handled correctly w/o knowing son nodes
578        get_father()->swap_node_info(other->get_father(), false);
579        fixKeeledOrientation();
580        other->fixKeeledOrientation();
581    }
582    else if (this == other) {
583        gb_assert(keeledOver && other->keeledOver);
584        inverseLeft = !inverseLeft;
585    }
586    else {
587        std::swap(name, other->name);
588        std::swap(gb_node, other->gb_node);
589        std::swap(keeledOver, other->keeledOver);
590    }
591}
592
593GB_ERROR AP_tree::swap_group_with(AP_tree *dest, bool swapKeeled) {
594    GB_ERROR error = NULp;
595    if (swapKeeled) {
596        ap_assert(!is_leaf() && is_keeled_group());
597
598        AP_tree *parent      = get_father();
599        AP_tree *dest_parent = dest->get_father();
600
601        if (parent != dest_parent && dest_parent->has_group_info() && dest_parent->keelTarget() != dest ) {
602            error = GBS_global_string("cannot move group '%s' (would create partial overlap with '%s')", parent->name, dest_parent->name);
603        }
604        if (!error && dest_parent->is_root_node()) {
605            error = "invalid move of keeled group to tree-root";
606        }
607        if (!error) {
608            swap_node_info(dest, true);
609            parent->load_node_info();
610            dest_parent->load_node_info();
611        }
612    }
613    else {
614        ap_assert(!is_leaf() && is_normal_group());
615        if (dest->has_group_info() && !dest->is_normal_group()) {
616            error = GBS_global_string("cannot move group '%s' (would create partial overlap with '%s')", name, dest->name);
617        }
618        if (!error) {
619            swap_node_info(dest, false);
620            load_node_info();
621            dest->load_node_info();
622        }
623    }
624    return error;
625}
626
627GB_ERROR AP_tree::move_group_to(AP_tree *dest) {
628    GB_ERROR error = NULp;
629
630    bool src_normal = !is_leaf() && is_normal_group();
631    bool src_keeled = !is_leaf() && is_keeled_group();
632
633    if (!src_normal && !src_keeled) {
634        error = "Please select a valid source group";
635    }
636    else {
637        if (dest->is_leaf()) {
638            error = GBS_global_string("'%s' is no valid destination for a group", dest->name);
639        }
640        else {
641            if (src_keeled) {
642                error = swap_group_with(dest, true);
643                if (error && src_normal) {
644                    GB_ERROR error1  = error;
645                    error            = swap_group_with(dest, false);
646                    if (error) error = GBS_global_string("Neighter keeled nor normal group can be moved that way:\n%s\n%s", error1, error);
647                }
648            }
649            else {
650                error = swap_group_with(dest, false);
651            }
652
653            if (!error) {
654                GBDATA *gb_retax = NULp;
655
656                if (!gb_retax && src_normal) gb_retax = find_group_name_entry(dest);
657                if (!gb_retax && src_keeled) gb_retax = find_group_name_entry(dest->get_father());
658                if (!gb_retax) gb_retax               = find_group_name_entry(this);
659                if (!gb_retax) gb_retax               = find_group_name_entry(get_father());
660
661                if (gb_retax) GB_touch(gb_retax); // force taxonomy reload
662                ap_assert(gb_retax); // should normally always find at least one name
663            }
664        }
665    }
666    return error;
667}
668
669#if defined(ASSERTION_USED) || defined(UNIT_TESTS)
670bool AP_tree::has_correct_mark_flags() const {
671    if (is_leaf()) return true;
672    if (!get_leftson() ->has_correct_mark_flags()) return false;
673    if (!get_rightson()->has_correct_mark_flags()) return false;
674
675    const AP_tree_members& left  = get_leftson()->gr;
676    const AP_tree_members& right = get_rightson()->gr;
677
678    unsigned wanted_mark_sum = left.mark_sum + right.mark_sum;
679    return gr.mark_sum == wanted_mark_sum;
680}
681#endif
682
683class AP_DefaultTreeShader: public AP_TreeShader {
684    // default tree shader (as used by unit-tests and ARB_PARSIMONY)
685
686    static void default_shader_never_shades() { ap_assert(0); }
687public:
688    AP_DefaultTreeShader() {}
689    void init() OVERRIDE {}
690    void update_settings() OVERRIDE {
691        colorize_marked = true;
692        colorize_groups = AW_color_groups_active();
693        shade_species   = false;
694    }
695
696    ShadedValue calc_shaded_leaf_GC(GBDATA */*gb_node*/) const OVERRIDE { default_shader_never_shades(); return NULp; }
697    ShadedValue calc_shaded_inner_GC(const ShadedValue& /*left*/, float /*left_ratio*/, const ShadedValue& /*right*/) const OVERRIDE { default_shader_never_shades(); return NULp; }
698    int to_GC(const ShadedValue& /*val*/) const OVERRIDE { default_shader_never_shades(); return -1; }
699};
700
701
702AP_TreeShader *AP_tree::shader = new AP_DefaultTreeShader;
703
704void AP_tree::set_tree_shader(AP_TreeShader *new_shader) {
705    delete shader;
706    shader = new_shader;
707    shader->init();
708}
709
710inline void AP_tree::recalc_hidden() {
711    // gr.hidden of father needs to be up-to-date!
712    gr.hidden = get_father() && (get_father()->gr.hidden || get_father()->gr.grouped);
713}
714
715inline void AP_tree::recalc_view_sum(const group_scaling& gscaling) {
716    // childs need to be up-to-date
717    // gr.leaf_sum needs to be up-to-date
718
719    ap_assert(!is_leaf());
720
721    if (is_folded_group()) {
722        ap_assert(gscaling.has_been_set());
723
724        const unsigned MIN_GROUP_SIZE = 2U;
725        unsigned       squared_size   = unsigned(pow(double(gr.leaf_sum), gscaling.pow)  * gscaling.linear);
726
727        gr.view_sum = std::max(squared_size, MIN_GROUP_SIZE);
728        gr.view_sum = std::min(gr.leaf_sum, gr.view_sum); // folded group will never use more space than unfolded
729    }
730    else {
731        gr.view_sum = get_leftson()->gr.view_sum + get_rightson()->gr.view_sum;
732    }
733}
734
735GB_ERROR AP_tree::update_and_write_folding(GBDATA *gb_tree, const group_scaling& gscaling) {
736    // Warning: only use if you know what you are doing!
737    //
738    // recalculates gr.hidden and gr.view_sum (after gr.grouped was modified)
739    // writes changed gr.grouped to database
740
741    GB_ERROR error = NULp;
742    recalc_hidden();
743    if (!is_leaf()) {
744        error = get_leftson()->update_and_write_folding(gb_tree, gscaling);
745        if (!error) error = get_rightson()->update_and_write_folding(gb_tree, gscaling);
746        if (!error) {
747            recalc_view_sum(gscaling);
748            error = tree_write_byte(gb_tree, this, gr.grouped, "grouped", 0);
749        }
750    }
751    return error;
752}
753
754void AP_tree::recompute_and_write_folding() {
755    // Warning: only use if you know what you are doing!
756    //
757    // Usable only if: folding changed (but nothing else!)
758    // Call with root-node of subtree for which folding shall be recalculated
759    // (needs to be the root of ALL folding changes!).
760    // Need to do resize afterwards.
761    //
762    // Note: gr.hidden is assumed to be correct for parent nodes!
763
764    AP_tree_root        *troot    = get_tree_root();
765    GBDATA              *gb_tree  = troot->get_gb_tree();
766    const group_scaling *gscaling = troot->get_group_scaling();
767
768    ap_assert(gb_tree);
769    ap_assert(gscaling);
770
771    GB_ERROR error = update_and_write_folding(gb_tree, *gscaling);
772    if (error) aw_message(GBS_global_string("Error in folding-update: %s", error));
773
774    // update view_sum of parent-nodes
775    {
776        AP_tree *fath = get_father();
777        while (fath) {
778            fath->recalc_view_sum(*gscaling);
779            fath = fath->get_father();
780        }
781    }
782}
783
784
785template<>
786ShadedValue AP_tree::update_subtree_information(const group_scaling& gscaling) {
787    ShadedValue val;
788    recalc_hidden();
789
790    if (is_leaf()) {
791        gr.view_sum = 1;
792        gr.leaf_sum = 1;
793
794        gr.max_tree_depth = 0.0;
795        gr.min_tree_depth = 0.0;
796
797        bool is_marked = gb_node && GB_read_flag(gb_node);
798
799        gr.mark_sum = int(is_marked);
800
801        gr.gc = shader->calc_leaf_GC(gb_node, is_marked);
802        if (shader->does_shade()) {
803            val = shader->calc_shaded_leaf_GC(gb_node);
804            if (gr.gc == AWT_GC_NONE_MARKED) {
805                gr.gc = shader->to_GC(val);
806            }
807        }
808    }
809    else {
810        ShadedValue leftVal  = get_leftson()->update_subtree_information<ShadedValue>(gscaling);
811        ShadedValue rightVal = get_rightson()->update_subtree_information<ShadedValue>(gscaling);
812
813        {
814            AP_tree_members& left  = get_leftson()->gr;
815            AP_tree_members& right = get_rightson()->gr;
816
817            gr.leaf_sum = left.leaf_sum + right.leaf_sum;
818
819            recalc_view_sum(gscaling);
820
821            gr.min_tree_depth = std::min(leftlen+left.min_tree_depth, rightlen+right.min_tree_depth);
822            gr.max_tree_depth = std::max(leftlen+left.max_tree_depth, rightlen+right.max_tree_depth);
823
824            gr.mark_sum = left.mark_sum + right.mark_sum;
825
826            gr.gc = shader->calc_inner_GC(left.gc, right.gc);
827            if (shader->does_shade()) {
828                float left_weight = left.leaf_sum / float(gr.leaf_sum);
829                val               = shader->calc_shaded_inner_GC(leftVal, left_weight, rightVal);
830                if (gr.gc == AWT_GC_NONE_MARKED) {
831                    gr.gc = shader->to_GC(val);
832                }
833            }
834        }
835    }
836    ap_assert(implicated(shader->does_shade(), val.isSet())); // expect we have shaded value (if shading is performed)
837    return val;
838}
839
840unsigned AP_tree::count_leafs() const {
841    return is_leaf()
842        ? 1
843        : get_leftson()->count_leafs() + get_rightson()->count_leafs();
844}
845
846int AP_tree::colorize(GB_HASH *hashptr) { // currently only used for multiprobes
847    // colorizes the tree according to 'hashptr'
848    // hashkey   = species id
849    // hashvalue = gc
850
851    int res;
852    if (is_leaf()) {
853        if (gb_node) {
854            res = GBS_read_hash(hashptr, name);
855            if (!res && GB_read_flag(gb_node)) { // marked but not in hash -> black
856                res = AWT_GC_BLACK;
857            }
858        }
859        else {
860            res = AWT_GC_ONLY_ZOMBIES;
861        }
862    }
863    else {
864        int l = get_leftson()->colorize(hashptr);
865        int r = get_rightson()->colorize(hashptr);
866
867        if      (l == r)                   res = l;
868        else if (l == AWT_GC_ONLY_ZOMBIES) res = r;
869        else if (r == AWT_GC_ONLY_ZOMBIES) res = l;
870        else                               res = AWT_GC_SOME_MARKED;
871    }
872    gr.gc = res;
873    return res;
874}
875
876void AP_tree::compute_tree() {
877#if defined(DEVEL_RALF) && 0
878    fputs(" - AP_tree::compute_tree() called\n", stderr);
879#endif
880    AP_tree_root        *troot    = get_tree_root();
881    const group_scaling *gscaling = troot->get_group_scaling();
882    ap_assert(gscaling && gscaling->has_been_set());
883
884    {
885        GB_transaction ta(troot->get_gb_main());
886        shader->update_settings();
887        update_subtree_information<ShadedValue>(*gscaling);
888    }
889}
890
891GB_ERROR AP_tree_root::loadFromDB(const char *name) {
892    GB_ERROR error = ARB_seqtree_root::loadFromDB(name);
893    if (!error) {
894        get_root_node()->load_subtree_info();
895    }
896    update_timers(); // maybe after link() ? // @@@ really do if error?
897    return error;
898}
899
900GB_ERROR AP_tree::relink() {
901    GB_transaction ta(get_tree_root()->get_gb_main()); // open close a transaction
902    GB_ERROR error = GBT_link_tree(this, get_tree_root()->get_gb_main(), false, NULp, NULp); // no status
903    get_tree_root()->update_timers();
904    return error;
905}
906
907AP_UPDATE_FLAGS AP_tree_root::check_update() {
908    GBDATA *gb_main = get_gb_main();
909    if (!gb_main) {
910        return AP_UPDATE_RELOADED;
911    }
912    else {
913        GB_transaction ta(gb_main);
914
915        if (is_tree_updated()) return AP_UPDATE_RELOADED;
916        if (is_species_updated()) return AP_UPDATE_RELINKED;
917        return AP_UPDATE_OK;
918    }
919}
920
921void AP_tree::buildLeafList_rek(AP_tree **list, long& num) {
922    // builds a list of all species
923    if (!is_leaf()) {
924        get_leftson()->buildLeafList_rek(list, num);
925        get_rightson()->buildLeafList_rek(list, num);
926    }
927    else {
928        list[num] = this;
929        num++;
930    }
931}
932
933void AP_tree::buildLeafList(AP_tree **&list, long &num) {
934    num        = count_leafs();
935    list       = new AP_tree *[num+1];
936    list[num]  = NULp;
937    long count = 0;
938
939    buildLeafList_rek(list, count);
940
941    ap_assert(count == num);
942}
943
944long AP_tree_root::remove_leafs(AWT_RemoveType awt_remove_type) {
945    // may remove the complete tree
946
947    ASSERT_VALID_TREE(get_root_node());
948
949    AP_tree **list;
950    long      count;
951    get_root_node()->buildLeafList(list, count);
952
953    GB_transaction ta(get_gb_main());
954    long removed = 0;
955
956    for (long i=0; i<count; i++) {
957        bool     removeNode = false;
958        AP_tree *leaf       = list[i];
959
960        if (leaf->gb_node) {
961            if ((awt_remove_type & AWT_REMOVE_NO_SEQUENCE) && !leaf->get_seq()) {
962                removeNode = true;
963            }
964            else if (awt_remove_type & (AWT_REMOVE_MARKED|AWT_REMOVE_UNMARKED)) {
965                long flag  = GB_read_flag(list[i]->gb_node);
966                removeNode = (flag && (awt_remove_type&AWT_REMOVE_MARKED)) || (!flag && (awt_remove_type&AWT_REMOVE_UNMARKED));
967            }
968        }
969        else {
970            if (awt_remove_type & AWT_REMOVE_ZOMBIES) {
971                removeNode = true;
972            }
973        }
974
975        if (removeNode) {
976            destroyNode(list[i]->REMOVE());
977            removed++;
978            if (!get_root_node()) {
979                break; // tree has been deleted
980            }
981        }
982        ASSERT_VALID_TREE(get_root_node());
983    }
984    delete [] list;
985
986    ASSERT_VALID_TREE_OR_NULL(get_root_node());
987    return removed;
988}
989
990// --------------------------------------------------------------------------------
991
992template <typename T>
993class ValueCounter {
994    T   min, max, sum;
995    int count;
996
997    char *mean_min_max_impl() const;
998    char *mean_min_max_percent_impl() const;
999
1000    mutable char *buf;
1001    const char *set_buf(char *content) const { freeset(buf, content); return buf; }
1002
1003public:
1004    ValueCounter()
1005        : min(INT_MAX),
1006          max(INT_MIN),
1007          sum(0),
1008          count(0),
1009          buf(NULp)
1010    {}
1011    ValueCounter(const ValueCounter<T>& other)
1012        : min(other.min),
1013          max(other.max),
1014          sum(other.sum),
1015          count(other.count),
1016          buf(NULp)
1017    {}
1018    ~ValueCounter() { free(buf); }
1019
1020    DECLARE_ASSIGNMENT_OPERATOR(ValueCounter<T>);
1021
1022    void count_value(T val) {
1023        count++;
1024        min  = std::min(min, val);
1025        max  = std::max(max, val);
1026        sum += val;
1027    }
1028
1029    int get_count() const { return count; }
1030    T get_min() const { return min; }
1031    T get_max() const { return max; }
1032    double get_mean() const { return sum/double(count); }
1033
1034    const char *mean_min_max()         const { return count ? set_buf(mean_min_max_impl())         : "<not available>"; }
1035    const char *mean_min_max_percent() const { return count ? set_buf(mean_min_max_percent_impl()) : "<not available>"; }
1036
1037    ValueCounter<T>& operator += (const T& inc) {
1038        min += inc;
1039        max += inc;
1040        sum += inc*count;
1041        return *this;
1042    }
1043    ValueCounter<T>& operator += (const ValueCounter<T>& other) {
1044        min    = std::min(min, other.min);
1045        max    = std::max(max, other.max);
1046        sum   += other.sum;
1047        count += other.count;
1048        return *this;
1049    }
1050};
1051
1052template<typename T>
1053inline ValueCounter<T> operator+(const ValueCounter<T>& c1, const ValueCounter<T>& c2) {
1054    return ValueCounter<T>(c1) += c2;
1055}
1056template<typename T>
1057inline ValueCounter<T> operator+(const ValueCounter<T>& c, const T& inc) {
1058    return ValueCounter<T>(c) += inc;
1059}
1060
1061template<> char *ValueCounter<int>::mean_min_max_impl() const {
1062    return GBS_global_string_copy("%.2f (range: %i .. %i)", get_mean(), get_min(), get_max());
1063}
1064template<> char *ValueCounter<double>::mean_min_max_impl() const {
1065    return GBS_global_string_copy("%.2f (range: %.2f .. %.2f)", get_mean(), get_min(), get_max());
1066}
1067template<> char *ValueCounter<double>::mean_min_max_percent_impl() const {
1068    return GBS_global_string_copy("%.2f%% (range: %.2f%% .. %.2f%%)", get_mean()*100.0, get_min()*100.0, get_max()*100.0);
1069}
1070
1071class LongBranchMarker {
1072    double min_rel_diff;
1073    double min_abs_diff;
1074
1075    int leafs;
1076    int nonzeroleafs;
1077    int multifurcs;
1078
1079    ValueCounter<double> absdiff;
1080    ValueCounter<double> reldiff;
1081    ValueCounter<double> absdiff_marked;
1082    ValueCounter<double> reldiff_marked;
1083
1084    double perform_marking(AP_tree *at, bool& marked) {
1085        marked = false;
1086        if (at->is_leaf()) {
1087            if (at->get_branchlength() != 0.0) {
1088                nonzeroleafs++;
1089            }
1090            leafs++;
1091            return 0.0;
1092        }
1093
1094        if (!at->is_root_node()) {
1095            if (at->get_branchlength() == 0.0) { // is multifurcation
1096                if (!at->get_father()->is_root_node() || at->is_leftson()) { // do not count two multifurcs @ sons of root
1097                    multifurcs++;
1098                }
1099            }
1100        }
1101
1102        bool   marked_left;
1103        bool   marked_right;
1104        double max = perform_marking(at->get_leftson(),  marked_left)  + at->leftlen;
1105        double min = perform_marking(at->get_rightson(), marked_right) + at->rightlen;
1106
1107        bool max_is_left = true;
1108        if (max<min) {
1109            double h = max; max = min; min = h;
1110            max_is_left = false;
1111        }
1112
1113        double abs_diff = max-min;
1114        absdiff.count_value(abs_diff);
1115
1116        double rel_diff = (max == 0.0) ? 0.0 : abs_diff/max;
1117        reldiff.count_value(rel_diff);
1118
1119        if (abs_diff>min_abs_diff && rel_diff>min_rel_diff) {
1120            if (max_is_left) {
1121                if (!marked_left) {
1122                    at->get_leftson()->mark_subtree();
1123                    marked = true;
1124                }
1125            }
1126            else {
1127                if (!marked_right) {
1128                    at->get_rightson()->mark_subtree();
1129                    marked = true;
1130                }
1131            }
1132        }
1133
1134        if (marked) { // just marked one of my subtrees
1135            absdiff_marked.count_value(abs_diff);
1136            reldiff_marked.count_value(rel_diff);
1137        }
1138        else {
1139            marked = marked_left||marked_right;
1140        }
1141
1142        return min; // use minimal distance for whole subtree
1143    }
1144
1145    static char *meanDiffs(const ValueCounter<double>& abs, const ValueCounter<double>& rel) {
1146        return GBS_global_string_copy(
1147            "Mean absolute diff: %s\n"
1148            "Mean relative diff: %s",
1149            abs.mean_min_max(), 
1150            rel.mean_min_max_percent());
1151    }
1152
1153public:
1154    LongBranchMarker(AP_tree *root, double min_rel_diff_, double min_abs_diff_)
1155        : min_rel_diff(min_rel_diff_),
1156          min_abs_diff(min_abs_diff_),
1157          leafs(0),
1158          nonzeroleafs(0),
1159          multifurcs(0)
1160    {
1161        bool UNUSED;
1162        perform_marking(root, UNUSED);
1163    }
1164
1165    const char *get_report() const {
1166        char *diffs_all    = meanDiffs(absdiff, reldiff);
1167        char *diffs_marked = meanDiffs(absdiff_marked, reldiff_marked);
1168
1169        int nodes     = leafs_2_nodes(leafs, UNROOTED);
1170        int edges     = nodes_2_edges(nodes);
1171        int zeroleafs = leafs-nonzeroleafs;
1172        int zeroedges = multifurcs+zeroleafs;
1173        int realedges = edges-zeroedges;
1174        int furcs     = nodes-leafs;    // = inner nodes
1175        int realfurcs = furcs-multifurcs;
1176
1177        int node_digits = calc_digits(nodes);
1178
1179        ap_assert(zeroleafs<=leafs);
1180        ap_assert(zeroedges<=edges);
1181        ap_assert(realedges<=edges);
1182        ap_assert(multifurcs<=furcs);
1183        ap_assert(realfurcs<=furcs);
1184
1185        const char *msg = GBS_global_string(
1186            "Unrooted tree contains %*i furcations,\n"
1187            "              of which %*i are multifurcations,\n"
1188            "                  i.e. %*i are \"real\" furcations.\n"
1189            "\n"
1190            "Unrooted tree contains %*i edges,\n"
1191            "              of which %*i have a length > zero.\n"
1192            "\n"
1193            "%s\n"
1194            "\n"
1195            "%i subtrees have been marked:\n"
1196            "%s\n"
1197            "\n",
1198            node_digits, furcs,
1199            node_digits, multifurcs,
1200            node_digits, realfurcs,
1201            node_digits, edges,
1202            node_digits, realedges,
1203            diffs_all,
1204            absdiff_marked.get_count(),
1205            diffs_marked);
1206
1207        free(diffs_all);
1208        free(diffs_marked);
1209
1210        return msg;
1211    }
1212
1213    double get_max_abs_diff() const { return absdiff.get_max(); }
1214};
1215
1216struct DepthMarker {
1217    // limits (marked if depth and dist are above)
1218    int    min_depth;
1219    double min_rootdist;
1220
1221    // current values (for recursion)
1222    int    depth;
1223    double dist;
1224
1225    // results
1226    ValueCounter<int>    depths,    depths_marked;
1227    ValueCounter<double> distances, distances_marked;
1228
1229    void perform_marking(AP_tree *at, AP_FLOAT atLen) {
1230        int depthInc = atLen == 0.0 ? 0 : 1;  // do NOT increase depth at multifurcations
1231
1232        depth += depthInc;
1233        dist  += atLen;
1234
1235        if (at->is_leaf()) {
1236            depths.count_value(depth);
1237            distances.count_value(dist);
1238
1239            int mark = depth >= min_depth && dist >= min_rootdist;
1240            if (at->gb_node) {
1241                GB_write_flag(at->gb_node, mark);
1242                if (mark) {
1243                    depths_marked.count_value(depth);
1244                    distances_marked.count_value(dist);
1245                }
1246            }
1247        }
1248        else {
1249            perform_marking(at->get_leftson(), at->leftlen);
1250            perform_marking(at->get_rightson(), at->rightlen);
1251        }
1252
1253        depth -= depthInc;
1254        dist  -= atLen;
1255    }
1256
1257public:
1258    DepthMarker(AP_tree *root, int min_depth_, double min_rootdist_)
1259        : min_depth(min_depth_),
1260          min_rootdist(min_rootdist_),
1261          depth(0),
1262          dist(0.0)
1263    {
1264        perform_marking(root, 0.0);
1265    }
1266
1267    const char *get_report() const {
1268        int    leafs          = depths.get_count();
1269        int    marked         = depths_marked.get_count();
1270        double balanced_depth = log10(leafs) / log10(2);
1271
1272        const char *msg = GBS_global_string(
1273            "The optimal mean depth of a tree with %i leafs\n"
1274            "      would be %.2f\n"
1275            "\n"
1276            "Your tree:\n"
1277            "mean depth:    %s\n"
1278            "mean distance: %s\n"
1279            "\n"
1280            "%i species (%.2f%%) have been marked:\n"
1281            "mean depth:    %s\n"
1282            "mean distance: %s\n"
1283            ,
1284            leafs,
1285            balanced_depth,
1286            depths.mean_min_max(),
1287            distances.mean_min_max(), 
1288            marked, marked/double(leafs)*100.0,
1289            depths_marked.mean_min_max(), 
1290            distances_marked.mean_min_max()
1291            );
1292        return msg;
1293    }
1294
1295    int get_max_depth() const { return depths.get_max(); }
1296    double get_max_rootdist() const { return distances.get_max(); }
1297};
1298
1299const char *AP_tree::mark_long_branches(double min_rel_diff, double min_abs_diff, double& found_max_abs_diff) {
1300    // look for asymmetric parts of the tree and mark all species with long branches
1301    LongBranchMarker lmarker(this, min_rel_diff, min_abs_diff);
1302    found_max_abs_diff = lmarker.get_max_abs_diff();
1303    return lmarker.get_report();
1304}
1305const char *AP_tree::mark_deep_leafs(int min_depth, double min_rootdist, int& found_max_depth, double& found_max_rootdist) {
1306    // mark all leafs with min_depth and min_rootdist
1307    DepthMarker dmarker(this, min_depth, min_rootdist);
1308    found_max_depth    = dmarker.get_max_depth();
1309    found_max_rootdist = dmarker.get_max_rootdist();
1310    return dmarker.get_report();
1311}
1312
1313// --------------------------------------------------------------------------------
1314
1315typedef ValueCounter<double> Distance;
1316
1317class DistanceCounter {
1318    Distance min, max, mean;
1319public:
1320
1321    void count_distance(const Distance& d) {
1322        mean.count_value(d.get_mean());
1323        min.count_value(d.get_min());
1324        max.count_value(d.get_max());
1325    }
1326
1327    int get_count() const { return mean.get_count(); }
1328
1329    char *get_report() const {
1330        return GBS_global_string_copy(
1331            "Mean mean distance: %s\n"
1332            "Mean min. distance: %s\n"
1333            "Mean max. distance: %s",
1334            mean.mean_min_max(), 
1335            min.mean_min_max(), 
1336            max.mean_min_max()
1337            );
1338    }
1339};
1340
1341class EdgeDistances {
1342    typedef map<AP_tree*, Distance> DistanceMap;
1343
1344    DistanceMap downdist; // inclusive length of branch itself
1345    DistanceMap updist;   // inclusive length of branch itself
1346
1347    GBT_LEN distSum; // of all distances in tree
1348
1349    arb_progress progress;
1350
1351    const Distance& calc_downdist(AP_tree *at, AP_FLOAT len) {
1352        if (at->is_leaf()) {
1353            Distance d;
1354            d.count_value(len);
1355            downdist[at] = d;
1356
1357            progress.inc();
1358        }
1359        else {
1360            downdist[at] =
1361                calc_downdist(at->get_leftson(), at->leftlen) +
1362                calc_downdist(at->get_rightson(), at->rightlen) +
1363                len;
1364        }
1365        return downdist[at];
1366    }
1367
1368    const Distance& calc_updist(AP_tree *at, AP_FLOAT len) {
1369        ap_assert(at->father); // impossible - root has no updist!
1370
1371        AP_tree *father  = at->get_father();
1372        AP_tree *brother = at->get_brother();
1373
1374        if (father->father) {
1375            ap_assert(updist.find(father) != updist.end());
1376            ap_assert(downdist.find(brother) != downdist.end());
1377
1378            updist[at] = updist[father] + downdist[brother] + len;
1379        }
1380        else {
1381            ap_assert(downdist.find(brother) != downdist.end());
1382
1383            updist[at] = downdist[brother]+len;
1384        }
1385
1386        if (!at->is_leaf()) {
1387            calc_updist(at->get_leftson(), at->leftlen);
1388            calc_updist(at->get_rightson(), at->rightlen);
1389        }
1390        else {
1391            progress.inc();
1392        }
1393       
1394        return updist[at];
1395    }
1396
1397    DistanceCounter alldists, markeddists;
1398
1399    void calc_distance_stats(AP_tree *at) {
1400        if (at->is_leaf()) {
1401            ap_assert(updist.find(at) != updist.end());
1402
1403            const Distance& upwards = updist[at];
1404
1405            alldists.count_distance(upwards);
1406            if (at->gb_node && GB_read_flag(at->gb_node)) {
1407                markeddists.count_distance(upwards);
1408            }
1409
1410            progress.inc();
1411        }
1412        else {
1413            calc_distance_stats(at->get_leftson());
1414            calc_distance_stats(at->get_rightson());
1415        }
1416    }
1417
1418public:
1419
1420    EdgeDistances(AP_tree *root)
1421        : distSum(root->sum_child_lengths()),
1422          progress("Analysing distances", root->count_leafs()*3L)
1423    {
1424        calc_downdist(root->get_leftson(),  root->leftlen);
1425        calc_downdist(root->get_rightson(), root->rightlen);
1426
1427        calc_updist(root->get_leftson(),  root->leftlen);
1428        calc_updist(root->get_rightson(), root->rightlen);
1429
1430        calc_distance_stats(root);
1431    }
1432
1433    const char *get_report() const {
1434        char *alldists_report    = alldists.get_report();
1435        char *markeddists_report = markeddists.get_report();
1436
1437        const char *msg = GBS_global_string(
1438            "Overall in-tree-distance (ITD): %.3f\n"
1439            "    per-species-distance (PSD): %.6f\n"
1440            "\n"
1441            "Distance statistic for %i leafs:\n"
1442            "(each leaf to all other leafs)\n"
1443            "\n"
1444            "%s\n"
1445            "\n"
1446            "Distance statistic for %i marked leafs:\n"
1447            "\n"
1448            "%s\n",
1449            distSum,
1450            distSum / alldists.get_count(),
1451            alldists.get_count(), alldists_report,
1452            markeddists.get_count(), markeddists_report);
1453
1454        free(markeddists_report);
1455        free(alldists_report);
1456
1457        return msg;
1458    }
1459};
1460
1461const char *AP_tree::analyse_distances() { return EdgeDistances(this).get_report(); }
1462
1463// --------------------------------------------------------------------------------
1464
1465static int ap_mark_degenerated(AP_tree *at, double degeneration_factor, double& max_degeneration) {
1466    // returns number of species in subtree
1467
1468    if (at->is_leaf()) return 1;
1469
1470    int lSons = ap_mark_degenerated(at->get_leftson(), degeneration_factor, max_degeneration);
1471    int rSons = ap_mark_degenerated(at->get_rightson(), degeneration_factor, max_degeneration);
1472
1473    double this_degeneration = 0;
1474
1475    if (lSons<rSons) {
1476        this_degeneration = rSons/double(lSons);
1477        if (this_degeneration >= degeneration_factor) {
1478            at->get_leftson()->mark_subtree();
1479        }
1480
1481    }
1482    else if (rSons<lSons) {
1483        this_degeneration = lSons/double(rSons);
1484        if (this_degeneration >= degeneration_factor) {
1485            at->get_rightson()->mark_subtree();
1486        }
1487    }
1488
1489    if (this_degeneration >= max_degeneration) {
1490        max_degeneration = this_degeneration;
1491    }
1492
1493    return lSons+rSons;
1494}
1495
1496double AP_tree::mark_degenerated_branches(double degeneration_factor) {
1497    // marks all species in degenerated branches.
1498    // For all nodes, where one branch contains 'degeneration_factor' more species than the
1499    // other branch, the smaller branch is considered degenerated.
1500
1501    double max_degeneration = 0;
1502    ap_mark_degenerated(this, degeneration_factor, max_degeneration);
1503    return max_degeneration;
1504}
1505
1506static int ap_mark_duplicates_rek(AP_tree *at, GB_HASH *seen_species) {
1507    if (at->is_leaf()) {
1508        if (at->name) {
1509            if (GBS_read_hash(seen_species, at->name)) { // already seen -> mark species
1510                if (at->gb_node) {
1511                    GB_write_flag(at->gb_node, 1);
1512                }
1513                else { // duplicated zombie
1514                    return 1;
1515                }
1516            }
1517            else { // first occurrence
1518                GBS_write_hash(seen_species, at->name, 1);
1519            }
1520        }
1521    }
1522    else {
1523        return
1524            ap_mark_duplicates_rek(at->get_leftson(), seen_species) +
1525            ap_mark_duplicates_rek(at->get_rightson(), seen_species);
1526    }
1527    return 0;
1528}
1529
1530void AP_tree::mark_duplicates() {
1531    GB_HASH *seen_species = GBS_create_hash(gr.leaf_sum, GB_IGNORE_CASE);
1532
1533    int dup_zombies = ap_mark_duplicates_rek(this, seen_species);
1534    if (dup_zombies) {
1535        aw_message(GBS_global_string("Warning: Detected %i duplicated zombies (can't mark them)", dup_zombies));
1536    }
1537    GBS_free_hash(seen_species);
1538}
1539
1540static double ap_just_tree_rek(AP_tree *at) {
1541    if (at->is_leaf()) {
1542        return 0.0;
1543    }
1544    else {
1545        double bl = ap_just_tree_rek(at->get_leftson());
1546        double br = ap_just_tree_rek(at->get_rightson());
1547
1548        double l = at->leftlen + at->rightlen;
1549        double diff = fabs(bl - br);
1550        if (l < diff * 1.1) l = diff * 1.1;
1551        double go = (bl + br + l) * .5;
1552        at->leftlen = go - bl;
1553        at->rightlen = go - br;
1554        return go;
1555    }
1556}
1557
1558
1559void AP_tree::justify_branch_lenghs(GBDATA *gb_main) {
1560    // shift branches to create a symmetric looking tree
1561    GB_transaction ta(gb_main);
1562    ap_just_tree_rek(this);
1563}
1564
1565static void relink_tree_rek(AP_tree *node, void (*relinker)(GBDATA *&ref_gb_node, char *&ref_name, GB_HASH *organism_hash), GB_HASH *organism_hash) {
1566    if (node->is_leaf()) {
1567        relinker(node->gb_node, node->name, organism_hash);
1568    }
1569    else {
1570        relink_tree_rek(node->get_leftson(), relinker, organism_hash);
1571        relink_tree_rek(node->get_rightson(), relinker, organism_hash);
1572    }
1573}
1574
1575void AP_tree::relink_tree(GBDATA *gb_main, void (*relinker)(GBDATA *&ref_gb_node, char *&ref_name, GB_HASH *organism_hash), GB_HASH *organism_hash) {
1576    // relinks the tree using a relinker-function
1577    // every node in tree is passed to relinker, relinker might modify
1578    // these values (ref_gb_node and ref_name) and the modified values are written back into tree
1579
1580    GB_transaction  ta(gb_main);
1581    relink_tree_rek(this, relinker, organism_hash);
1582}
1583
1584
1585void AP_tree::reset_child_angles() {
1586    if (!is_leaf()) {
1587        gr.reset_both_child_angles();
1588        get_leftson()->reset_child_angles();
1589        get_rightson()->reset_child_angles();
1590    }
1591}
1592
1593void AP_tree::reset_child_linewidths() {
1594    if (!is_leaf()) {
1595        gr.reset_both_child_linewidths();
1596        get_leftson()->reset_child_linewidths();
1597        get_rightson()->reset_child_linewidths();
1598    }
1599}
1600
1601void AP_tree::set_linewidth_recursive(int width) {
1602    set_linewidth(width);
1603    if (!is_leaf()) {
1604        get_leftson()->set_linewidth_recursive(width);
1605        get_rightson()->set_linewidth_recursive(width);
1606    }
1607}
1608
1609void AP_tree::reset_child_layout() {
1610    if (!is_leaf()) {
1611        gr.reset_child_spread();
1612        gr.reset_both_child_angles();
1613        gr.reset_both_child_linewidths();
1614        get_leftson()->reset_child_layout();
1615        get_rightson()->reset_child_layout();
1616    }
1617}
1618
1619void AP_tree::reset_subtree_spreads() {
1620    gr.reset_child_spread();
1621    if (!is_leaf()) {
1622        get_leftson()->reset_subtree_spreads();
1623        get_rightson()->reset_subtree_spreads();
1624    }
1625}
1626void AP_tree::reset_subtree_angles() {
1627    reset_angle();
1628    if (!is_leaf()) reset_child_angles();
1629}
1630void AP_tree::reset_subtree_linewidths() {
1631    reset_linewidth();
1632    if (!is_leaf()) reset_child_linewidths();
1633}
1634void AP_tree::reset_subtree_layout() {
1635    reset_linewidth();
1636    reset_angle();
1637    if (!is_leaf()) reset_child_layout();
1638}
1639
1640bool AP_tree::is_inside_folded_group() const {
1641    if (!is_leaf() && is_folded_group()) return true;
1642    if (!father) return false;
1643    return get_father()->is_inside_folded_group();
1644}
1645
Note: See TracBrowser for help on using the repository browser.