source: tags/old_import_filter/SL/AP_TREE/AP_Tree.cxx

Last change on this file was 9533, checked in by westram, 11 years ago
  • use explicit override
    • rest
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 55.5 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
13#include <AP_filter.hxx>
14#include <awt_attributes.hxx>
15#include <aw_msg.hxx>
16
17#include <math.h>
18#include <map>
19#include <climits>
20#include <arb_progress.h>
21
22using namespace std;
23
24/*!***************************************************************************************
25 ************           Rates                       **********
26 *****************************************************************************************/
27void AP_rates::print()
28{
29    AP_FLOAT max;
30    int i;
31
32    max = 0.0;
33    for (i=0; i<rate_len; i++) {
34        if (rates[i] > max) max = rates[i];
35    }
36    printf("rates:");
37    for (i=0; i<rate_len; i++) {
38        putchar('0' + (int)(rates[i]/max*9.9));
39    }
40    printf("\n");
41}
42
43AP_rates::AP_rates() {
44    memset ((char *)this, 0, sizeof(AP_rates));
45}
46
47char *AP_rates::init(AP_filter *fil) {
48    if (fil->get_timestamp() <= this->update) return 0;
49
50    rate_len = fil->get_filtered_length();
51    delete [] rates;
52    rates = new AP_FLOAT[rate_len];
53    for (int i=0; i<rate_len; i++) {
54        rates[i] = 1.0;
55    }
56    this->update = fil->get_timestamp();
57    return 0;
58}
59
60char *AP_rates::init(AP_FLOAT * ra, AP_filter *fil)
61{
62    if (fil->get_timestamp() <= this->update) return 0;
63
64    rate_len = fil->get_filtered_length();
65    delete [] rates;
66    rates = new AP_FLOAT[rate_len];
67    int i, j;
68    for (j=i=0; i<rate_len; j++) {
69        if (fil->use_position(j)) {
70            rates[i++] = ra[j];
71        }
72    }
73    this->update = fil->get_timestamp();
74    return 0;
75}
76
77AP_rates::~AP_rates() { delete [] rates; }
78
79
80/*!***************************************************************************************
81 ************           AP_tree_root                    **********
82 *****************************************************************************************/
83
84AP_tree_root::AP_tree_root(AliView *aliView, const AP_tree& tree_proto, AP_sequence *seq_proto, bool add_delete_callbacks)
85    : ARB_tree_root(aliView, tree_proto, seq_proto, add_delete_callbacks)
86    , root_changed_cb(NULL), root_changed_cd(NULL)
87    , node_deleted_cb(NULL), node_deleted_cd(NULL)
88    , gb_tree_gone(NULL)
89    , tree_timer(0)
90    , species_timer(0)
91    , table_timer(0)
92    , rates(NULL)
93{
94    GBDATA         *gb_main = get_gb_main();
95    GB_transaction  ta(gb_main);
96
97    gb_species_data = GBT_get_species_data(gb_main);
98    gb_table_data   = GB_search(gb_main, "table_data", GB_CREATE_CONTAINER);
99}
100
101AP_tree_root::~AP_tree_root() {
102    delete get_root_node();
103    ap_assert(!get_root_node());
104}
105
106bool AP_tree_root::is_tree_updated() {
107    GBDATA *gbtree = get_gb_tree();
108    if (gbtree) {
109        GB_transaction ta(gbtree);
110        return GB_read_clock(gbtree)>tree_timer;
111    }
112    return true;
113}
114
115bool AP_tree_root::is_species_updated() {
116    if (gb_species_data) {
117        GB_transaction ta(gb_species_data);
118        return
119            GB_read_clock(gb_species_data)>species_timer ||
120            GB_read_clock(gb_table_data)>table_timer;
121    }
122    return true;
123}
124
125void AP_tree_root::update_timers() {
126    if (gb_species_data) {
127        GB_transaction  dummy(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        table_timer            = GB_read_clock(gb_table_data);
132    }
133}
134
135/*!***************************************************************************************
136 ************           AP_tree                     **********
137 *****************************************************************************************/
138static void ap_tree_node_deleted(GBDATA *, int *cl, GB_CB_TYPE) {
139    AP_tree *THIS = (AP_tree *)cl;
140    THIS->gb_node = 0;
141}
142
143AP_tree::AP_tree(AP_tree_root *tree_rooti)
144    : ARB_tree(tree_rooti)
145    , stack_level(0)
146{
147    gr.clear();
148    gr.spread = 1.0;
149    br.clear();
150}
151
152AP_tree::~AP_tree() {
153    if (gr.callback_exists && gb_node) {
154        GB_remove_callback(gb_node, GB_CB_DELETE, ap_tree_node_deleted, (int *)this);
155    }
156
157    AP_tree_root *root = get_tree_root();
158    if (root) root->inform_about_delete(this);
159}
160
161AP_tree *AP_tree::dup() const {
162    return new AP_tree(const_cast<AP_tree_root*>(get_tree_root()));
163}
164
165#if defined(WARN_TODO)
166#warning move to ARB_tree ?
167#endif
168void AP_tree::replace_self(AP_tree *new_son) {
169    ap_assert(father);
170    if (father) {
171        if (is_leftson(father)) father->leftson  = new_son;
172        else                    father->rightson = new_son;
173    }
174}
175#if defined(WARN_TODO)
176#warning move to ARB_tree ?
177#endif
178void AP_tree::set_brother(AP_tree *new_son) {
179    ap_assert(father);
180    if (father) {
181        if (is_leftson(father)) father->rightson = new_son;
182        else                    father->leftson  = new_son;
183    }
184}
185
186void AP_tree::clear_branch_flags() {
187    br.clear();
188    if (!is_leaf) {
189        get_leftson()->clear_branch_flags();
190        get_rightson()->clear_branch_flags();
191    }
192}
193
194void AP_tree::insert(AP_tree *new_brother) {
195    ASSERT_VALID_TREE(this);
196    ASSERT_VALID_TREE(new_brother);
197
198    AP_tree  *new_tree = dup();
199    AP_FLOAT  laenge;
200
201    new_tree->leftson  = this;
202    new_tree->rightson = new_brother;
203    new_tree->father   = new_brother->father;
204    father             = new_tree;
205
206    if (new_brother->father) {
207        if (new_brother->father->leftson == new_brother) {
208            laenge = new_brother->father->leftlen = new_brother->father->leftlen*.5;
209            new_brother->father->leftson = new_tree;
210        }
211        else {
212            laenge = new_brother->father->rightlen = new_brother->father->rightlen*.5;
213            new_brother->father->rightson = new_tree;
214        }
215    }
216    else {
217        laenge = 0.5;
218    }
219    new_tree->leftlen   = laenge;
220    new_tree->rightlen  = laenge;
221    new_brother->father = new_tree;
222
223    AP_tree_root *troot = new_brother->get_tree_root();
224    if (troot) {
225        if (!new_tree->father) troot->change_root(new_brother, new_tree);
226        else new_tree->set_tree_root(troot);
227
228        ASSERT_VALID_TREE(troot->get_root_node());
229    }
230    else { // loose nodes (not in a tree)
231        ASSERT_VALID_TREE(new_tree);
232    }
233}
234
235#if defined(WARN_TODO)
236#warning move to ARB_tree ?
237#endif
238void AP_tree_root::change_root(ARB_tree *oldroot, ARB_tree *newroot) {
239    if (root_changed_cb) {
240        root_changed_cb(root_changed_cd, DOWNCAST(AP_tree*, oldroot), DOWNCAST(AP_tree*, newroot));
241    }
242    if (!oldroot) {
243        ap_assert(newroot);
244        if (gb_tree_gone) {                         // when tree was temporarily deleted (e.g. by 'Remove & add all')
245            set_gb_tree(gb_tree_gone);              // re-use previous DB entry
246            gb_tree_gone = NULL;
247        }
248    }
249    if (!newroot) {                                 // tree empty
250        GBDATA *gbtree = get_gb_tree();
251        if (gbtree) {
252            ap_assert(gb_tree_gone == NULL);        // no tree should be remembered yet
253            gb_tree_gone = gbtree;                  // remember for deletion (done in AP_tree::save)
254        }
255    }
256    ARB_tree_root::change_root(oldroot, newroot);
257}
258
259void AP_tree_root::inform_about_delete(AP_tree *del) {
260    if (node_deleted_cb) node_deleted_cb(node_deleted_cd, del);
261}
262
263void AP_tree_root::set_root_changed_callback(AP_rootChangedCb cb, void *cd) {
264    root_changed_cb = cb;
265    root_changed_cd = cd;
266}
267
268void AP_tree_root::set_node_deleted_callback(AP_nodeDelCb cb, void *cd) {
269    node_deleted_cb = cb;
270    node_deleted_cd = cd;
271}
272
273
274void AP_tree::remove() {
275    // remove this + father from tree
276    // Note: does not delete this or father!
277
278    ASSERT_VALID_TREE(this);
279    if (father == 0) {
280        get_tree_root()->change_root(this, NULL); // wot ?
281    }
282    else {
283        AP_tree      *brother     = get_brother();  // brother remains in tree
284        GBT_LEN       brothersLen = brother->get_branchlength();
285        AP_tree      *fath        = get_father();   // fath of this is removed as well
286        ARB_tree     *grandfather = fath->father;
287        AP_tree_root *troot       = get_tree_root();
288
289        if (fath->gb_node) {                      // move inner information to remaining subtree
290            if (!brother->gb_node && !brother->is_leaf) {
291                brother->gb_node = fath->gb_node;
292                fath->gb_node    = 0;
293            }
294        }
295
296        if (grandfather) {
297            brother->unlink_from_father();
298
299            bool isLeftSon = fath->is_leftson(grandfather);
300            fath->unlink_from_father();
301
302            if (isLeftSon) {
303                ap_assert(!grandfather->leftson);
304                grandfather->leftlen += brothersLen;
305                grandfather->leftson  = brother;
306            }
307            else {
308                ap_assert(!grandfather->rightson);
309                grandfather->rightlen += brothersLen;
310                grandfather->rightson  = brother;
311            }
312            brother->father = grandfather;
313        }
314        else {                                      // father is root, make brother the new root
315            if (brother->is_leaf) {
316                troot->change_root(fath, NULL);     // erase tree from root
317            }
318            else {
319                brother->unlink_from_father();
320                troot->change_root(fath, brother);
321            }
322        }
323
324        ap_assert(fath == father);
325
326#if defined(CHECK_TREE_STRUCTURE)
327        if (troot->get_root_node()) ASSERT_VALID_TREE(troot->get_root_node());
328#endif // CHECK_TREE_STRUCTURE
329
330        troot->inform_about_delete(fath);
331        troot->inform_about_delete(this);
332    }
333}
334
335GB_ERROR AP_tree::cantMoveTo(AP_tree *new_brother) {
336    GB_ERROR error = 0;
337
338    if (!father)                                error = "Can't move the root of the tree";
339    else if (!new_brother->father)              error = "Can't move to the root of the tree";
340    else if (new_brother->father == father)     error = "Already there";
341    else if (!father->father)                   error = "Can't move son of root";
342    else if (new_brother->is_inside(this))      error = "Can't move a subtree into itself";
343
344    return error;
345}
346
347void AP_tree::moveTo(AP_tree *new_brother, AP_FLOAT rel_pos) {
348    // rel_pos == 0.0 -> at father
349    //         == 1.0 -> at brother
350
351    ap_assert(father);
352    ap_assert(new_brother->father);
353    ap_assert(new_brother->father != father);       // already there
354    ap_assert(!new_brother->is_inside(this));       // can't move tree into itself
355
356    if (father->leftson != this) get_father()->swap_sons();
357
358    if (father->father == 0) {
359        get_brother()->father = 0;
360        get_tree_root()->change_root(get_father(), get_brother());
361    }
362    else {
363        ARB_tree *grandfather = father->father;
364        if (father == new_brother) {    // just pull branches !!
365            new_brother  = get_brother();
366            if (grandfather->leftson == father) {
367                rel_pos *= grandfather->leftlen / (father->rightlen+grandfather->leftlen);
368            }
369            else {
370                rel_pos *= grandfather->rightlen / (father->rightlen+grandfather->rightlen);
371            }
372        }
373        else if (new_brother->father == father) { // just pull branches !!
374            rel_pos =
375                1.0 + (rel_pos-1.0) * father->rightlen
376                /
377                (father->rightlen + (grandfather->leftson == father ? grandfather->leftlen : grandfather->rightlen));
378        }
379
380        if (grandfather->leftson == father) {
381            grandfather->leftlen     += father->rightlen;
382            grandfather->leftson      = father->rightson;
383            father->rightson->father  = grandfather;
384        }
385        else {
386            grandfather->rightlen    += father->rightlen;
387            grandfather->rightson     = father->rightson;
388            father->rightson->father  = grandfather;
389        }
390    }
391
392    ARB_tree *new_tree          = father;
393    ARB_tree *brother_father    = new_brother->father;
394    AP_FLOAT  laenge;
395
396    if (brother_father->leftson == new_brother) {
397        laenge  = brother_father->leftlen;
398        laenge -= brother_father->leftlen = laenge * rel_pos;
399        new_brother->father->leftson = new_tree;
400    }
401    else {
402        laenge  = brother_father->rightlen;
403        laenge -= brother_father->rightlen = laenge * rel_pos;
404        brother_father->rightson = new_tree;
405    }
406
407    new_tree->rightlen  = laenge;
408    new_brother->father = new_tree;
409    new_tree->rightson  = new_brother;
410    new_tree->father    = brother_father;
411}
412
413void AP_tree::swap_sons() {
414    ARB_tree *h_at = this->leftson;
415    this->leftson = this->rightson;
416    this->rightson = h_at;
417
418    double h = this->leftlen;
419    this->leftlen = this->rightlen;
420    this->rightlen = h;
421}
422
423void AP_tree::swap_assymetric(AP_TREE_SIDE mode) {
424    // mode AP_LEFT exchanges lefts with brother
425    // mode AP_RIGHT exchanges rights with brother
426
427    ap_assert(!is_leaf);                           // cannot swap leafs
428    ap_assert(father);                             // cannot swap root (has no brother)
429    ap_assert(mode == AP_LEFT || mode == AP_RIGHT); // illegal mode
430
431    ARB_tree *pntr;
432
433    if (father->father == 0) { // father is root
434        AP_tree *pntr_brother = get_brother();
435        if (!pntr_brother->is_leaf) {
436            if (mode == AP_LEFT) {
437                pntr_brother->leftson->father = this;
438                pntr                          = pntr_brother->leftson;
439                pntr_brother->leftson         = leftson;
440                leftson->father               = pntr_brother;
441                leftson                       = pntr;
442            }
443            else {
444                pntr_brother->leftson->father = this;
445                rightson->father              = pntr_brother;
446                pntr                          = pntr_brother->leftson;
447                pntr_brother->leftson         = rightson;
448                rightson                      = pntr;
449            }
450        }
451    }
452    else {
453        if (mode == AP_LEFT) { // swap leftson with brother
454            if (father->leftson == this) {
455                father->rightson->father = this;
456                leftson->father          = father;
457                pntr                     = father->rightson;
458                AP_FLOAT help_len        = father->rightlen;
459                father->rightlen         = leftlen;
460                leftlen                  = help_len;
461                father->rightson         = leftson;
462                leftson                  = pntr;
463            }
464            else {
465                father->leftson->father = this;
466                leftson->father         = father;
467                pntr                    = father->leftson;
468                AP_FLOAT help_len       = father->leftlen;
469                father->leftlen         = leftlen;
470                leftlen                 = help_len;
471                father->leftson         = leftson;
472                leftson                 = pntr;
473            }
474        }
475        else { // swap rightson with brother
476            if (father->leftson == this) {
477                father->rightson->father = this;
478                rightson->father         = father;
479                pntr                     = father->rightson;
480                AP_FLOAT help_len        = father->rightlen;
481                father->rightlen         = rightlen;
482                rightlen                 = help_len;
483                father->rightson         = rightson;
484                rightson                 = pntr;
485            }
486            else {
487                father->leftson->father = this;
488                rightson->father        = father;
489                pntr                    = father->leftson;
490                AP_FLOAT help_len       = father->leftlen;
491                father->leftson         = rightson;
492                father->leftlen         = rightlen;
493                rightlen                = help_len;
494                rightson                = pntr;
495            }
496        }
497    }
498}
499
500void AP_tree::set_root() {
501    if (!father) return;                            // already root
502    if (!father->father) return;                    // already root?
503
504    AP_tree *old_root    = 0;
505    AP_tree *old_brother = 0;
506    {
507        AP_branch_members  br1 = br;
508        AP_tree           *pntr;
509
510        for  (pntr = get_father(); pntr->father; pntr = pntr->get_father()) {
511            AP_branch_members br2 = pntr->br;
512            pntr->br              = br1;
513            br1                   = br2;
514            old_brother           = pntr;
515        }
516        if (pntr->leftson == old_brother) {
517            pntr->get_rightson()->br = br1;
518        }
519        old_root = pntr;
520    }
521    if (old_brother) old_brother = old_brother->get_brother();
522
523    {
524        // move remark branches to top
525        ARB_tree *node;
526        char     *remark = nulldup(remark_branch);
527
528        for (node = this; node->father; node = node->father) {
529            char *sh            = node->remark_branch;
530            node->remark_branch = remark;
531            remark              = sh;
532        }
533        delete remark;
534    }
535    AP_FLOAT old_root_len = old_root->leftlen + old_root->rightlen;
536
537    // new node & this init
538
539    old_root->leftson  = this;
540    old_root->rightson = father;                    // will be set later
541
542    if (father->leftson == this) {
543        old_root->leftlen = old_root->rightlen = father->leftlen*.5;
544    }
545    else {
546        old_root->leftlen = old_root->rightlen = father->rightlen*.5;
547    }
548
549    ARB_tree *next = father->father;
550    ARB_tree *prev = old_root;
551    ARB_tree *pntr = father;
552
553    if (father->leftson == this)    father->leftson = old_root; // to set the flag correctly
554
555    // loop from father to son of root, rotate tree
556    while  (next->father) {
557        double len = (next->leftson == pntr) ? next->leftlen : next->rightlen;
558
559        if (pntr->leftson == prev) {
560            pntr->leftson = next;
561            pntr->leftlen = len;
562        }
563        else {
564            pntr->rightson = next;
565            pntr->rightlen = len;
566        }
567
568        pntr->father = prev;
569        prev         = pntr;
570        pntr         = next;
571        next         = next->father;
572    }
573    // now next points to the old root, which has been destroyed above
574    //
575    // pointer at oldroot
576    // pntr == brother before old root == next
577
578    if (pntr->leftson == prev) {
579        pntr->leftlen = old_root_len;
580        pntr->leftson = old_brother;
581    }
582    else {
583        pntr->rightlen = old_root_len;
584        pntr->rightson = old_brother;
585    }
586
587    old_brother->father = pntr;
588    pntr->father        = prev;
589    father              = old_root;
590}
591
592
593inline short tree_read_byte(GBDATA *tree, const char *key, int init) {
594    GBDATA *gbd;
595    if (!tree) return init;
596    gbd = GB_entry(tree, key);
597    if (!gbd) return init;
598    return (short)GB_read_byte(gbd);
599}
600
601inline float tree_read_float(GBDATA *tree, const char *key, double init) {
602    GBDATA *gbd;
603    if (!tree) return (float)init;
604    gbd = GB_entry(tree, key);
605    if (!gbd) return (float)init;
606    return (float)GB_read_float(gbd);
607}
608
609
610
611//! moves all node/leaf information from struct GBT_TREE to AP_tree
612void AP_tree::load_node_info() {
613    gr.spread          = tree_read_float(gb_node, "spread",          1.0);
614    gr.left_angle      = tree_read_float(gb_node, "left_angle",      0.0);
615    gr.right_angle     = tree_read_float(gb_node, "right_angle",     0.0);
616    gr.left_linewidth  = tree_read_byte (gb_node, "left_linewidth",  0);
617    gr.right_linewidth = tree_read_byte (gb_node, "right_linewidth", 0);
618    gr.grouped         = tree_read_byte (gb_node, "grouped",         0);
619}
620
621void AP_tree::move_gbt_info(GBT_TREE *tree) {
622    ARB_tree::move_gbt_info(tree);
623    load_node_info();
624}
625
626
627
628#if defined(DEBUG)
629#if defined(DEVEL_RALF)
630#define DEBUG_tree_write_byte
631#endif // DEVEL_RALF
632#endif // DEBUG
633
634
635static GB_ERROR tree_write_byte(GBDATA *gb_tree, AP_tree *node, short i, const char *key, int init) {
636    GBDATA *gbd;
637    GB_ERROR error = 0;
638    if (i==init) {
639        if (node->gb_node) {
640            gbd = GB_entry(node->gb_node, key);
641            if (gbd) {
642#if defined(DEBUG_tree_write_byte)
643                printf("[tree_write_byte] deleting db entry %p\n", gbd);
644#endif // DEBUG_tree_write_byte
645                GB_delete(gbd);
646            }
647        }
648    }
649    else {
650        if (!node->gb_node) {
651            node->gb_node = GB_create_container(gb_tree, "node");
652#if defined(DEBUG_tree_write_byte)
653            printf("[tree_write_byte] created node-container %p\n", node->gb_node);
654#endif // DEBUG_tree_write_byte
655        }
656        gbd = GB_entry(node->gb_node, key);
657        if (!gbd) {
658            gbd = GB_create(node->gb_node, key, GB_BYTE);
659#if defined(DEBUG_tree_write_byte)
660            printf("[tree_write_byte] created db entry %p\n", gbd);
661#endif // DEBUG_tree_write_byte
662        }
663        error = GB_write_byte(gbd, i);
664    }
665    return error;
666}
667
668static GB_ERROR tree_write_float(GBDATA *gb_tree, AP_tree *node, float f, const char *key, float init) {
669    GB_ERROR error = NULL;
670    if (f==init) {
671        if (node->gb_node) {
672            GBDATA *gbd    = GB_entry(node->gb_node, key);
673            if (gbd) error = GB_delete(gbd);
674        }
675    }
676    else {
677        if (!node->gb_node) {
678            node->gb_node             = GB_create_container(gb_tree, "node");
679            if (!node->gb_node) error = GB_await_error();
680        }
681        if (!error) error = GBT_write_float(node->gb_node, key, f);
682    }
683    return error;
684}
685
686GB_ERROR AP_tree::tree_write_tree_rek(GBDATA *gb_tree) {
687    GB_ERROR error = NULL;
688    if (!is_leaf) {
689        error             = get_leftson()->tree_write_tree_rek(gb_tree);
690        if (!error) error = get_rightson()->tree_write_tree_rek(gb_tree);
691
692        if (!error) error = tree_write_float(gb_tree, this, gr.spread,          "spread",          1.0);
693        if (!error) error = tree_write_float(gb_tree, this, gr.left_angle,      "left_angle",      0.0);
694        if (!error) error = tree_write_float(gb_tree, this, gr.right_angle,     "right_angle",     0.0);
695        if (!error) error = tree_write_byte (gb_tree, this, gr.left_linewidth,  "left_linewidth",  0);
696        if (!error) error = tree_write_byte (gb_tree, this, gr.right_linewidth, "right_linewidth", 0);
697        if (!error) error = tree_write_byte (gb_tree, this, gr.grouped,         "grouped",         0);
698    }
699    return error;
700}
701
702GB_ERROR AP_tree_root::saveToDB() {
703    GB_ERROR  error  = GB_push_transaction(get_gb_main());
704    if (get_gb_tree()) {
705        error = get_root_node()->tree_write_tree_rek(get_gb_tree());
706    }
707    else {
708        ap_assert(!gb_tree_gone);      // should have been handled by caller (e.g. in AWT_graphic_tree::save)
709    }
710    if (!error) error = ARB_tree_root::saveToDB();
711    if (!error) update_timers();
712
713    return GB_end_transaction(get_gb_main(), error);
714}
715
716
717GB_ERROR AP_tree::move_group_info(AP_tree *new_group) {
718    GB_ERROR error = 0;
719    if (is_leaf || !name) {
720        error = "Please select a valid group";
721    }
722    else if (!gb_node) {
723        error = "Internal Error: group with name is missing DB-entry";
724    }
725    else if (new_group->is_leaf) {
726        if (new_group->name) {
727            error = GBS_global_string("'%s' is not a valid target for group information of '%s'.", new_group->name, name);
728        }
729        else if (new_group->gb_node) {
730            error = "Internal Error: Target node already has a database entry (but no name)";
731        }
732    }
733
734    if (!error) {
735        if (new_group->name) {
736            if (!new_group->gb_node) {
737                error = "Internal Error: Target node has a database entry (but no name)";
738            }
739            else { // exchange information of two groups
740                GBDATA *tmp_node   = new_group->gb_node;
741                char   *tmp_name   = new_group->name;
742                new_group->gb_node = gb_node;
743                new_group->name    = name;
744                name               = tmp_name;
745                gb_node            = tmp_node;
746            }
747        }
748        else { // move group info
749            new_group->gb_node = this->gb_node;
750            new_group->name    = this->name;
751            this->name         = 0;
752            this->gb_node      = 0;
753        }
754
755        this->load_node_info();
756        new_group->load_node_info();
757
758        {
759            GBDATA *gb_group_name;
760            gb_group_name = GB_entry(new_group->gb_node, "group_name");
761            if (gb_group_name) GB_touch(gb_group_name); // force taxonomy reload
762        }
763    }
764    return error;
765}
766
767void AP_tree::update() {
768    GB_transaction dummy(get_tree_root()->get_gb_main());
769    get_tree_root()->update_timers();
770}
771
772GBT_LEN AP_tree::arb_tree_deep()
773{
774    GBT_LEN l, r;
775    if (is_leaf) return 0.0;
776    l = leftlen + get_leftson()->arb_tree_deep();
777    r = rightlen + get_rightson()->arb_tree_deep();
778    if (l<r) l=r;
779    gr.tree_depth = l;
780    return l;
781}
782
783GBT_LEN AP_tree::arb_tree_min_deep()
784{
785    GBT_LEN l, r;
786    if (is_leaf) return 0.0;
787    l = leftlen + get_leftson()->arb_tree_min_deep();
788    r = rightlen + get_rightson()->arb_tree_min_deep();
789    if (l>r) l=r;
790    gr.min_tree_depth = l;
791    return l;
792}
793
794int AP_tree::arb_tree_set_leafsum_viewsum() // count all visible leafs
795{
796    int l, r;
797    if (is_leaf) {
798        gr.view_sum = 1;
799        gr.leaf_sum = 1;
800        return 1;
801    }
802    l =  get_leftson()->arb_tree_set_leafsum_viewsum();
803    r =  get_rightson()->arb_tree_set_leafsum_viewsum();
804    gr.leaf_sum = r+l;
805    gr.view_sum = get_leftson()->gr.view_sum + get_rightson()->gr.view_sum;
806    if (gr.grouped) {
807        gr.view_sum = (int)pow((double)(gr.leaf_sum - GROUPED_SUM + 9), .33);
808    }
809    return gr.leaf_sum;
810}
811
812int AP_tree::arb_tree_leafsum2()    // count all leafs
813{
814    if (is_leaf) return 1;
815    return get_leftson()->arb_tree_leafsum2() + get_rightson()->arb_tree_leafsum2();
816}
817
818void AP_tree::calc_hidden_flag(int father_is_hidden) {
819    gr.hidden = father_is_hidden;
820    if (!is_leaf) {
821        if (gr.grouped) father_is_hidden = 1;
822        get_leftson()->calc_hidden_flag(father_is_hidden);
823        get_rightson()->calc_hidden_flag(father_is_hidden);
824    }
825}
826
827int AP_tree::calc_color() {
828    int res;
829    if (is_leaf) {
830        if (gb_node) {
831            if (GB_read_flag(gb_node)) {
832                res = AWT_GC_SELECTED;
833            }
834            else {
835                // check for user color
836                long color_group = AWT_species_get_dominant_color(gb_node);
837                if (color_group == 0) {
838                    res = AWT_GC_NSELECTED;
839                }
840                else {
841                    res = AWT_GC_FIRST_COLOR_GROUP+color_group-1;
842                }
843            }
844        }
845        else {
846            res = AWT_GC_SOME_MISMATCHES;
847        }
848    }
849    else {
850        int l = get_leftson()->calc_color();
851        int r = get_rightson()->calc_color();
852
853        if (l == r) res = l;
854
855        else if (l == AWT_GC_SELECTED && r != AWT_GC_SELECTED) res = AWT_GC_UNDIFF;
856        else if (l != AWT_GC_SELECTED && r == AWT_GC_SELECTED) res = AWT_GC_UNDIFF;
857
858        else if (l == AWT_GC_SOME_MISMATCHES) res = r;
859        else if (r == AWT_GC_SOME_MISMATCHES) res = l;
860
861        else if (l == AWT_GC_UNDIFF || r == AWT_GC_UNDIFF) res = AWT_GC_UNDIFF;
862
863        else {
864            ap_assert(l != AWT_GC_SELECTED && r != AWT_GC_SELECTED);
865            ap_assert(l != AWT_GC_UNDIFF && r != AWT_GC_UNDIFF);
866            res = AWT_GC_NSELECTED;
867        }
868    }
869
870    gr.gc = res;
871    if (res == AWT_GC_NSELECTED) {
872        gr.has_marked_children = 0;
873    }
874    else {
875        gr.has_marked_children = 1;
876    }
877    return res;
878}
879
880// Diese Funktion nimmt eine Hashtabelle mit Bakteriennamen und
881// faerbt Bakterien die darin vorkommen mit den entsprechenden Farben
882// in der Hashtabelle ist eine Struktur aus Bak.namen und Farben(GC's)
883int AP_tree::calc_color_probes(GB_HASH *hashptr) {
884    int res;
885
886    if (is_leaf) {
887        if (gb_node) {
888            res = GBS_read_hash(hashptr, name);
889            if (!res && GB_read_flag(gb_node)) { // marked but not in hash -> black
890                res = AWT_GC_BLACK;
891            }
892        }
893        else {
894            res = AWT_GC_SOME_MISMATCHES;
895        }
896    }
897    else {
898        int l = get_leftson()->calc_color_probes(hashptr);
899        int r = get_rightson()->calc_color_probes(hashptr);
900
901        if      (l == r)                      res = l;
902        else if (l == AWT_GC_SOME_MISMATCHES) res = r;
903        else if (r == AWT_GC_SOME_MISMATCHES) res = l;
904        else                                  res = AWT_GC_UNDIFF;
905    }
906    gr.gc = res;
907    return res;
908}
909
910int AP_tree::compute_tree(GBDATA *gb_main)
911{
912    GB_transaction dummy(gb_main);
913    arb_tree_deep();
914    arb_tree_min_deep();
915    arb_tree_set_leafsum_viewsum();
916    calc_color();
917    calc_hidden_flag(0);
918    return 0;
919}
920
921GB_ERROR AP_tree_root::loadFromDB(const char *name) {
922    GB_ERROR error = ARB_tree_root::loadFromDB(name);
923    update_timers();                                // maybe after link() ?
924    return error;
925}
926
927GB_ERROR AP_tree::relink() {
928    GB_transaction dummy(get_tree_root()->get_gb_main()); // open close a transaction
929    GB_ERROR error = GBT_link_tree(get_gbt_tree(), get_tree_root()->get_gb_main(), false, 0, 0); // no status
930    get_tree_root()->update_timers();
931    return error;
932}
933
934AP_UPDATE_FLAGS AP_tree::check_update() {
935    AP_tree_root *troot   = get_tree_root();
936    GBDATA       *gb_main = troot->get_gb_main();
937
938    if (!gb_main) {
939        return AP_UPDATE_RELOADED;
940    }
941    else {
942        GB_transaction dummy(gb_main);
943
944        if (troot->is_tree_updated()) return AP_UPDATE_RELOADED;
945        if (troot->is_species_updated()) return AP_UPDATE_RELINKED;
946        return AP_UPDATE_OK;
947    }
948}
949
950#if defined(WARN_TODO)
951#warning buildLeafList, buildNodeList and buildBranchList should return a AP_tree_list (new class!)
952#endif
953
954void AP_tree::buildLeafList_rek(AP_tree **list, long& num) {
955    // builds a list of all species
956    if (!is_leaf) {
957        get_leftson()->buildLeafList_rek(list, num);
958        get_rightson()->buildLeafList_rek(list, num);
959    }
960    else {
961        list[num] = this;
962        num++;
963    }
964}
965
966void AP_tree::buildLeafList(AP_tree **&list, long &num) {
967    num        = arb_tree_leafsum2();
968    list       = new AP_tree *[num+1];
969    list[num]  = 0;
970    long count = 0;
971
972    buildLeafList_rek(list, count);
973
974    ap_assert(count == num);
975}
976
977void AP_tree::buildNodeList_rek(AP_tree **list, long& num) {
978    // builds a list of all inner nodes (w/o root node)
979    if (!is_leaf) {
980        if (father) list[num++] = this;
981        get_leftson()->buildNodeList_rek(list, num);
982        get_rightson()->buildNodeList_rek(list, num);
983    }
984}
985
986void AP_tree::buildNodeList(AP_tree **&list, long &num) {
987    num = this->arb_tree_leafsum2()-1;
988    list = new AP_tree *[num+1];
989    list[num] = 0;
990    num  = 0;
991    buildNodeList_rek(list, num);
992}
993
994void AP_tree::buildBranchList_rek(AP_tree **list, long& num, bool create_terminal_branches, int deep) {
995    // builds a list of all species
996    // (returns pairs of leafs/father and nodes/father)
997
998    if (deep) {
999        if (father && (create_terminal_branches || !is_leaf)) {
1000            if (father->father) {
1001                list[num++] = this;
1002                list[num++] = get_father();
1003            }
1004            else {                  // root
1005                if (father->leftson == this) {
1006                    list[num++] = this;
1007                    list[num++] = get_brother();
1008                }
1009            }
1010        }
1011        if (!is_leaf) {
1012            get_leftson() ->buildBranchList_rek(list, num, create_terminal_branches, deep-1);
1013            get_rightson()->buildBranchList_rek(list, num, create_terminal_branches, deep-1);
1014        }
1015    }
1016}
1017
1018void AP_tree::buildBranchList(AP_tree **&list, long &num, bool create_terminal_branches, int deep) {
1019    if (deep>=0) {
1020        num = 2;
1021        for (int i=0; i<deep; i++) num *= 2;
1022    }
1023    else {
1024        num = arb_tree_leafsum2() * (create_terminal_branches ? 2 : 1);
1025    }
1026
1027    ap_assert(num >= 0);
1028
1029    list = new AP_tree *[num*2+4];
1030
1031    if (num) {
1032        long count = 0;
1033
1034        buildBranchList_rek(list, count, create_terminal_branches, deep);
1035        list[count] = 0;
1036        num         = count/2;
1037    }
1038}
1039
1040
1041long AP_tree_root::remove_leafs(int awt_remove_type) {
1042    // may remove the complete tree (if awt_remove_type does not contain AWT_REMOVE_BUT_DONT_FREE)
1043
1044    ASSERT_VALID_TREE(get_root_node());
1045
1046    AP_tree **list;
1047    long      count;
1048    get_root_node()->buildLeafList(list, count);
1049
1050    GB_transaction ta(get_gb_main());
1051    long removed = 0;
1052
1053    for (long i=0; i<count; i++) {
1054        bool     removeNode = false;
1055        AP_tree *leaf       = list[i];
1056
1057        if (leaf->gb_node) {
1058            if ((awt_remove_type & AWT_REMOVE_NO_SEQUENCE) && !leaf->get_seq()) {
1059                removeNode = true;
1060            }
1061            else if (awt_remove_type & (AWT_REMOVE_MARKED|AWT_REMOVE_NOT_MARKED)) {
1062                long flag  = GB_read_flag(list[i]->gb_node);
1063                removeNode = (flag && (awt_remove_type&AWT_REMOVE_MARKED)) || (!flag && (awt_remove_type&AWT_REMOVE_NOT_MARKED));
1064            }
1065        }
1066        else {
1067            if (awt_remove_type & AWT_REMOVE_DELETED) {
1068                removeNode = true;
1069            }
1070        }
1071
1072        if (removeNode) {
1073            list[i]->remove();
1074            removed++;
1075            if (!(awt_remove_type & AWT_REMOVE_BUT_DONT_FREE)) {
1076                delete list[i]->father;
1077            }
1078            if (!get_root_node()) {
1079                break; // tree has been deleted
1080            }
1081        }
1082        ASSERT_VALID_TREE(get_root_node());
1083    }
1084    delete [] list;
1085
1086#if defined(CHECK_TREE_STRUCTURE)
1087    if (get_root_node()) ASSERT_VALID_TREE(get_root_node());
1088#endif // CHECK_TREE_STRUCTURE
1089    return removed;
1090}
1091
1092// ----------------------------
1093//      find_innermost_edge
1094
1095class NodeLeafDistance {
1096    AP_FLOAT  downdist, updist;
1097    enum { NLD_NODIST = 0, NLD_DOWNDIST, NLD_BOTHDIST } state;
1098
1099public:
1100
1101    NodeLeafDistance()
1102        : downdist(-1.0)
1103        , updist(-1.0)
1104        , state(NLD_NODIST)
1105    {}
1106
1107    AP_FLOAT get_downdist() const { ap_assert(state >= NLD_DOWNDIST); return downdist; }
1108    void set_downdist(AP_FLOAT DownDist) {
1109        if (state < NLD_DOWNDIST) state = NLD_DOWNDIST;
1110        downdist = DownDist;
1111    }
1112
1113    AP_FLOAT get_updist() const { ap_assert(state >= NLD_BOTHDIST); return updist; }
1114    void set_updist(AP_FLOAT UpDist) {
1115        if (state < NLD_BOTHDIST) state = NLD_BOTHDIST;
1116        updist = UpDist;
1117    }
1118
1119};
1120
1121class EdgeFinder {
1122    map<AP_tree*, NodeLeafDistance> data;
1123
1124    ARB_edge innermost;
1125    AP_FLOAT min_maxDist;
1126
1127    void insert_tree(AP_tree *node) {
1128        if (node->is_leaf) {
1129            data[node].set_downdist(0.0);
1130        }
1131        else {
1132            insert_tree(node->get_leftson());
1133            insert_tree(node->get_rightson());
1134
1135            data[node].set_downdist(max(data[node->get_leftson()].get_downdist()+node->leftlen,
1136                                        data[node->get_rightson()].get_downdist()+node->rightlen));
1137        }
1138    }
1139
1140    void findBetterEdge_sub(AP_tree *node) {
1141        AP_tree *father  = node->get_father();
1142        AP_tree *brother = node->get_brother();
1143
1144        AP_FLOAT len      = node->get_branchlength();
1145        AP_FLOAT brothLen = brother->get_branchlength();
1146
1147        AP_FLOAT upDist   = max(data[father].get_updist(), data[brother].get_downdist()+brothLen);
1148        AP_FLOAT downDist = data[node].get_downdist();
1149
1150        AP_FLOAT edgedist = max(upDist, downDist)+len/2;
1151
1152        if (edgedist<min_maxDist) { // found better edge
1153            innermost   = ARB_edge(node, father);
1154            min_maxDist = edgedist;
1155        }
1156
1157        data[node].set_updist(upDist+len);
1158
1159        if (!node->is_leaf) {
1160            findBetterEdge_sub(node->get_leftson());
1161            findBetterEdge_sub(node->get_rightson());
1162        }
1163    }
1164
1165    void findBetterEdge(AP_tree *node) {
1166        if (!node->is_leaf) {
1167            findBetterEdge_sub(node->get_leftson());
1168            findBetterEdge_sub(node->get_rightson());
1169        }
1170    }
1171
1172public:
1173    EdgeFinder(AP_tree *rootNode)
1174        : innermost(rootNode->get_leftson(), rootNode->get_rightson()) // root-edge
1175    {
1176        insert_tree(rootNode);
1177
1178        AP_tree *lson = rootNode->get_leftson();
1179        AP_tree *rson = rootNode->get_rightson();
1180
1181        AP_FLOAT rootEdgeLen = rootNode->leftlen + rootNode->rightlen;
1182
1183        AP_FLOAT lddist = data[lson].get_downdist();
1184        AP_FLOAT rddist = data[rson].get_downdist();
1185
1186        data[lson].set_updist(rddist+rootEdgeLen);
1187        data[rson].set_updist(lddist+rootEdgeLen);
1188
1189        min_maxDist = max(lddist, rddist)+rootEdgeLen/2;
1190
1191        findBetterEdge(lson);
1192        findBetterEdge(rson);
1193    }
1194
1195    const ARB_edge& innermost_edge() const { return innermost; }
1196};
1197
1198ARB_edge AP_tree_root::find_innermost_edge() {
1199    EdgeFinder edgeFinder(get_root_node());
1200    return edgeFinder.innermost_edge();
1201}
1202
1203// ----------------------------------------
1204
1205void AP_tree::remove_bootstrap() {
1206    freenull(remark_branch);
1207    if (!is_leaf) {
1208        get_leftson()->remove_bootstrap();
1209        get_rightson()->remove_bootstrap();
1210    }
1211}
1212void AP_tree::reset_branchlengths() {
1213    if (!is_leaf) {
1214        leftlen = rightlen = 0.1;
1215
1216        get_leftson()->reset_branchlengths();
1217        get_rightson()->reset_branchlengths();
1218    }
1219}
1220
1221void AP_tree::scale_branchlengths(double factor) {
1222    if (!is_leaf) {
1223        leftlen  *= factor;
1224        rightlen *= factor;
1225
1226        get_leftson()->scale_branchlengths(factor);
1227        get_rightson()->scale_branchlengths(factor);
1228    }
1229}
1230
1231void AP_tree::bootstrap2branchlen() { // copy bootstraps to branchlengths
1232    if (is_leaf) {
1233        set_branchlength(0.1);
1234    }
1235    else {
1236        if (remark_branch && father) {
1237            int    bootstrap = atoi(remark_branch);
1238            double len       = bootstrap/100.0;
1239            set_branchlength(len);
1240        }
1241        get_leftson()->bootstrap2branchlen();
1242        get_rightson()->bootstrap2branchlen();
1243    }
1244}
1245
1246void AP_tree::branchlen2bootstrap() {               // copy branchlengths to bootstraps
1247    freenull(remark_branch);
1248    if (!is_leaf) {
1249        if (!is_root_node()) {
1250            remark_branch = GBS_global_string_copy("%i%%", int(get_branchlength()*100.0 + .5));
1251        }
1252
1253        get_leftson()->branchlen2bootstrap();
1254        get_rightson()->branchlen2bootstrap();
1255    }
1256}
1257
1258
1259AP_tree ** AP_tree::getRandomNodes(int anzahl) {
1260    // function returns a random constructed tree
1261    // root is tree with species (needed to build a list of species)
1262
1263    AP_tree **retlist = NULL;
1264    if (anzahl) {
1265        AP_tree **list; 
1266        long      sumnodes;
1267        buildNodeList(list, sumnodes);
1268
1269        if (sumnodes) {
1270            retlist = (AP_tree **)calloc(anzahl, sizeof(AP_tree *));
1271
1272            long count = sumnodes;
1273            for (int i=0; i< anzahl; i++) {
1274                long num = GB_random(count);
1275
1276                retlist[i] = list[num]; // export node
1277                count--;                // exclude node
1278
1279                list[num]   = list[count];
1280                list[count] = retlist[i];
1281
1282                if (count == 0) count = sumnodes; // restart it
1283            }
1284        }
1285        delete [] list;
1286    }
1287    return retlist;
1288}
1289
1290// --------------------------------------------------------------------------------
1291
1292template <typename T>
1293class ValueCounter {
1294    T   min, max, sum;
1295    int count;
1296
1297    char *mean_min_max_impl() const;
1298    char *mean_min_max_percent_impl() const;
1299
1300    mutable char *buf;
1301    const char *set_buf(char *content) const { freeset(buf, content); return buf; }
1302
1303public:
1304    ValueCounter()
1305        : min(INT_MAX),
1306          max(INT_MIN),
1307          sum(0),
1308          count(0),
1309          buf(NULL)
1310    {}
1311    ValueCounter(const ValueCounter<T>& other)
1312        : min(other.min),
1313          max(other.max),
1314          sum(other.sum),
1315          count(other.count),
1316          buf(NULL)
1317    {}
1318    ~ValueCounter() { free(buf); }
1319
1320    DECLARE_ASSIGNMENT_OPERATOR(ValueCounter<T>);
1321
1322    void count_value(T val) {
1323        count++;
1324        min  = std::min(min, val);
1325        max  = std::max(max, val);
1326        sum += val;
1327    }
1328
1329    int get_count() const { return count; }
1330    T get_min() const { return min; }
1331    T get_max() const { return max; }
1332    double get_mean() const { return sum/double(count); }
1333
1334    const char *mean_min_max()         const { return count ? set_buf(mean_min_max_impl())         : "<not available>"; }
1335    const char *mean_min_max_percent() const { return count ? set_buf(mean_min_max_percent_impl()) : "<not available>"; }
1336
1337    ValueCounter<T>& operator += (const T& inc) {
1338        min += inc;
1339        max += inc;
1340        sum += inc*count;
1341        return *this;
1342    }
1343    ValueCounter<T>& operator += (const ValueCounter<T>& other) {
1344        min    = std::min(min, other.min);
1345        max    = std::max(max, other.max);
1346        sum   += other.sum;
1347        count += other.count;
1348        return *this;
1349    }
1350};
1351
1352template<typename T>
1353inline ValueCounter<T> operator+(const ValueCounter<T>& c1, const ValueCounter<T>& c2) {
1354    return ValueCounter<T>(c1) += c2;
1355}
1356template<typename T>
1357inline ValueCounter<T> operator+(const ValueCounter<T>& c, const T& inc) {
1358    return ValueCounter<T>(c) += inc;
1359}
1360
1361template<> char *ValueCounter<int>::mean_min_max_impl() const {
1362    return GBS_global_string_copy("%.2f (range: %i .. %i)", get_mean(), get_min(), get_max());
1363}
1364template<> char *ValueCounter<double>::mean_min_max_impl() const {
1365    return GBS_global_string_copy("%.2f (range: %.2f .. %.2f)", get_mean(), get_min(), get_max());
1366}
1367template<> char *ValueCounter<double>::mean_min_max_percent_impl() const {
1368    return GBS_global_string_copy("%.2f%% (range: %.2f%% .. %.2f%%)", get_mean()*100.0, get_min()*100.0, get_max()*100.0);
1369}
1370
1371class LongBranchMarker {
1372    double min_rel_diff;
1373    double min_abs_diff;
1374
1375    int furcs;
1376
1377    ValueCounter<double> absdiff;
1378    ValueCounter<double> reldiff;
1379    ValueCounter<double> absdiff_marked;
1380    ValueCounter<double> reldiff_marked;
1381
1382    double perform_marking(AP_tree *at, bool& marked) {
1383        furcs++;
1384
1385        marked = false;
1386        if (at->is_leaf) return 0.0;
1387
1388        bool marked_left;
1389        bool marked_right;
1390
1391        double max = perform_marking(at->get_leftson(),  marked_left)  + at->leftlen;
1392        double min = perform_marking(at->get_rightson(), marked_right) + at->rightlen;
1393
1394        bool max_is_left = true;
1395        if (max<min) {
1396            double h = max; max = min; min = h;
1397            max_is_left = false;
1398        }
1399
1400        double abs_diff = max-min;
1401        absdiff.count_value(abs_diff);
1402
1403        double rel_diff = (max == 0.0) ? 0.0 : abs_diff/max;
1404        reldiff.count_value(rel_diff);
1405
1406        if (abs_diff>min_abs_diff && rel_diff>min_rel_diff) {
1407            if (max_is_left) {
1408                if (!marked_left) {
1409                    at->leftson->mark_subtree();
1410                    marked = true;
1411                }
1412            }
1413            else {
1414                if (!marked_right) {
1415                    at->rightson->mark_subtree();
1416                    marked = true;
1417                }
1418            }
1419        }
1420
1421        if (marked) { // just marked one of my subtrees
1422            absdiff_marked.count_value(abs_diff);
1423            reldiff_marked.count_value(rel_diff);
1424        }
1425        else {
1426            marked = marked_left||marked_right;
1427        }
1428
1429        return min; // use minimal distance for whole subtree
1430    }
1431
1432    static char *meanDiffs(const ValueCounter<double>& abs, const ValueCounter<double>& rel) {
1433        return GBS_global_string_copy(
1434            "Mean absolute diff: %s\n"
1435            "Mean relative diff: %s",
1436            abs.mean_min_max(), 
1437            rel.mean_min_max_percent());
1438    }
1439
1440public:
1441    LongBranchMarker(AP_tree *root, double min_rel_diff_, double min_abs_diff_)
1442        : min_rel_diff(min_rel_diff_),
1443          min_abs_diff(min_abs_diff_),
1444          furcs(0)
1445    {
1446        bool UNUSED;
1447        perform_marking(root, UNUSED);
1448    }
1449
1450    const char *get_report() const {
1451        char *diffs_all    = meanDiffs(absdiff, reldiff);
1452        char *diffs_marked = meanDiffs(absdiff_marked, reldiff_marked);
1453
1454        const char *msg = GBS_global_string(
1455            "Tree contains %i furcations.\n"
1456            "\n"
1457            "%s\n"
1458            "\n"
1459            "%i subtrees have been marked:\n"
1460            "%s\n"
1461            "\n",
1462            furcs-1, // no furcation at root!
1463            diffs_all,
1464            absdiff_marked.get_count(), 
1465            diffs_marked);
1466
1467        free(diffs_all);
1468        free(diffs_marked);
1469
1470        return msg;
1471    }
1472};
1473
1474struct DepthMarker {
1475    // limits (marked if depth and dist are above)
1476    int    min_depth;
1477    double min_rootdist;
1478
1479    // current values (for recursion)
1480    int    depth;
1481    double dist;
1482
1483    // results
1484    ValueCounter<int>    depths,    depths_marked;
1485    ValueCounter<double> distances, distances_marked;
1486
1487    void perform_marking(AP_tree *at, AP_FLOAT atLen) {
1488        int depthInc = atLen == 0.0 ? 0 : 1;  // do NOT increase depth at multifurcations
1489
1490        depth += depthInc;
1491        dist  += atLen;
1492
1493        if (at->is_leaf) {
1494            depths.count_value(depth);
1495            distances.count_value(dist);
1496
1497            int mark = depth >= min_depth && dist >= min_rootdist;
1498            if (at->gb_node) {
1499                GB_write_flag(at->gb_node, mark);
1500                if (mark) {
1501                    depths_marked.count_value(depth);
1502                    distances_marked.count_value(dist);
1503                }
1504            }
1505        }
1506        else {
1507            perform_marking(at->get_leftson(), at->leftlen);
1508            perform_marking(at->get_rightson(), at->rightlen);
1509        }
1510
1511        depth -= depthInc;
1512        dist  -= atLen;
1513    }
1514
1515public:
1516    DepthMarker(AP_tree *root, int min_depth_, double min_rootdist_)
1517        : min_depth(min_depth_),
1518          min_rootdist(min_rootdist_),
1519          depth(0),
1520          dist(0.0)
1521    {
1522        perform_marking(root, 0.0);
1523    }
1524
1525    const char *get_report() const {
1526        int    leafs          = depths.get_count();
1527        int    marked         = depths_marked.get_count();
1528        double balanced_depth = log10(leafs) / log10(2);
1529
1530        const char *msg = GBS_global_string(
1531            "The optimal mean depth of a tree with %i leafs\n"
1532            "      would be %.2f\n"
1533            "\n"
1534            "Your tree:\n"
1535            "mean depth:    %s\n"
1536            "mean distance: %s\n"
1537            "\n"
1538            "%i species (%.2f%%) have been marked:\n"
1539            "mean depth:    %s\n"
1540            "mean distance: %s\n"
1541            ,
1542            leafs,
1543            balanced_depth,
1544            depths.mean_min_max(),
1545            distances.mean_min_max(), 
1546            marked, marked/double(leafs)*100.0,
1547            depths_marked.mean_min_max(), 
1548            distances_marked.mean_min_max()
1549            );
1550        return msg;
1551    }
1552};
1553
1554const char *AP_tree::mark_long_branches(double min_rel_diff, double min_abs_diff) {
1555    // look for asymmetric parts of the tree and mark all species with long branches
1556    return LongBranchMarker(this, min_rel_diff, min_abs_diff).get_report();
1557}
1558const char *AP_tree::mark_deep_leafs(int min_depth, double min_rootdist) {
1559    // mark all leafs with min_depth and min_rootdist
1560    return DepthMarker(this, min_depth, min_rootdist).get_report();
1561}
1562
1563// --------------------------------------------------------------------------------
1564
1565typedef ValueCounter<double> Distance;
1566
1567class DistanceCounter {
1568    ValueCounter<double> min, max, mean;
1569public:
1570
1571    void count_distance(const Distance& d) {
1572        mean.count_value(d.get_mean());
1573        min.count_value(d.get_min());
1574        max.count_value(d.get_max());
1575    }
1576
1577    int get_count() const { return mean.get_count(); }
1578
1579    char *get_report() const {
1580        return GBS_global_string_copy(
1581            "Mean mean distance: %s\n"
1582            "Mean min. distance: %s\n"
1583            "Mean max. distance: %s",
1584            mean.mean_min_max(), 
1585            min.mean_min_max(), 
1586            max.mean_min_max()
1587            );
1588    }
1589};
1590
1591class EdgeDistances {
1592    typedef map<AP_tree*, Distance> DistanceMap;
1593
1594    DistanceMap downdist; // inclusive length of branch itself
1595    DistanceMap updist;   // inclusive length of branch itself
1596
1597    arb_progress progress;
1598
1599    const Distance& calc_downdist(AP_tree *at, AP_FLOAT len) {
1600        if (at->is_leaf) {
1601            Distance d;
1602            d.count_value(len);
1603            downdist[at] = d;
1604
1605            progress.inc();
1606        }
1607        else {
1608            downdist[at] =
1609                calc_downdist(at->get_leftson(), at->leftlen) +
1610                calc_downdist(at->get_rightson(), at->rightlen) +
1611                len;
1612        }
1613        return downdist[at];
1614    }
1615
1616    const Distance& calc_updist(AP_tree *at, AP_FLOAT len) {
1617        ap_assert(at->father); // impossible - root has no updist!
1618
1619        AP_tree *father  = at->get_father();
1620        AP_tree *brother = at->get_brother();
1621
1622        if (father->father) {
1623            ap_assert(updist.find(father) != updist.end());
1624            ap_assert(downdist.find(brother) != downdist.end());
1625
1626            updist[at] = updist[father] + downdist[brother] + len;
1627        }
1628        else {
1629            ap_assert(downdist.find(brother) != downdist.end());
1630
1631            updist[at] = downdist[brother]+len;
1632        }
1633
1634        if (!at->is_leaf) {
1635            calc_updist(at->get_leftson(), at->leftlen);
1636            calc_updist(at->get_rightson(), at->rightlen);
1637        }
1638        else {
1639            progress.inc();
1640        }
1641       
1642        return updist[at];
1643    }
1644
1645    DistanceCounter alldists, markeddists;
1646
1647    void calc_distance_stats(AP_tree *at) {
1648        if (at->is_leaf) {
1649            ap_assert(updist.find(at) != updist.end());
1650
1651            const Distance& upwards = updist[at];
1652
1653            alldists.count_distance(upwards);
1654            if (at->gb_node && GB_read_flag(at->gb_node)) {
1655                markeddists.count_distance(upwards);
1656            }
1657
1658            progress.inc();
1659        }
1660        else {
1661            calc_distance_stats(at->get_leftson());
1662            calc_distance_stats(at->get_rightson());
1663        }
1664    }
1665
1666public:
1667
1668    EdgeDistances(AP_tree *root)
1669        : progress("Analysing distances", root->arb_tree_leafsum2()*3)
1670    {
1671        calc_downdist(root->get_leftson(),  root->leftlen);
1672        calc_downdist(root->get_rightson(), root->rightlen);
1673
1674        calc_updist(root->get_leftson(),  root->leftlen);
1675        calc_updist(root->get_rightson(), root->rightlen);
1676
1677        calc_distance_stats(root);
1678    }
1679
1680    const char *get_report() const {
1681        char *alldists_report    = alldists.get_report();
1682        char *markeddists_report = markeddists.get_report();
1683
1684        const char *msg = GBS_global_string(
1685            "Distance statistic for %i leafs:\n"
1686            "(each leaf to all other leafs)\n"
1687            "\n"
1688            "%s\n"
1689            "\n"
1690            "Distance statistic for %i marked leafs:\n"
1691            "\n"
1692            "%s\n", 
1693            alldists.get_count(), alldists_report,
1694            markeddists.get_count(), markeddists_report);
1695        free(markeddists_report);
1696        free(alldists_report);
1697        return msg;
1698    }
1699};
1700
1701const char *AP_tree::analyse_distances() { return EdgeDistances(this).get_report(); }
1702
1703// --------------------------------------------------------------------------------
1704
1705static int ap_mark_degenerated(AP_tree *at, double degeneration_factor, double& max_degeneration) {
1706    // returns number of species in subtree
1707
1708    if (at->is_leaf) return 1;
1709
1710    int lSons = ap_mark_degenerated(at->get_leftson(), degeneration_factor, max_degeneration);
1711    int rSons = ap_mark_degenerated(at->get_rightson(), degeneration_factor, max_degeneration);
1712
1713    double this_degeneration = 0;
1714
1715    if (lSons<rSons) {
1716        this_degeneration = rSons/double(lSons);
1717        if (this_degeneration >= degeneration_factor) {
1718            at->leftson->mark_subtree();
1719        }
1720
1721    }
1722    else if (rSons<lSons) {
1723        this_degeneration = lSons/double(rSons);
1724        if (this_degeneration >= degeneration_factor) {
1725            at->rightson->mark_subtree();
1726        }
1727    }
1728
1729    if (this_degeneration >= max_degeneration) {
1730        max_degeneration = this_degeneration;
1731    }
1732
1733    return lSons+rSons;
1734}
1735
1736const char *AP_tree::mark_degenerated_branches(double degeneration_factor) {
1737    // marks all species in degenerated branches.
1738    // For all nodes, where one branch contains 'degeneration_factor' more species than the
1739    // other branch, the smaller branch is considered degenerated.
1740
1741    double max_degeneration = 0;
1742    ap_mark_degenerated(this, degeneration_factor, max_degeneration);
1743    return GBS_global_string("Maximum degeneration = %.2f", max_degeneration);
1744}
1745
1746static int ap_mark_duplicates_rek(AP_tree *at, GB_HASH *seen_species) {
1747    if (at->is_leaf) {
1748        if (at->name) {
1749            if (GBS_read_hash(seen_species, at->name)) { // already seen -> mark species
1750                if (at->gb_node) {
1751                    GB_write_flag(at->gb_node, 1);
1752                }
1753                else { // duplicated zombie
1754                    return 1;
1755                }
1756            }
1757            else { // first occurrence
1758                GBS_write_hash(seen_species, at->name, 1);
1759            }
1760        }
1761    }
1762    else {
1763        return
1764            ap_mark_duplicates_rek(at->get_leftson(), seen_species) +
1765            ap_mark_duplicates_rek(at->get_rightson(), seen_species);
1766    }
1767    return 0;
1768}
1769
1770void AP_tree::mark_duplicates() {
1771    GB_HASH *seen_species = GBS_create_hash(gr.leaf_sum, GB_IGNORE_CASE);
1772
1773    int dup_zombies = ap_mark_duplicates_rek(this, seen_species);
1774    if (dup_zombies) {
1775        aw_message(GBS_global_string("Warning: Detected %i duplicated zombies (can't mark them)", dup_zombies));
1776    }
1777    GBS_free_hash(seen_species);
1778}
1779
1780static double ap_just_tree_rek(AP_tree *at) {
1781    if (at->is_leaf) {
1782        return 0.0;
1783    }
1784    else {
1785        double bl = ap_just_tree_rek(at->get_leftson());
1786        double br = ap_just_tree_rek(at->get_rightson());
1787
1788        double l = at->leftlen + at->rightlen;
1789        double diff = fabs(bl - br);
1790        if (l < diff * 1.1) l = diff * 1.1;
1791        double go = (bl + br + l) * .5;
1792        at->leftlen = go - bl;
1793        at->rightlen = go - br;
1794        return go;
1795    }
1796}
1797
1798
1799void AP_tree::justify_branch_lenghs(GBDATA *gb_main) {
1800    // shift branches to create a symmetric looking tree
1801    GB_transaction dummy(gb_main);
1802    ap_just_tree_rek(this);
1803}
1804
1805static void relink_tree_rek(AP_tree *node, void (*relinker)(GBDATA *&ref_gb_node, char *&ref_name, GB_HASH *organism_hash), GB_HASH *organism_hash) {
1806    if (node->is_leaf) {
1807        relinker(node->gb_node, node->name, organism_hash);
1808    }
1809    else {
1810        relink_tree_rek(node->get_leftson(), relinker, organism_hash);
1811        relink_tree_rek(node->get_rightson(), relinker, organism_hash);
1812    }
1813}
1814
1815void AP_tree::relink_tree(GBDATA *gb_main, void (*relinker)(GBDATA *&ref_gb_node, char *&ref_name, GB_HASH *organism_hash), GB_HASH *organism_hash) {
1816    // relinks the tree using a relinker-function
1817    // every node in tree is passed to relinker, relinker might modify
1818    // these values (ref_gb_node and ref_name) and the modified values are written back into tree
1819
1820    GB_transaction  dummy(gb_main);
1821    relink_tree_rek(this, relinker, organism_hash);
1822}
1823
1824void AP_tree::reset_spread() {
1825    gr.spread = 1.0;
1826    if (!is_leaf) {
1827        get_leftson()->reset_spread();
1828        get_rightson()->reset_spread();
1829    }
1830}
1831
1832void AP_tree::reset_rotation() {
1833    if (!is_leaf) {
1834        gr.left_angle  = 0.0;
1835        gr.right_angle = 0.0;
1836        get_leftson()->reset_rotation();
1837        get_rightson()->reset_rotation();
1838    }
1839}
1840
1841void AP_tree::reset_child_linewidths() {
1842    if (!is_leaf) {
1843        gr.left_linewidth  = 0;
1844        gr.right_linewidth = 0;
1845
1846        get_leftson()->reset_child_linewidths();
1847        get_rightson()->reset_child_linewidths();
1848    }
1849}
1850
Note: See TracBrowser for help on using the repository browser.