source: branches/stable/CONSENSUS_TREE/CT_part.hxx

Last change on this file was 11488, checked in by westram, 7 years ago
  • reintegrates 'tree' into 'trunk'
    • implements #417 (multifurcate tree)
    • tree display
      • adds MULTIFURC MODE
      • reordered modes (synchronizes NTREE and PARSIMONY)
    • branch analysis
      • display number of multifurcations in 'mark long branches'
      • display "in-tree-distance" and "per-species-distance"
    • added function to toggle '100%' bootstraps
    • document bug in GBT_remove_leafs (#452)
  • adds:
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.3 KB
Line 
1// ============================================================= //
2//                                                               //
3//   File      : CT_part.hxx                                     //
4//   Purpose   :                                                 //
5//                                                               //
6//   Institute of Microbiology (Technical University Munich)     //
7//   http://www.arb-home.de/                                     //
8//                                                               //
9// ============================================================= //
10
11#ifndef CT_PART_HXX
12#define CT_PART_HXX
13
14#ifndef ARBTOOLS_H
15#include <arbtools.h>
16#endif
17#ifndef ARBDB_BASE_H
18#include <arbdb_base.h>
19#endif
20#ifndef CT_DEF_HXX
21#include "CT_def.hxx"
22#endif
23#ifndef CT_MEM_HXX
24#include "CT_mem.hxx"
25#endif
26#ifndef ARB_ASSERT_H
27#include <arb_assert.h>
28#endif
29
30
31typedef unsigned int PELEM;
32class PART;
33
34#if defined(ASSERTION_USED)
35inline bool is_similar(double d1, double d2, double eps) {
36    double diff = d1-d2;
37    if (diff<0) diff = -diff;
38    return diff<eps;
39}
40#endif
41
42class PartitionSize {
43    PELEM cutmask;  // this mask is used to zero unused bits from the last long (in PART::p)
44    int   longs;    // number of longs per part
45    int   bits;     // number of bits per part ( = overall number of species)
46   
47    mutable size_t id;      // unique id (depending on part creation order, which depends on the added trees)
48
49public:
50    PartitionSize(int len);
51
52    int get_longs() const { return longs; }
53    int get_bits() const { return bits; }
54    size_t get_unique_id() const { id++; return id; }
55    PELEM get_cutmask() const { return cutmask; }
56
57    PART *create_root() const;
58    PELEM *alloc_mem() const { return (PELEM*)getmem(get_longs()*sizeof(PELEM)); }
59};
60
61class PART {
62    const PartitionSize *info;
63
64    PELEM   *p;
65    GBT_LEN  len;    // length between two nodes (weighted by weight)
66    double   weight; // sum of weights
67    size_t   id;
68    int      members;
69
70    PART(const PART& other)
71        : info(other.info),
72          p(info->alloc_mem()),
73          len(other.len),
74          weight(other.weight),
75          id(info->get_unique_id()),
76          members(other.members)
77    {
78        int longs = get_longs();
79        for (int i = 0; i<longs; ++i) {
80            p[i] = other.p[i];
81        }
82        arb_assert(is_valid());
83    }
84    DECLARE_ASSIGNMENT_OPERATOR(PART);
85
86    int count_members() const;
87
88    void set_weight(double pc) {
89        double lenScale  = pc/weight;
90        weight           = pc;
91        len             *= lenScale;
92    }
93
94    static int byte_pos(int pos) { return pos / sizeof(PELEM) / 8; }
95    static PELEM bit_mask(int pos) { return 1 << (pos % (sizeof(PELEM)*8)); }
96
97public:
98
99    PART(const PartitionSize* size_info, double weight_)
100        : info(size_info),
101          p(info->alloc_mem()),
102          len(0),
103          weight(weight_),
104          id(info->get_unique_id()),
105          members(0)
106    {
107        arb_assert(is_valid());
108    }
109    ~PART() {
110        free(p);
111    }
112
113    PART *clone() const { return new PART(*this); }
114
115    int get_longs() const { return info->get_longs(); }
116    int get_maxsize() const { return info->get_bits(); }
117    PELEM get_cutmask() const { return info->get_cutmask(); }
118
119    bool is_valid() const {
120        PELEM last = p[get_longs()-1];
121
122        bool only_valid_bits    = (last&get_cutmask()) == last;
123        bool valid_member_count = members >= 0 && members <= get_maxsize();
124        bool valid_weight       = weight >= 0.0;
125
126        return only_valid_bits && valid_member_count && valid_weight;
127    }
128
129    void setbit(int pos) {
130        /*! set the bit at position 'pos'
131         */
132        arb_assert(is_valid());
133
134        int   idx = byte_pos(pos);
135        PELEM bit = bit_mask(pos);
136
137        if (!(p[idx]&bit)) {
138            p[idx] |= bit;
139            members++;
140        }
141
142        arb_assert(is_valid());
143    }
144    bool bit_is_set(int pos) const {
145        /*! return true if the bit at the position 'pos' is set
146         */
147        arb_assert(is_valid());
148
149        int   idx = byte_pos(pos);
150        PELEM bit = bit_mask(pos);
151
152        return p[idx] & bit;
153    }
154
155    void set_len(GBT_LEN length) {
156        arb_assert(is_valid());
157        len = length*weight;
158    }
159    GBT_LEN get_len() const {
160        if (weight == 0.0) {
161            // if weight is 0.0 = > branch never occurred!
162            return 1.0; // return worst distance between these nodes (which are disjunct trees)
163        }
164        return len/weight;
165    }
166
167    double get_weight() const { return weight; }
168
169    void addWeightAndLength(const PART *other) {
170        weight += other->weight;
171        len    += other->len;
172    }
173
174    void takeMean(double overall_weight) {
175        set_weight(get_weight() / overall_weight);
176        arb_assert(is_valid());
177        arb_assert(weight <= 1.0);
178    }
179
180    void set_faked_weight(double w) { set_weight(w); }
181
182    void invert();
183    void invertInSuperset(const PART *superset);
184
185    bool is_standardized() const;
186    void standardize();
187
188    void add_members_from(const PART *source);
189
190    bool is_subset_of(const PART *other) const {
191        int longs = get_longs();
192        for (int i=0; i<longs; i++) {
193            if ((p[i] & other->p[i]) != p[i]) {
194                return false;
195            }
196        }
197        return true;
198    }
199    bool is_real_son_of(const PART *father) const { return is_subset_of(father) && differs(father); }
200
201    bool overlaps_with(const PART *other) const;
202    bool disjunct_from(const PART *other) const { return !overlaps_with(other); }
203
204    bool equals(const PART *other) const;
205    bool differs(const PART *other) const { return !equals(other); }
206
207    int index() const;
208    unsigned key() const;
209
210    int get_members() const { return members; }
211    int get_nonmembers() const { return get_maxsize() - get_members(); }
212
213    bool is_leaf_edge() const { int size = members; return size == 1 || size == (get_maxsize()-1); }
214    int distance_to_tree_center() const { return abs(get_maxsize()/2 - members); }
215
216    int insertionOrder_cmp(const PART *other) const;
217    int topological_cmp(const PART *other) const;
218
219#if defined(NTREE_DEBUG_FUNCTIONS)
220    void print() const;
221    static void start_pretty_printing(const class CharPtrArray& names_);
222    static void stop_pretty_printing();
223#endif
224};
225
226#else
227#error CT_part.hxx included twice
228#endif // CT_PART_HXX
Note: See TracBrowser for help on using the repository browser.