source: branches/port5/AWT/AWT_tree.cxx

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