← Back to branch summary

~ubuntu-branches/debian/jessie/arb/jessie

« back to all changes in this revision

Viewing changes to STAT/st_ml.hxx

  • Committer: Package Import Robot
  • Author(s): Elmar Pruesse, Andreas Tille, Elmar Pruesse
  • Date: 2014-09-02 15:15:06 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20140902151506-jihq58b3iz342wif
Tags: 6.0.2-1
http://bugs.debian.org/741890
http://bugs.debian.org/757497
[ Andreas Tille ]
* New upstream version
  Closes: #741890
* debian/upstream -> debian/upstream/metadata
* debian/control:
   - Build-Depends: added libglib2.0-dev
   - Depends: added mafft, mrbayes
* debian/rules
   - Add explicite --remove-section=.comment option to manual strip call
* cme fix dpkg-control
* arb-common.dirs: Do not create unneeded lintian dir
* Add turkish debconf translation (thanks for the patch to Mert Dirik
  <mertdirik@gmail.com>)
  Closes: #757497

[ Elmar Pruesse ]
* patches removed:
   - 10_config.makefiles.patch,
     80_no_GL.patch
       removed in favor of creating file from config.makefile.template via 
       sed in debian/control
   - 20_Makefile_main.patch
       merged upstream
   - 21_Makefiles.patch
       no longer needed
   - 30_tmpfile_CVE-2008-5378.patch: 
       merged upstream
   - 50_fix_gcc-4.8.patch:
       merged upstream
   - 40_add_libGLU.patch:
       libGLU not needed for arb_ntree)
   - 60_use_debian_packaged_raxml.patch:
       merged upstream
   - 70_hardening.patch
       merged upstream
   - 72_add_math_lib_to_linker.patch
       does not appear to be needed
* patches added:
   - 10_upstream_r12793__show_db_load_progress:
       backported patch showing progress while ARB is loading a database
       (needed as indicator/splash screen while ARB is launching)
   - 20_upstream_r12794__socket_permissions:
       backported security fix
   - 30_upstream_r12814__desktop_keywords:
       backported add keywords to desktop (fixes lintian warning)
   - 40_upstream_r12815__lintian_spelling:
       backported fix for lintian reported spelling errors
   - 50_private_nameservers
       change configuration to put nameservers into users home dirs
       (avoids need for shared writeable directory)
   - 60_use_debian_phyml
       use phyml from debian package for both interfaces in ARB
* debian/rules:
   - create config.makefile from override_dh_configure target
   - use "make tarfile" in override_dh_install
   - remove extra cleaning not needed for ARB 6
   - use "dh_install --list-missing" to avoid missing files
   - added override_dh_fixperms target
* debian/control:
   - added libarb-dev package
   - Depends: added phyml, xdg-utils
   - Suggests: removed phyml
   - fix lintian duplicate-short-description (new descriptions)
* debian/*.install:
   - "unrolled" confusing globbing to select files
   - pick files from debian/tmp
   - moved all config files to /etc/arb
* debian/arb-common.templates: updated
* scripts:
   - removed arb-add-pt-server
   - launch-wrapper: 
     - only add demo.arb to newly created $ARBUSERDATA
     - pass commandline arguments through bin/arb wrapper
   - preinst: removing old PT server index files on upgrade from 5.5*
   - postinst: set setgid on shared PT dir
* rewrote arb.1 manfile
* added file icon for ARB databases
* using upstream arb_tcp.dat

Show diffs side-by-side

added added

removed removed

Lines of Context:
STAT/st_ml.hxx
1
 
#ifndef ARB_ASSERT_H
2
 
#include <arb_assert.h>
3
 
#endif
4
 
#define st_assert(bed) arb_assert(bed)
5
 
 
6
 
enum AWT_dna_base {
 
1
// =============================================================== //
 
2
//                                                                 //
 
3
//   File      : st_ml.hxx                                         //
 
4
//   Purpose   :                                                   //
 
5
//                                                                 //
 
6
//   Institute of Microbiology (Technical University Munich)       //
 
7
//   http://www.arb-home.de/                                       //
 
8
//                                                                 //
 
9
// =============================================================== //
 
10
 
 
11
#ifndef ST_ML_HXX
 
12
#define ST_ML_HXX
 
13
 
 
14
#ifndef ARBDB_BASE_H
 
15
#include <arbdb_base.h>
 
16
#endif
 
17
#ifndef _GLIBCXX_ALGORITHM
 
18
#include <algorithm>
 
19
#endif
 
20
#ifndef ATTRIBUTES_H
 
21
#include <attributes.h>
 
22
#endif
 
23
#ifndef ARBTOOLS_H
 
24
#include <arbtools.h>
 
25
#endif
 
26
#ifndef CB_H
 
27
#include <cb.h>
 
28
#endif
 
29
 
 
30
#define st_assert(cond) arb_assert(cond)
 
31
 
 
32
class ColumnStat;
 
33
class MostLikelySeq;
 
34
 
 
35
class AW_root;
 
36
class AW_awar;
 
37
class AW_window;
 
38
 
 
39
class AP_tree;
 
40
class AP_tree_root;
 
41
class GBT_TREE;
 
42
class AP_filter;
 
43
class WeightedFilter;
 
44
 
 
45
typedef unsigned char ST_ML_Color;
 
46
 
 
47
typedef float ST_FLOAT;
 
48
// typedef double ST_FLOAT; // careful, using double has quite an impact on the memory footprint
 
49
 
 
50
enum DNA_Base {
7
51
    ST_A,
8
52
    ST_C,
9
53
    ST_G,
10
54
    ST_T,
11
55
    ST_GAP,
12
 
    ST_MAX_BASE,
 
56
    ST_MAX_BASE,                                    // [can't be changed w/o recoding several parts]
13
57
    ST_UNKNOWN = -1
14
58
};
15
59
 
16
 
extern class AWT_dna_table {
17
 
    char char_to_enum_table[256];
18
 
public:
19
 
    AWT_dna_base char_to_enum(char i) {
20
 
        return (AWT_dna_base)char_to_enum_table[(unsigned char)i];
21
 
    }
22
 
    AWT_dna_table();
23
 
} awt_dna_table;
24
 
 
25
 
typedef unsigned char ST_ML_Color;
26
 
 
27
 
const int ST_MAX_SEQ_PART = 256;
28
 
// should be greater than the editor width
29
 
// otherwise extrem performance penalties
30
 
const int ST_BUCKET_SIZE = 16;
31
 
// at minumum ST_BUCKET_SIZE characters are calculated per call
32
 
const int LD_BUCKET_SIZE = 4; // log dualis of ST_BUCKET_SIZE
33
 
 
34
 
class ST_base_vector {
35
 
public:
36
 
    float b[ST_MAX_BASE]; // acgt-
37
 
    int ld_lik;
38
 
    float lik; // likelihood  = 2^ld_lik * lik * (b[0] + b[1] + b[2] ..)
39
 
    void set(char base, ST_base_vector *frequencies);
40
 
    inline void mult(ST_base_vector *other);
 
60
struct ST_base_vector {
 
61
    ST_FLOAT b[ST_MAX_BASE];                        // acgt-
 
62
    int      ld_lik;
 
63
    ST_FLOAT lik;                                   // likelihood = 2^ld_lik * lik * (b[0] + b[1] + b[2] ..)
 
64
 
 
65
    void setTo(ST_FLOAT freq) {
 
66
        for (int base = ST_A; base<ST_MAX_BASE; ++base) {
 
67
            b[base] = freq;
 
68
        }
 
69
    }
 
70
    void multiplyWith(ST_FLOAT factor) {
 
71
        for (int base = ST_A; base<ST_MAX_BASE; ++base) {
 
72
            b[base] *= factor;
 
73
        }
 
74
    }
 
75
    void increaseBy(ST_FLOAT inc) {
 
76
        for (int base = ST_A; base<ST_MAX_BASE; ++base) {
 
77
            b[base] += inc;
 
78
        }
 
79
    }
 
80
    void makeInverseOf(const ST_base_vector& other, ST_FLOAT other_min_freq) {
 
81
        for (int base = ST_A; base<ST_MAX_BASE; ++base) {
 
82
            b[base] = other_min_freq / other.b[base];
 
83
        }
 
84
    }
 
85
 
 
86
    void setBase(const ST_base_vector& frequencies, char base);
 
87
 
41
88
    void check_overflow();
42
89
    void print();
 
90
 
 
91
    ST_FLOAT summarize() const { return b[ST_A] + b[ST_C] + b[ST_G] + b[ST_T] + b[ST_GAP]; }
 
92
    ST_FLOAT min_frequency() { return std::min(std::min(std::min(b[ST_A], b[ST_C]), std::min(b[ST_G], b[ST_T])), b[ST_GAP]); }
 
93
    ST_FLOAT max_frequency() { return std::max(std::max(std::max(b[ST_A], b[ST_C]), std::max(b[ST_G], b[ST_T])), b[ST_GAP]); }
 
94
 
 
95
    inline ST_base_vector& operator *= (const ST_base_vector& other);
43
96
};
44
97
 
45
98
class ST_rate_matrix {
46
 
    float m[ST_MAX_BASE][ST_MAX_BASE];
 
99
    // symmetric matrix containing only two different values
 
100
 
 
101
    ST_FLOAT diag;                                  // value for indices [0][0], [1][1]...
 
102
    ST_FLOAT rest;                                  // other indices
 
103
 
47
104
public:
 
105
 
 
106
    ST_rate_matrix() : diag(0.0), rest(0.0) {}
48
107
    void set(double dist, double TT_ratio);
49
 
    inline void mult(ST_base_vector *in, ST_base_vector *out);
50
 
    void print();
51
 
};
52
 
 
53
 
class ST_ML;
54
 
class AWT_csp;
55
 
 
56
 
/** Note: Because we have only limited memory we split the
57
 
 sequence into ST_MAX_SEQ_PART long parts */
58
 
class ST_sequence_ml : private AP_sequence {
59
 
    friend class ST_ML;
60
 
public:
61
 
 
62
 
    GBDATA *gb_data; // the sequence
63
 
    static ST_base_vector *tmp_out; // len = alignment length
64
 
 
65
 
protected:
66
 
 
67
 
    ST_ML *st_ml; // link to a global ST object
68
 
    ST_base_vector *sequence; // A part of the sequence
69
 
    int last_updated;
70
 
    ST_ML_Color *color_out;
71
 
    int *color_out_valid_till; // color_out is valid up to
72
 
 
73
 
public:
74
 
    void delete_sequence(); // remove link to database
75
 
    void sequence_change(); // sequence has changed in db
76
 
    AP_FLOAT combine(const AP_sequence* lefts, const AP_sequence *rights);
77
 
    void partial_match(const AP_sequence* part, long *overlap, long *penalty) const;
78
 
    ST_sequence_ml(AP_tree_root *rooti, ST_ML *st_ml);
79
 
    ~ST_sequence_ml();
80
 
    AP_sequence *dup(void);
81
 
 
82
 
    void set(const char *sequence);
83
 
    void set_gb(GBDATA *gbd);
84
 
 
85
 
    void set_sequence(); // start at st_ml->base
86
 
 
87
 
    void go(const ST_sequence_ml *lefts, double leftl,
88
 
            const ST_sequence_ml *rights, double rightl);
89
 
    void ungo(); // undo go
90
 
 
91
 
    void calc_out(ST_sequence_ml *sequence_of_brother, double dist);
92
 
    void print();
93
 
};
94
 
 
95
 
class AW_window;
96
 
typedef void (*AW_CB0)(AW_window*);
97
 
 
98
 
class ST_ML {
99
 
    char *alignment_name;
100
 
    friend AP_tree *st_ml_convert_species_name_to_node(ST_ML *st_ml,
101
 
            const char *species_name);
102
 
    GB_HASH *hash_2_ap_tree; // hash table to get from name to tree_node
103
 
    GB_HASH *keep_species_hash; // temporary hash to find
104
 
    int refresh_n;
105
 
    int *not_valid; // which columns are valid
106
 
 
107
 
    ST_sequence_ml *do_tree(AP_tree *node);
108
 
    void undo_tree(AP_tree *node); //opposite of do_tree
 
108
    inline void transform(const ST_base_vector& in, ST_base_vector& out) const;
 
109
    void print();
 
110
};
 
111
 
 
112
class ST_ML : virtual Noncopyable {
 
113
 
 
114
    /* Note: Because we have only limited memory we split the
 
115
     * sequence into several ST_MAX_SEQ_PART-sized parts during stat-calculation
 
116
     */
 
117
 
 
118
    char    *alignment_name;
 
119
    GB_HASH *hash_2_ap_tree;                        // hash table to get from name to tree_node
 
120
    GB_HASH *keep_species_hash;                     // temporary hash to find ??? 
 
121
    int      refresh_n;
 
122
    int     *not_valid;                             // which columns are valid
 
123
 
 
124
    AP_tree_root *tree_root;
 
125
    int           latest_modification;              // last mod
 
126
    size_t        first_pos;                        // first..
 
127
    size_t        last_pos;                         // .. and last alignment position of current slice
 
128
 
 
129
    AW_CB0     postcalc_cb;
 
130
    AW_window *cb_window;
 
131
 
 
132
    GBDATA *gb_main;
 
133
 
 
134
    ColumnStat *column_stat;
 
135
    // if column_stat == 0 -> rates and ttratio are owned by ST_ML,
 
136
    // otherwise they belong to column_stat
 
137
    const float *rates;                             // column independent
 
138
    const float *ttratio;                           // column independent
 
139
 
 
140
    ST_base_vector *base_frequencies;               // column independent
 
141
    ST_base_vector *inv_base_frequencies;           // column independent
 
142
 
 
143
    double max_dist;                                // max_dist for rate_matrices
 
144
    double step_size;                               // max_dist/step_size matrices
 
145
 
 
146
    int             max_rate_matrices;              // number of rate_matrices
 
147
    ST_rate_matrix *rate_matrices;                  // one matrix for each distance
 
148
 
 
149
    bool is_initialized;
 
150
 
 
151
    size_t distance_2_matrices_index(int distance) {
 
152
        if (distance<0) distance = -distance;       // same matrix for neg/pos distance
 
153
 
 
154
        return distance >= max_rate_matrices
 
155
            ? (max_rate_matrices-1)
 
156
            : size_t(distance);
 
157
    }
 
158
 
 
159
    MostLikelySeq *getOrCreate_seq(AP_tree *node);
 
160
 
 
161
    const MostLikelySeq *get_mostlikely_sequence(AP_tree *node);
 
162
    void                 undo_tree(AP_tree *node);  // opposite of get_mostlikely_sequence()
 
163
 
109
164
    void insert_tree_into_hash_rek(AP_tree *node);
110
 
    void create_matrizes(double max_disti, int nmatrizes);
 
165
    void create_matrices(double max_disti, int nmatrices);
111
166
    void create_frequencies();
 
167
    
112
168
    static long delete_species(const char *key, long val, void *cd_st_ml);
 
169
 
113
170
public:
114
 
    AP_tree_root *tree_root;
115
 
    int latest_modification; // last mod;
116
 
    int base;
117
 
    int to;
118
 
    AW_CB0 refresh_func;
119
 
    AW_window *aw_window;
120
 
 
121
 
    GBDATA *gb_main;
122
 
    float *ttratio; // column independent
123
 
    ST_base_vector *base_frequencies; // column independent
124
 
    ST_base_vector *inv_base_frequencies; // column independent
125
 
    float *rates; // column independent
126
 
    double max_dist; // max_dist for rate_matrizes
127
 
    double step_size; // max_dist/step_size matrizes
128
 
    int max_matr;
129
 
    ST_rate_matrix *rate_matrizes; // for each distance a new matrix
130
 
    long alignment_len;
131
 
    AWT_csp *awt_csp;
 
171
 
 
172
    ST_ML(GBDATA *gb_main);
 
173
    ~ST_ML();
 
174
 
 
175
    void create_column_statistic(AW_root *awr, const char *awarname, AW_awar *awar_default_alignment);
 
176
    ColumnStat *get_column_statistic() { return column_stat; }
 
177
 
 
178
    GB_ERROR calc_st_ml(const char           *tree_name,
 
179
                        const char           *alignment_name,
 
180
                        const char           *species_names, // 0 -> all [marked] species (else species_names is a (char)1 separated list of species)
 
181
                        int                   marked_only,
 
182
                        ColumnStat           *colstat,
 
183
                        const WeightedFilter *weighted_filter) __ATTR__USERESULT;
 
184
 
 
185
    void cleanup();
 
186
 
 
187
    const AP_filter *get_filter() const;
 
188
    size_t get_filtered_length() const;
 
189
    size_t get_alignment_length() const;
 
190
 
 
191
    ST_rate_matrix& get_matrix_for(int distance) {
 
192
        return rate_matrices[distance_2_matrices_index(distance)];
 
193
    }
 
194
 
 
195
    GBDATA *get_gb_main() const { return gb_main; }
 
196
    const GBT_TREE *get_gbt_tree() const;
 
197
    size_t count_species_in_tree() const;
 
198
 
 
199
    AP_tree *find_node_by_name(const char *species_name);
 
200
 
 
201
    void set_postcalc_callback(AW_CB0 postcalc_cb_, AW_window *cb_window_) {
 
202
        postcalc_cb = postcalc_cb_;
 
203
        cb_window   = cb_window_;
 
204
    }
 
205
    void do_postcalc_callback() { if (postcalc_cb) postcalc_cb(cb_window); }
 
206
 
 
207
    size_t get_first_pos() const { return first_pos; }
 
208
    size_t get_last_pos() const { return last_pos; }
 
209
 
 
210
    double get_step_size() const { return step_size; }
 
211
 
 
212
    const ST_base_vector *get_inv_base_frequencies() const { return inv_base_frequencies; }
 
213
    const ST_base_vector& get_base_frequency_at(size_t pos) const { return base_frequencies[pos]; }
 
214
    float get_rate_at(size_t pos) const { return rates[pos]; }
 
215
 
132
216
    void set_modified(int *what = 0);
133
 
    void set_refresh(); // set flag for refresh
134
 
 
135
 
    ~ST_ML();
136
 
    ST_ML(GBDATA *gb_main);
137
 
    void print();
138
 
    int is_inited;
139
 
 
140
 
    GB_ERROR init(const char *tree_name, const char *alignment_name,
141
 
            const char *species_names, int marked_only,
142
 
            const char *filter_string, AWT_csp *awt_csp);
143
 
    // species_names is 0 -> all [marked] species (else species_names is a (char)1 seperated list of species)
144
 
    // filter_string==0 -> no filter
145
 
 
146
 
    void clear_all(); // delete all caches
147
 
 
148
 
 
149
 
    ST_sequence_ml *get_ml_vectors(char *species_name, AP_tree *node,
150
 
            int start_ali_pos, int end_ali_pos);
151
 
    ST_ML_Color *get_color_string(char *species_name, AP_tree *node,
152
 
            int start_ali_pos, int end_ali_pos);
153
 
 
154
 
    int update_ml_likelihood(char *result[4], int *latest_update,
155
 
            char *species_name, AP_tree *node);
 
217
    void set_refresh();                             // set flag for refresh
 
218
 
 
219
    void clear_all();                               // delete all caches
 
220
 
 
221
    MostLikelySeq *get_ml_vectors(const char *species_name, AP_tree *node, size_t start_ali_pos, size_t end_ali_pos);
 
222
    ST_ML_Color *get_color_string(const char *species_name, AP_tree *node, size_t start_ali_pos, size_t end_ali_pos);
 
223
 
 
224
    bool update_ml_likelihood(char *result[4], int& latest_update, const char *species_name, AP_tree *node);
156
225
 
157
226
    int refresh_needed();
158
227
};
 
228
 
 
229
#else
 
230
#error st_ml.hxx included twice
 
231
#endif // ST_ML_HXX