source: branches/lib/DIST/DI_matr.cxx

Last change on this file was 19432, checked in by westram, 2 years ago
  • GBT_get_alignment_len
    • now also reports error if alignment length is zero
      • this case often was unhandled and did easily lead to allocation problems.
    • catch error in case of zero alignment length ⇒ fixes alloc-size-larger-than-warning (in NT_count_different_chars).
    • check + fix callers.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 68.8 KB
Line 
1// =============================================================== //
2//                                                                 //
3//   File      : DI_matr.cxx                                       //
4//   Purpose   :                                                   //
5//                                                                 //
6//   Institute of Microbiology (Technical University Munich)       //
7//   http://www.arb-home.de/                                       //
8//                                                                 //
9// =============================================================== //
10
11#include "di_protdist.hxx"
12#include "di_clusters.hxx"
13#include "di_view_matrix.hxx"
14#include "di_awars.hxx"
15
16#include <neighbourjoin.hxx>
17#include <AP_filter.hxx>
18#include <CT_ctree.hxx>
19#include <ColumnStat.hxx>
20
21#include <awt.hxx>
22#include <awt_sel_boxes.hxx>
23#include <awt_filter.hxx>
24
25#include <aw_preset.hxx>
26#include <aw_awars.hxx>
27#include <aw_file.hxx>
28#include <aw_msg.hxx>
29#include <aw_root.hxx>
30
31#include <gui_aliview.hxx>
32
33#include <climits>
34#include <ctime>
35#include <cmath>
36#include <arb_sort.h>
37#include <arb_global_defs.h>
38#include <macros.hxx>
39#include <ad_cb.h>
40#include <awt_TreeAwars.hxx>
41#include <arb_defs.h>
42
43using std::string;
44
45// --------------------------------------------------------------------------------
46
47#define AWAR_DIST_BOOTSTRAP_COUNT    AWAR_DIST_PREFIX "bootstrap/count"
48#define AWAR_DIST_CANCEL_CHARS       AWAR_DIST_PREFIX "cancel/chars"
49#define AWAR_DIST_MATRIX_AUTO_RECALC AWAR_DIST_PREFIX "recalc"
50
51#define AWAR_DIST_COLUMN_STAT_NAME AWAR_DIST_COLUMN_STAT_PREFIX "name"
52
53#define AWAR_DIST_TREE_SORT_NAME        "tmp/" AWAR_DIST_TREE_PREFIX "sort_tree_name"
54#define AWAR_DIST_TREE_COMP_NAME        "tmp/" AWAR_DIST_TREE_PREFIX "compr_tree_name"
55#define AWAR_DIST_TREE_STD_NAME         AWAR_DIST_TREE_PREFIX "tree_name"
56#define AWAR_DIST_MATRIX_AUTO_CALC_TREE AWAR_DIST_TREE_PREFIX "autocalc"
57
58#define AWAR_DIST_SAVE_MATRIX_TYPE     AWAR_DIST_SAVE_MATRIX_BASE "/type"
59#define AWAR_DIST_SAVE_MATRIX_FILENAME AWAR_DIST_SAVE_MATRIX_BASE "/file_name"
60
61// --------------------------------------------------------------------------------
62
63DI_GLOBAL_MATRIX GLOBAL_MATRIX;
64
65static MatrixDisplay     *matrixDisplay = NULp;
66static AP_userdef_matrix  userdef_DNA_matrix(AP_MAX, AWAR_DIST_MATRIX_DNA_BASE);
67
68static AP_matrix *get_user_matrix() {
69    AW_root *awr = AW_root::SINGLETON;
70    if (!awr->awar(AWAR_DIST_MATRIX_DNA_ENABLED)->read_int()) {
71        return NULp;
72    }
73    userdef_DNA_matrix.update_from_awars(awr);
74    userdef_DNA_matrix.normize();
75    return &userdef_DNA_matrix;
76}
77
78class BoundWindowCallback : virtual Noncopyable {
79    AW_window      *aww;
80    WindowCallback  cb;
81public:
82    BoundWindowCallback(AW_window *aww_, const WindowCallback& cb_)
83        : aww(aww_),
84          cb(cb_)
85    {}
86    void operator()() { cb(aww); }
87};
88
89static SmartPtr<BoundWindowCallback> recalculate_matrix_cb;
90static SmartPtr<BoundWindowCallback> recalculate_tree_cb;
91
92static GB_ERROR last_matrix_calculation_error = NULp;
93
94static void matrix_changed_cb() {
95    if (matrixDisplay) {
96        matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
97        matrixDisplay->update_display();
98    }
99}
100
101struct RecalcNeeded {
102    bool matrix;
103    bool tree;
104
105    RecalcNeeded() : matrix(false), tree(false) { }
106};
107
108static RecalcNeeded need_recalc;
109
110CONSTEXPR unsigned UPDATE_DELAY = 200;
111
112static unsigned update_cb(AW_root *aw_root);
113inline void add_update_cb() {
114    AW_root::SINGLETON->add_timed_callback(UPDATE_DELAY, makeTimedCallback(update_cb));
115}
116
117inline void matrix_needs_recalc_cb() {
118    need_recalc.matrix = true;
119    add_update_cb();
120}
121inline void tree_needs_recalc_cb() {
122    need_recalc.tree = true;
123    add_update_cb();
124}
125
126inline void compressed_matrix_needs_recalc_cb() {
127    if (GLOBAL_MATRIX.has_type(DI_MATRIX_COMPRESSED)) {
128        matrix_needs_recalc_cb();
129    }
130}
131
132static unsigned update_cb(AW_root *aw_root) {
133    if (need_recalc.matrix) {
134        GLOBAL_MATRIX.forget();
135        need_recalc.matrix = false; // because it's forgotten
136
137        int matrix_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->read_int();
138        if (matrix_autocalc) {
139            bool recalc_now    = true;
140            int  tree_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
141            if (!tree_autocalc) recalc_now = matrixDisplay ? matrixDisplay->willShow() : false;
142
143            if (recalc_now) {
144                di_assert(recalculate_matrix_cb.isSet());
145                (*recalculate_matrix_cb)();
146                di_assert(need_recalc.tree == true);
147            }
148        }
149        di_assert(need_recalc.matrix == false);
150    }
151
152    if (need_recalc.tree) {
153        int tree_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
154        if (tree_autocalc) {
155            di_assert(recalculate_tree_cb.isSet());
156            (*recalculate_tree_cb)();
157            need_recalc.matrix = false; // otherwise endless loop, e.g. if output-tree is used for sorting
158        }
159    }
160
161    return 0; // do not call again
162}
163
164static void auto_calc_changed_cb(AW_root *aw_root) {
165    int matrix_autocalc = aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->read_int();
166    int tree_autocalc   = aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->read_int();
167
168    if (matrix_autocalc && !GLOBAL_MATRIX.exists()) matrix_needs_recalc_cb();
169    if (tree_autocalc && (matrix_autocalc || GLOBAL_MATRIX.exists())) tree_needs_recalc_cb();
170}
171
172static AW_window *create_dna_matrix_window(AW_root *aw_root) {
173    AW_window_simple *aws = new AW_window_simple;
174    aws->init(aw_root, "SET_DNA_MATRIX", "SET MATRIX");
175    aws->auto_increment(50, 50);
176    aws->button_length(10);
177    aws->callback(AW_POPDOWN);
178    aws->create_button("CLOSE", "CLOSE");
179
180    aws->callback(makeHelpCallback("user_matrix.hlp"));
181    aws->create_button("HELP", "HELP");
182
183    aws->at_newline();
184
185    userdef_DNA_matrix.create_input_fields(aws);
186    aws->window_fit();
187    return aws;
188}
189
190static void selected_tree_changed_cb() {
191    if (AW_root::SINGLETON->awar(AWAR_DIST_CORR_TRANS)->read_int() == DI_TRANSFORMATION_FROM_TREE) {
192        matrix_needs_recalc_cb();
193    }
194}
195
196GB_ERROR DI_create_matrix_variables(AW_root *aw_root, AW_default def, AW_default db) {
197    GB_transaction ta(db);
198
199    userdef_DNA_matrix.set_descriptions(AP_A,   "A");
200    userdef_DNA_matrix.set_descriptions(AP_C,   "C");
201    userdef_DNA_matrix.set_descriptions(AP_G,   "G");
202    userdef_DNA_matrix.set_descriptions(AP_T,   "TU");
203    userdef_DNA_matrix.set_descriptions(AP_GAP, "GAP");
204
205    userdef_DNA_matrix.create_awars(aw_root);
206
207    RootCallback matrix_needs_recalc_callback = makeRootCallback(matrix_needs_recalc_cb);
208    aw_root->awar_int(AWAR_DIST_MATRIX_DNA_ENABLED, 0)->add_callback(matrix_needs_recalc_callback); // user matrix disabled by default
209    {
210        GBDATA *gbd = GB_search(AW_ROOT_DEFAULT, AWAR_DIST_MATRIX_DNA_BASE, GB_FIND);
211        GB_add_callback(gbd, GB_CB_CHANGED, makeDatabaseCallback(matrix_needs_recalc_cb));
212    }
213
214    aw_root->awar_string(AWAR_DIST_WHICH_SPECIES, "marked", def)->add_callback(matrix_needs_recalc_callback);
215    {
216        char *default_ali = GBT_get_default_alignment(db);
217        if (!default_ali) return GB_await_error();
218
219        aw_root->awar_string(AWAR_DIST_ALIGNMENT, "", def)->add_callback(matrix_needs_recalc_callback)->write_string(default_ali);
220        free(default_ali);
221    }
222    aw_root->awar_string(AWAR_DIST_FILTER_ALIGNMENT, "none", def);
223    aw_root->awar_string(AWAR_DIST_FILTER_NAME,      "none", def);
224    aw_root->awar_string(AWAR_DIST_FILTER_FILTER,    "",     def)->add_callback(matrix_needs_recalc_callback);
225    aw_root->awar_int   (AWAR_DIST_FILTER_SIMPLIFY,  0,      def)->add_callback(matrix_needs_recalc_callback);
226
227    aw_root->awar_string(AWAR_DIST_CANCEL_CHARS, ".", def)->add_callback(matrix_needs_recalc_callback);
228    aw_root->awar_int(AWAR_DIST_CORR_TRANS, (int)DI_TRANSFORMATION_SIMILARITY, def)->add_callback(matrix_needs_recalc_callback)->set_minmax(0, DI_TRANSFORMATION_COUNT-1);
229
230    aw_root->awar(AWAR_DIST_FILTER_ALIGNMENT)->map(AWAR_DIST_ALIGNMENT);
231
232    AW_create_fileselection_awars(aw_root, AWAR_DIST_SAVE_MATRIX_BASE, ".", "", "infile");
233    aw_root->awar_int(AWAR_DIST_SAVE_MATRIX_TYPE, 0, def);
234
235    enum treetype { CURR, SORT, COMPRESS, TREEAWARCOUNT };
236    AW_awar *tree_awar[TREEAWARCOUNT] = { NULp, NULp, NULp };
237
238    aw_root->awar_string(AWAR_DIST_TREE_STD_NAME,  "tree_nj", def);
239    {
240        char *currentTree = aw_root->awar_string(AWAR_TREE, "", db)->read_string();
241        tree_awar[CURR]   = aw_root->awar_string(AWAR_DIST_TREE_CURR_NAME, currentTree, def);
242        tree_awar[SORT]   = aw_root->awar_string(AWAR_DIST_TREE_SORT_NAME, currentTree, def);
243        free(currentTree);
244    }
245    tree_awar[COMPRESS] = aw_root->awar_string(AWAR_DIST_TREE_COMP_NAME, NO_TREE_SELECTED, def);
246
247    aw_root->awar_int(AWAR_DIST_BOOTSTRAP_COUNT, 1000, def);
248    aw_root->awar_int(AWAR_DIST_MATRIX_AUTO_RECALC, 0, def)->add_callback(auto_calc_changed_cb);
249    aw_root->awar_int(AWAR_DIST_MATRIX_AUTO_CALC_TREE, 0, def)->add_callback(auto_calc_changed_cb);
250
251    aw_root->awar_float(AWAR_DIST_MIN_DIST, 0.0)->set_minmax(0.0, 4.0);
252    aw_root->awar_float(AWAR_DIST_MAX_DIST, 0.0)->set_minmax(0.0, 4.0);
253
254    aw_root->awar_string(AWAR_SPECIES_NAME, "", db);
255
256    DI_create_cluster_awars(aw_root, def, db);
257
258#if defined(DEBUG)
259    AWT_create_db_browser_awars(aw_root, def);
260#endif // DEBUG
261
262    {
263        GB_push_transaction(db);
264
265        GBDATA *gb_species_data = GBT_get_species_data(db);
266        GB_add_callback(gb_species_data, GB_CB_CHANGED, makeDatabaseCallback(matrix_needs_recalc_cb));
267
268        GB_pop_transaction(db);
269    }
270
271    AWT_registerTreeAwarCallback(tree_awar[CURR],     makeTreeAwarCallback(selected_tree_changed_cb),          true);
272    AWT_registerTreeAwarCallback(tree_awar[SORT],     makeTreeAwarCallback(matrix_needs_recalc_cb),            true);
273    AWT_registerTreeAwarCallback(tree_awar[COMPRESS], makeTreeAwarCallback(compressed_matrix_needs_recalc_cb), true);
274
275    auto_calc_changed_cb(aw_root);
276
277    return NULp;
278}
279
280DI_ENTRY::DI_ENTRY(GBDATA *gbd, DI_MATRIX *phmatrix_)
281    : phmatrix(phmatrix_),
282      full_name(NULp),
283      sequence(NULp),
284      name(NULp),
285      group_nr(0)
286{
287    GBDATA *gb_ali = GB_entry(gbd, phmatrix->get_aliname());
288    if (gb_ali) {
289        GBDATA *gb_data = GB_entry(gb_ali, "data");
290        if (gb_data) {
291            if (phmatrix->is_AA) {
292                sequence = new AP_sequence_simple_protein(phmatrix->get_aliview());
293            }
294            else {
295                sequence = new AP_sequence_parsimony(phmatrix->get_aliview());
296            }
297            sequence->bind_to_species(gbd);
298            sequence->lazy_load_sequence(); // load sequence
299
300            name      = GBT_read_string(gbd, "name");
301            full_name = GBT_read_string(gbd, "full_name");
302        }
303    }
304}
305
306DI_ENTRY::DI_ENTRY(const char *name_, DI_MATRIX *phmatrix_)
307    : phmatrix(phmatrix_),
308      full_name(NULp),
309      sequence(NULp),
310      name(ARB_strdup(name_)),
311      group_nr(0)
312{}
313
314DI_ENTRY::~DI_ENTRY() {
315    delete sequence;
316    free(name);
317    free(full_name);
318}
319
320DI_MATRIX::DI_MATRIX(const AliView& aliview_)
321    : gb_species_data(NULp),
322      seq_len(0),
323      allocated_entries(0),
324      aliview(new AliView(aliview_)),
325      is_AA(false),
326      entries(NULp),
327      nentries(0),
328      matrix(NULp),
329      matrix_type(DI_MATRIX_FULL)
330{
331    memset(cancel_columns, 0, sizeof(cancel_columns));
332}
333
334char *DI_MATRIX::unload() {
335    for (size_t i=0; i<nentries; i++) {
336        delete entries[i];
337    }
338    freenull(entries);
339    nentries = 0;
340    return NULp;
341}
342
343DI_MATRIX::~DI_MATRIX() {
344    unload();
345    delete matrix;
346    delete aliview;
347}
348
349struct TreeOrderedSpecies {
350    GBDATA  *gbd;
351    int      order_index;
352
353    TreeOrderedSpecies(const MatrixOrder& order, GBDATA *gb_spec)
354        : gbd(gb_spec),
355          order_index(order.get_index(GBT_get_name_or_description(gbd)))
356    {}
357};
358
359MatrixOrder::MatrixOrder(GBDATA *gb_main, GB_CSTR sort_tree_name)
360    : name2pos(NULp),
361      leafs(0)
362{
363    if (sort_tree_name) {
364        int       size;
365        TreeNode *sort_tree = GBT_read_tree_and_size(gb_main, sort_tree_name, new SimpleRoot, &size);
366
367        if (sort_tree) {
368            leafs    = size+1;
369            name2pos = GBS_create_hash(leafs, GB_IGNORE_CASE);
370
371            IF_ASSERTION_USED(int leafsLoaded = leafs);
372            leafs = 0;
373            insert_in_hash(sort_tree);
374
375            arb_assert(leafsLoaded == leafs);
376            destroy(sort_tree);
377        }
378        else {
379            GB_clear_error();
380        }
381    }
382}
383static int TreeOrderedSpecies_cmp(const void *p1, const void *p2, void *) {
384    TreeOrderedSpecies *s1 = (TreeOrderedSpecies*)p1;
385    TreeOrderedSpecies *s2 = (TreeOrderedSpecies*)p2;
386
387    return s2->order_index - s1->order_index;
388}
389
390void MatrixOrder::applyTo(TreeOrderedSpecies **species_array, size_t array_size) const {
391    GB_sort((void**)species_array, 0, array_size, TreeOrderedSpecies_cmp, NULp);
392}
393
394static size_t get_load_count(LoadWhat what, GBDATA *gb_main, GBDATA **species_list) {
395    size_t no_of_species = -1U;
396
397    switch (what) {
398        case DI_LOAD_ALL:
399            no_of_species = GBT_get_species_count(gb_main);
400            break;
401        case DI_LOAD_MARKED:
402            no_of_species = GBT_count_marked_species(gb_main);
403            break;
404        case DI_LOAD_LIST:
405            di_assert(species_list);
406            for (no_of_species = 0; species_list[no_of_species]; ++no_of_species) ;
407            break;
408    }
409
410    di_assert(no_of_species != -1U);
411
412    return no_of_species;
413}
414
415GB_ERROR DI_MATRIX::load(LoadWhat what, const MatrixOrder& order, bool show_warnings, GBDATA **species_list) {
416    GBDATA     *gb_main = get_gb_main();
417    const char *use     = get_aliname();
418
419    GB_transaction ta(gb_main);
420
421    seq_len          = GBT_get_alignment_len(gb_main, use); arb_assert(seq_len>0);
422    is_AA            = GBT_is_alignment_protein(gb_main, use);
423    gb_species_data  = GBT_get_species_data(gb_main);
424
425    allocated_entries = 1000;
426    ARB_calloc(entries, allocated_entries);
427
428    nentries = 0;
429
430    size_t no_of_species = get_load_count(what, gb_main, species_list);
431    if (no_of_species<2) {
432        return GBS_global_string("Not enough input species (%zu)", no_of_species);
433    }
434
435    TreeOrderedSpecies *species_to_load[no_of_species];
436
437    {
438        size_t i = 0;
439        switch (what) {
440            case DI_LOAD_ALL: {
441                for (GBDATA *gb_species = GBT_first_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_species(gb_species), ++i) {
442                    species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
443                }
444                break;
445            }
446            case DI_LOAD_MARKED: {
447                for (GBDATA *gb_species = GBT_first_marked_species_rel_species_data(gb_species_data); gb_species; gb_species = GBT_next_marked_species(gb_species), ++i) {
448                    species_to_load[i] = new TreeOrderedSpecies(order, gb_species);
449                }
450                break;
451            }
452            case DI_LOAD_LIST: {
453                for (i = 0; species_list[i]; ++i) {
454                    species_to_load[i] = new TreeOrderedSpecies(order, species_list[i]);
455                }
456                break;
457            }
458        }
459        arb_assert(i == no_of_species);
460    }
461
462    if (order.defined()) {
463        order.applyTo(species_to_load, no_of_species);
464        if (show_warnings) {
465            int species_not_in_sort_tree = 0;
466            for (size_t i = 0; i<no_of_species; ++i) {
467                if (!species_to_load[i]->order_index) {
468                    species_not_in_sort_tree++;
469                }
470            }
471            if (species_not_in_sort_tree) {
472                GBT_message(gb_main, GBS_global_string("Warning: %i of the affected species are not in sort-tree", species_not_in_sort_tree));
473            }
474        }
475    }
476    else {
477        if (show_warnings) {
478            static bool shown = false;
479            if (!shown) { // showing once is enough
480                GBT_message(gb_main, "Warning: No valid tree given to sort matrix (using default database order)");
481                shown = true;
482            }
483        }
484    }
485
486    if (no_of_species>allocated_entries) {
487        allocated_entries = no_of_species;
488        ARB_realloc(entries, allocated_entries);
489    }
490
491    GB_ERROR     error = NULp;
492    arb_progress progress("Preparing sequence data", no_of_species);
493    for (size_t i = 0; i<no_of_species && !error; ++i) {
494        DI_ENTRY *phentry = new DI_ENTRY(species_to_load[i]->gbd, this);
495        if (phentry->sequence) {    // a species found
496            arb_assert(nentries<allocated_entries);
497            entries[nentries++] = phentry;
498        }
499        else {
500            delete phentry;
501        }
502        delete species_to_load[i];
503        species_to_load[i] = NULp;
504
505        progress.inc_and_check_user_abort(error);
506    }
507
508    return error;
509}
510
511char *DI_MATRIX::calculate_overall_freqs(double rel_frequencies[AP_MAX], char *cancel) {
512    di_assert(is_AA == false);
513
514    long hits2[AP_MAX];
515    long sum   = 0;
516    int  i;
517    int  pos;
518    int  b;
519    long s_len = aliview->get_length();
520
521    memset((char *) &hits2[0], 0, sizeof(hits2));
522    for (size_t row = 0; row < nentries; row++) {
523        const char *seq1 = entries[row]->get_nucl_seq()->get_sequence();
524        for (pos = 0; pos < s_len; pos++) {
525            b = *(seq1++);
526            if (cancel[b]) continue;
527            hits2[b]++;
528        }
529    }
530    for (i = 0; i < AP_MAX; i++) { // LOOP_VECTORIZED
531        sum += hits2[i];
532    }
533    for (i = 0; i < AP_MAX; i++) {
534        rel_frequencies[i] = hits2[i] / (double) sum;
535    }
536    return NULp;
537}
538
539double DI_MATRIX::corr(double dist, double b, double & sigma) {
540    const double eps = 0.01;
541    double ar = 1.0 - dist/b;
542    sigma = 1000.0;
543    if (ar< eps) return 3.0;
544    sigma = b/ar;
545    return - b * log(1-dist/b);
546}
547
548GB_ERROR DI_MATRIX::calculate(const char *cancel, DI_TRANSFORMATION transformation, bool *aborted_flag, AP_matrix *userdef_matrix) {
549    di_assert(is_AA == false);
550
551    if (userdef_matrix) {
552        switch (transformation) {
553            case DI_TRANSFORMATION_NONE:
554            case DI_TRANSFORMATION_SIMILARITY:
555            case DI_TRANSFORMATION_JUKES_CANTOR:
556                break;
557            default:
558                aw_message("Sorry: this kind of distance correction does not support a user defined matrix - it will be ignored");
559                userdef_matrix = NULp;
560                break;
561        }
562    }
563
564    matrix = new AP_smatrix(nentries);
565
566    long   s_len = aliview->get_length();
567    long   hits[AP_MAX][AP_MAX];
568    size_t i;
569
570    if (nentries<=1) {
571        return "Not enough species selected to calculate matrix";
572    }
573    memset(&cancel_columns[0], 0, 256);
574
575    for (i=0; cancel[i]; i++) {
576        cancel_columns[safeCharIndex(AP_sequence_parsimony::table[safeCharIndex(cancel[i])])] = 1;
577        UNCOVERED(); // @@@ cover
578    }
579
580    long   columns;
581    double b;
582    long   frequencies[AP_MAX];
583    double rel_frequencies[AP_MAX];
584    double S_square = 0;
585
586    switch (transformation) {
587        case DI_TRANSFORMATION_FELSENSTEIN:
588            this->calculate_overall_freqs(rel_frequencies, cancel_columns);
589            S_square = 0.0;
590            for (i=0; i<AP_MAX; i++) { // LOOP_VECTORIZED[!<8.1]
591                S_square += rel_frequencies[i]*rel_frequencies[i];
592            }
593            break;
594        default:    break;
595    };
596
597    arb_progress progress("Calculating distance matrix", matrix_halfsize(nentries, true));
598    GB_ERROR     error = NULp;
599    for (size_t row = 0; row<nentries && !error; row++) {
600        for (size_t col=0; col<=row && !error; col++) {
601            columns = 0;
602
603            const unsigned char *seq1 = entries[row]->get_nucl_seq()->get_usequence();
604            const unsigned char *seq2 = entries[col]->get_nucl_seq()->get_usequence();
605
606            b = 0.0;
607            switch (transformation) {
608                case  DI_TRANSFORMATION_FROM_TREE:
609                    di_assert(0);
610                    break;
611                case  DI_TRANSFORMATION_JUKES_CANTOR:
612                    b = 0.75;
613                    // fall-through
614                case  DI_TRANSFORMATION_NONE:
615                case  DI_TRANSFORMATION_SIMILARITY: {
616                    double  dist = 0.0;
617                    if (userdef_matrix) {
618                        memset((char *)hits, 0, sizeof(long) * AP_MAX * AP_MAX);
619                        int pos;
620                        for (pos = s_len; pos >= 0; pos--) {
621                            hits[*(seq1++)][*(seq2++)]++;
622                        }
623                        int x, y;
624                        double diffsum = 0.0;
625                        double all_sum = 0.001;
626                        for (x = AP_A; x < AP_MAX; x*=2) {
627                            for (y = AP_A; y < AP_MAX; y*=2) {
628                                if (x==y) {
629                                    all_sum += hits[x][y];
630                                }
631                                else {
632                                    UNCOVERED(); // @@@ cover
633                                    diffsum += hits[x][y] * userdef_matrix->get(x, y);
634                                    all_sum += hits[x][y] * userdef_matrix->get(x, y);
635                                }
636                            }
637                        }
638                        dist = diffsum / all_sum;
639                    }
640                    else {
641                        for (int pos = s_len; pos >= 0; pos--) {
642                            int b1 = *(seq1++);
643                            int b2 = *(seq2++);
644                            if (cancel_columns[b1]) continue;
645                            if (cancel_columns[b2]) continue;
646                            columns++;
647                            if (b1&b2) continue;
648                            dist+=1.0;
649                        }
650                        if (columns == 0) columns = 1;
651                        dist /= columns;
652                    }
653                    if (transformation==DI_TRANSFORMATION_SIMILARITY) {
654                        dist =  (1.0-dist);
655                    }
656                    else if (b) {
657                        double sigma;
658                        dist = this->corr(dist, b, sigma);
659                    }
660                    matrix->set(row, col, dist);
661                    break;
662                }
663                case DI_TRANSFORMATION_KIMURA:
664                case DI_TRANSFORMATION_OLSEN:
665                case DI_TRANSFORMATION_FELSENSTEIN: {
666                    int    pos;
667                    double dist = 0.0;
668                    long   N, P, Q, M;
669                    double p, q;
670
671                    memset((char *)hits, 0, sizeof(long) * AP_MAX * AP_MAX);
672                    for (pos = s_len; pos >= 0; pos--) {
673                        hits[*(seq1++)][*(seq2++)]++;
674                    }
675                    switch (transformation) {
676                        case DI_TRANSFORMATION_KIMURA:
677                            P = hits[AP_A][AP_G] +
678                                hits[AP_G][AP_A] +
679                                hits[AP_C][AP_T] +
680                                hits[AP_T][AP_C];
681                            Q = hits[AP_A][AP_C] +
682                                hits[AP_A][AP_T] +
683                                hits[AP_C][AP_A] +
684                                hits[AP_T][AP_A] +
685                                hits[AP_G][AP_C] +
686                                hits[AP_G][AP_T] +
687                                hits[AP_C][AP_G] +
688                                hits[AP_T][AP_G];
689                            M = hits[AP_A][AP_A] +
690                                hits[AP_C][AP_C] +
691                                hits[AP_G][AP_G] +
692                                hits[AP_T][AP_T];
693                            N = P+Q+M;
694                            if (N==0) N=1;
695                            p = (double)P/(double)N;
696                            q = (double)Q/(double)N;
697                            dist = - .5 * log(
698                                              (1.0-2.0*p-q)*sqrt(1.0-2.0*q)
699                                              );
700                            break;
701
702                        case DI_TRANSFORMATION_OLSEN:
703                        case DI_TRANSFORMATION_FELSENSTEIN:
704
705                            memset((char *)frequencies, 0,
706                                   sizeof(long) * AP_MAX);
707
708                            N = 0;
709                            M = 0;
710
711                            for (i=0; i<AP_MAX; i++) {
712                                if (cancel_columns[i]) continue;
713                                unsigned int j;
714                                for (j=0; j<i; j++) {
715                                    if (cancel_columns[j]) continue;
716                                    frequencies[i] +=
717                                        hits[i][j]+
718                                        hits[j][i];
719                                }
720                                frequencies[i] += hits[i][i];
721                                N += frequencies[i];
722                                M += hits[i][i];
723                            }
724                            if (N==0) N=1;
725                            if (transformation == DI_TRANSFORMATION_OLSEN) { // Calc sum square freq individually for each line
726                                S_square = 0.0;
727                                for (i=0; i<AP_MAX; i++) S_square += frequencies[i]*frequencies[i];
728                                b = 1.0 - S_square/((double)N*(double)N);
729                            }
730                            else {
731                                b = 1.0 - S_square;
732                            }
733
734                            dist = ((double)(N-M)) / (double) N;
735                            double sigma;
736                            dist = this->corr(dist, b, sigma);
737                            break;
738
739                        default: return "Sorry: Transformation not implemented";
740                    }
741                    matrix->set(row, col, dist);
742                    break;
743                }
744                default:;
745            }   // switch
746            progress.inc_and_check_user_abort(error);
747        }
748    }
749    if (aborted_flag && progress.aborted()) *aborted_flag = true;
750    return error;
751}
752
753GB_ERROR DI_MATRIX::calculate_pro(DI_TRANSFORMATION transformation, bool *aborted_flag) {
754    di_assert(is_AA == true);
755
756    di_cattype catType;
757    switch (transformation) {
758        case DI_TRANSFORMATION_NONE:                catType = NONE;       break;
759        case DI_TRANSFORMATION_SIMILARITY:          catType = SIMILARITY; break;
760        case DI_TRANSFORMATION_KIMURA:              catType = KIMURA;     break;
761        case DI_TRANSFORMATION_PAM:                 catType = PAM;        break;
762        case DI_TRANSFORMATION_CATEGORIES_HALL:     catType = HALL;       break;
763        case DI_TRANSFORMATION_CATEGORIES_BARKER:   catType = GEORGE;     break;
764        case DI_TRANSFORMATION_CATEGORIES_CHEMICAL: catType = CHEMICAL;   break;
765        default:
766            return "This correction is not available for protein data";
767    }
768    matrix = new AP_smatrix(nentries);
769
770    di_protdist prodist(UNIVERSAL, catType, nentries, entries, aliview->get_length(), matrix);
771    return prodist.makedists(aborted_flag);
772}
773
774typedef std::map<const char*, TreeNode*, charpLess> NamedNodes;
775
776GB_ERROR link_to_tree(NamedNodes& named, TreeNode *node) {
777    GB_ERROR error = NULp;
778    if (node->is_leaf()) {
779        NamedNodes::iterator found = named.find(node->name);
780        if (found != named.end()) {
781            if (found->second) {
782                error = GBS_global_string("Invalid tree (two nodes named '%s')", node->name);
783            }
784            else {
785                found->second = node;
786            }
787        }
788        // otherwise, we do not care about the node (e.g. because it is not marked)
789    }
790    else {
791        error             = link_to_tree(named, node->get_leftson());
792        if (!error) error = link_to_tree(named, node->get_rightson());
793    }
794    return error;
795}
796
797#if defined(ASSERTION_USED)
798static TreeNode *findNode(TreeNode *node, const char *name) {
799    if (node->is_leaf()) {
800        return strcmp(node->name, name) == 0 ? node : NULp;
801    }
802
803    TreeNode *found   = findNode(node->get_leftson(), name);
804    if (!found) found = findNode(node->get_rightson(), name);
805    return found;
806}
807#endif
808
809static GB_ERROR init(NamedNodes& node, TreeNode *tree, const DI_ENTRY*const*const entries, size_t nentries) {
810    GB_ERROR error = NULp;
811    for (size_t n = 0; n<nentries; ++n) {
812        node[entries[n]->name] = NULp;
813    }
814    error = link_to_tree(node, tree);
815    if (!error) { // check for missing species (needed but not in tree)
816        size_t      missing     = 0;
817        const char *exampleName = NULp;
818
819        for (size_t n = 0; n<nentries; ++n) {
820            NamedNodes::iterator found = node.find(entries[n]->name);
821            if (found == node.end()) {
822                ++missing;
823                exampleName = entries[n]->name;
824            }
825            else {
826                di_assert(node[entries[n]->name] == findNode(tree, entries[n]->name));
827                if (!node[entries[n]->name]) {
828                    ++missing;
829                    exampleName = entries[n]->name;
830                }
831            }
832        }
833
834        if (missing) {
835            error = GBS_global_string("Tree is missing %zu required species (e.g. '%s')", missing, exampleName);
836        }
837    }
838    return error;
839}
840
841GB_ERROR DI_MATRIX::extract_from_tree(const char *treename, bool *aborted_flag) {
842    GB_ERROR error         = NULp;
843    if (nentries<=1) error = "Not enough species selected to calculate matrix";
844    else {
845        TreeNode *tree;
846        {
847            GB_transaction ta(get_gb_main());
848            tree = GBT_read_tree(get_gb_main(), treename, new SimpleRoot);
849        }
850        if (!tree) error = GB_await_error();
851        else {
852            arb_progress progress("Extracting distances from tree", matrix_halfsize(nentries, true));
853            NamedNodes   node;
854
855            error  = init(node, tree, entries, nentries);
856            matrix = new AP_smatrix(nentries);
857
858            for (size_t row = 0; row<nentries && !error; row++) {
859                TreeNode *rnode = node[entries[row]->name];
860                for (size_t col=0; col<=row && !error; col++) {
861                    double dist;
862                    if (col != row) {
863                        TreeNode *cnode = node[entries[col]->name];
864                        dist  = rnode->intree_distance_to(cnode);
865                    }
866                    else {
867                        dist = 0.0;
868                    }
869                    matrix->set(row, col, dist);
870                    progress.inc_and_check_user_abort(error);
871                }
872            }
873            UNCOVERED();
874            destroy(tree);
875            if (aborted_flag && progress.aborted()) *aborted_flag = true;
876            if (error) progress.done();
877        }
878    }
879    return error;
880}
881
882inline LoadWhat whatToLoad() {
883    const char *load_what = AW_root::SINGLETON->awar(AWAR_DIST_WHICH_SPECIES)->read_char_pntr();
884    return (strcmp(load_what, "all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
885}
886
887__ATTR__USERESULT static GB_ERROR di_calculate_matrix(AW_root *aw_root, const WeightedFilter *weighted_filter, bool bootstrap_flag, bool show_warnings, bool *aborted_flag) {
888    // sets 'aborted_flag' to true, if it is non-NULp and the calculation has been aborted
889    GB_ERROR  error   = NULp;
890
891    if (GLOBAL_MATRIX.exists()) {
892        di_assert(!need_recalc.matrix);
893    }
894    else {
895        GBDATA         *gb_main = weighted_filter->get_gb_main();
896        GB_transaction  ta(gb_main);
897
898        char *use     = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
899        long  ali_len = GBT_get_alignment_len(gb_main, use);
900
901        if (ali_len<=0) {
902            error = "Please select a valid alignment";
903            GB_clear_error();
904        }
905        else {
906            const LoadWhat what          = whatToLoad();
907            size_t         species_count = get_load_count(what, gb_main, NULp);
908            double         phase1_fraction;
909            {
910                // keep sync with .@MatrixSteps
911                double O_loadData   = 10*species_count;
912                double O_calcMatrix = triangular_number(species_count);
913                phase1_fraction     = O_loadData / (O_loadData+O_calcMatrix);
914            }
915
916            arb_progress progress(WEIGHTED, "Calculating matrix", phase1_fraction); // two steps: 1x load data + 1x calc matrix
917
918            char    *cancel  = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
919            AliView *aliview = weighted_filter->create_aliview(use, error);
920
921            if (!error) {
922                if (bootstrap_flag) aliview->get_filter()->enable_bootstrap();
923
924                const char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_char_pntr();
925                {
926                    DI_MATRIX *phm   = new DI_MATRIX(*aliview);
927                    phm->matrix_type = DI_MATRIX_FULL;
928
929                    static SmartCharPtr          last_sort_tree_name;
930                    static SmartPtr<MatrixOrder> last_order;
931
932                    if (last_sort_tree_name.isNull() || !sort_tree_name || strcmp(&*last_sort_tree_name, sort_tree_name) != 0) {
933                        last_sort_tree_name = nulldup(sort_tree_name);
934                        last_order = new MatrixOrder(gb_main, sort_tree_name);
935                    }
936                    di_assert(last_order.isSet());
937                    error = phm->load(what, *last_order, show_warnings, NULp); // step1 of progress
938                    progress.inc_and_check_user_abort(error);
939                    bool aborted = error ? progress.aborted() : false;
940                    error = ta.close(error);
941
942                    if (!error) {
943                        DI_TRANSFORMATION trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
944
945                        if (trans == DI_TRANSFORMATION_FROM_TREE) {
946                            const char *treename = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_char_pntr();
947                            error                = phm->extract_from_tree(treename, &aborted);
948                        }
949                        else {
950                            if (phm->is_AA) error = phm->calculate_pro(trans, &aborted);
951                            else            error = phm->calculate(cancel, trans, &aborted, get_user_matrix());
952                        }
953                    }
954
955                    if (aborted) {
956                        di_assert(error);
957                        if (aborted_flag) *aborted_flag = true;
958                    }
959                    if (error) {
960                        delete phm;
961                        GLOBAL_MATRIX.forget();
962                    }
963                    else {
964                        GLOBAL_MATRIX.replaceBy(phm);
965                        tree_needs_recalc_cb();
966                        need_recalc.matrix = false;
967                    }
968                }
969            }
970
971            free(cancel);
972            delete aliview;
973        }
974        free(use);
975
976        di_assert(contradicted(error, GLOBAL_MATRIX.exists()));
977    }
978    return error;
979}
980
981static void di_mark_by_distance(AW_window *aww, WeightedFilter *weighted_filter) {
982    AW_root *aw_root    = aww->get_root();
983    float    lowerBound = aw_root->awar(AWAR_DIST_MIN_DIST)->read_float();
984    float    upperBound = aw_root->awar(AWAR_DIST_MAX_DIST)->read_float();
985
986    GB_ERROR error = NULp;
987    if (lowerBound >= upperBound) {
988        error = GBS_global_string("Lower bound (%f) has to be smaller than upper bound (%f)", lowerBound, upperBound);
989    }
990    else if (lowerBound<0.0 || lowerBound > 1.0) {
991        error = GBS_global_string("Lower bound (%f) is not in allowed range [0.0 .. 1.0]", lowerBound);
992    }
993    else if (upperBound<0.0 || upperBound > 1.0) {
994        error = GBS_global_string("Upper bound (%f) is not in allowed range [0.0 .. 1.0]", upperBound);
995    }
996    else {
997        GBDATA         *gb_main = weighted_filter->get_gb_main();
998        GB_transaction  ta(gb_main);
999
1000        char *selected = aw_root->awar(AWAR_SPECIES_NAME)->read_string();
1001        if (!selected[0]) {
1002            error = "Please select a species";
1003        }
1004        else {
1005            GBDATA *gb_selected = GBT_find_species(gb_main, selected);
1006            if (!gb_selected) {
1007                error = GBS_global_string("Couldn't find species '%s'", selected);
1008            }
1009            else {
1010                char              *use     = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1011                char              *cancel  = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1012                DI_TRANSFORMATION  trans   = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1013                AliView           *aliview = weighted_filter->create_aliview(use, error);
1014
1015                if (!error) {
1016                    DI_MATRIX *prev_global = GLOBAL_MATRIX.swap(NULp);
1017
1018                    size_t speciesCount   = GBT_get_species_count(gb_main);
1019                    bool   markedSelected = false;
1020
1021                    arb_progress progress("Mark species by distance", speciesCount);
1022                    MatrixOrder  order(gb_main, NULp);
1023
1024                    for (GBDATA *gb_species = GBT_first_species(gb_main);
1025                         gb_species && !error;
1026                         gb_species = GBT_next_species(gb_species))
1027                    {
1028                        DI_MATRIX *phm         = new DI_MATRIX(*aliview);
1029                        phm->matrix_type       = DI_MATRIX_FULL;
1030                        GBDATA *species_pair[] = { gb_selected, gb_species, NULp };
1031
1032                        error = phm->load(DI_LOAD_LIST, order, false, species_pair);
1033
1034                        if (phm->nentries == 2) { // if species has no alignment -> nentries<2
1035                            if (!error) {
1036                                if (phm->is_AA) error = phm->calculate_pro(trans, NULp);
1037                                else            error = phm->calculate(cancel, trans, NULp, get_user_matrix());
1038                            }
1039
1040                            if (!error) {
1041                                double dist_value = phm->matrix->get(0, 1);                         // distance or conformance
1042                                bool   mark       = (lowerBound <= dist_value && dist_value <= upperBound);
1043                                GB_write_flag(gb_species, mark);
1044
1045                                if (!markedSelected) {
1046                                    dist_value = phm->matrix->get(0, 0);                                     // distance or conformance to self
1047                                    mark       = (lowerBound <= dist_value && dist_value <= upperBound);
1048                                    GB_write_flag(gb_selected, mark);
1049
1050                                    markedSelected = true;
1051                                }
1052                            }
1053                        }
1054
1055                        delete phm;
1056                        if (!error) progress.inc_and_check_user_abort(error);
1057                    }
1058
1059                    di_assert(!GLOBAL_MATRIX.exists());
1060                    ASSERT_RESULT(DI_MATRIX*, NULp, GLOBAL_MATRIX.swap(prev_global));
1061
1062                    if (error) progress.done();
1063                }
1064
1065                delete aliview;
1066                free(cancel);
1067                free(use);
1068            }
1069        }
1070
1071        free(selected);
1072        error = ta.close(error);
1073    }
1074
1075    if (error) {
1076        aw_message(error);
1077    }
1078}
1079
1080static GB_ERROR di_recalc_matrix() {
1081    // recalculate matrix
1082    last_matrix_calculation_error = NULp;
1083    if (need_recalc.matrix && GLOBAL_MATRIX.exists()) {
1084        GLOBAL_MATRIX.forget();
1085    }
1086    di_assert(recalculate_matrix_cb.isSet());
1087    (*recalculate_matrix_cb)();
1088    return last_matrix_calculation_error;
1089}
1090
1091static void di_view_matrix_cb(AW_window *aww, save_matrix_params *sparam) {
1092    GB_ERROR error = di_recalc_matrix();
1093    if (!error) {
1094        if (!matrixDisplay) matrixDisplay = new MatrixDisplay(sparam->weighted_filter->get_gb_main());
1095
1096        static AW_window *viewer = NULp;
1097        if (!viewer) viewer = DI_create_view_matrix_window(aww->get_root(), matrixDisplay, sparam);
1098
1099        matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
1100        matrixDisplay->update_display();
1101
1102        GLOBAL_MATRIX.set_changed_cb(matrix_changed_cb);
1103
1104        viewer->activate();
1105    }
1106}
1107
1108static void di_save_matrix_cb(AW_window *aww) {
1109    // save the matrix
1110    GB_ERROR error = di_recalc_matrix();
1111    if (!error) {
1112        char              *filename = aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_FILENAME)->read_string();
1113        enum DI_SAVE_TYPE  type     = (enum DI_SAVE_TYPE)aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_TYPE)->read_int();
1114
1115        GLOBAL_MATRIX.get()->save(filename, type);
1116        free(filename);
1117    }
1118    AW_refresh_fileselection(aww->get_root(), AWAR_DIST_SAVE_MATRIX_BASE);
1119    aww->hide_or_notify(error);
1120}
1121
1122AW_window *DI_create_save_matrix_window(AW_root *aw_root, save_matrix_params *save_params) {
1123    static AW_window_simple *aws = NULp;
1124    if (!aws) {
1125        aws = new AW_window_simple;
1126        aws->init(aw_root, "SAVE_MATRIX", "Save Matrix");
1127        aws->load_xfig("sel_box_user.fig");
1128
1129        aws->at("close");
1130        aws->callback(AW_POPDOWN);
1131        aws->create_button("CLOSE", "CANCEL", "C");
1132
1133
1134        aws->at("help"); aws->callback(makeHelpCallback("save_matrix.hlp"));
1135        aws->create_button("HELP", "HELP", "H");
1136
1137        aws->at("user");
1138        aws->create_option_menu(AWAR_DIST_SAVE_MATRIX_TYPE);
1139        aws->insert_default_option("Phylip Format (Lower Triangular Matrix)", "P", DI_SAVE_PHYLIP_COMP);
1140        aws->insert_option("Readable (using NDS)", "R", DI_SAVE_READABLE);
1141        aws->insert_option("Tabbed (using NDS)", "R", DI_SAVE_TABBED);
1142        aws->update_option_menu();
1143
1144        AW_create_standard_fileselection(aws, save_params->awar_base);
1145
1146        aws->at("save2");
1147        aws->callback(makeWindowCallback(di_save_matrix_cb));
1148        aws->create_button("SAVE", "SAVE", "S");
1149
1150        aws->at("cancel2");
1151        aws->callback(AW_POPDOWN);
1152        aws->create_button("CLOSE", "CANCEL", "C");
1153    }
1154    return aws;
1155}
1156
1157static AW_window *awt_create_select_cancel_window(AW_root *aw_root) {
1158    AW_window_simple *aws = new AW_window_simple;
1159    aws->init(aw_root, "SELECT_CHARS_TO_CANCEL_COLUMN", "CANCEL SELECT");
1160    aws->load_xfig("di_cancel.fig");
1161
1162    aws->at("close");
1163    aws->callback(AW_POPDOWN);
1164    aws->create_button("CLOSE", "CLOSE", "C");
1165
1166    aws->at("cancel");
1167    aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1168
1169    return aws;
1170}
1171
1172static const char *enum_trans_to_string[] = {
1173    "none",
1174    "similarity",
1175    "jukes_cantor",
1176    "felsenstein",
1177
1178    "pam",
1179    "hall",
1180    "barker",
1181    "chemical",
1182
1183    "kimura",
1184    "olsen",
1185    "felsenstein voigt",
1186    "olsen voigt",
1187    "max ml",
1188
1189    NULp, // treedist
1190};
1191
1192STATIC_ASSERT(ARRAY_ELEMS(enum_trans_to_string) == DI_TRANSFORMATION_COUNT);
1193
1194static void di_calculate_tree_cb(AW_window *aww, WeightedFilter *weighted_filter, bool bootstrap_flag) {
1195    recalculate_tree_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_tree_cb, weighted_filter, bootstrap_flag));
1196
1197    AW_root  *aw_root   = aww->get_root();
1198    GB_ERROR  error     = NULp;
1199    StrArray *all_names = NULp;
1200
1201    long loop_count      = 0;
1202    long bootstrap_count = aw_root->awar(AWAR_DIST_BOOTSTRAP_COUNT)->read_int();
1203
1204    {
1205        char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1206        error           = GBT_check_tree_name(tree_name);
1207        free(tree_name);
1208    }
1209
1210    SmartPtr<arb_progress>  progress, tree_progress;
1211    SmartPtr<ConsensusTree> ctree;
1212
1213    GBDATA *gb_main = weighted_filter->get_gb_main();
1214
1215    double phase1_fraction;
1216    {
1217        const size_t species_count = get_load_count(whatToLoad(), gb_main, NULp);
1218
1219        double O_calcMatrix = 60 * (triangular_number(species_count) + 10*species_count); // keep sync with .@MatrixSteps
1220        double O_joinTree   = 1  * sum_of_triangular_numbers(species_count);              // keep sync with ../SL/NEIGHBOURJOIN/NJ.cxx@NJsteps
1221
1222        phase1_fraction = O_calcMatrix / (O_calcMatrix+O_joinTree);
1223    }
1224
1225    if (!error) {
1226        if (bootstrap_flag) {
1227            if (bootstrap_count) {
1228                progress = new arb_progress("Calculating bootstrap trees", bootstrap_count+1); // 1 for each bootstrap-tree + 1 for final consensus tree
1229            }
1230            else {
1231                progress = new arb_progress("Calculating bootstrap trees (KILL to stop)", long(INT_MAX));
1232            }
1233            progress->auto_subtitles("tree");
1234        }
1235        else {
1236            progress = new arb_progress("Calculating tree");
1237        }
1238
1239        // @@@ handle case where matrix is already up-to-date?
1240        // @@@ if need_recalc.matrix -> use simple wrapping progress (using 1 step); do the same below at 2nd construction point
1241        tree_progress = new arb_progress(WEIGHTED, NULp, phase1_fraction); // each tree calculation consists of 2 steps (matrix+neighbourjoining)
1242
1243        if (bootstrap_flag) {
1244            GLOBAL_MATRIX.forget();
1245            GLOBAL_MATRIX.set_changed_cb(NULp); // otherwise matrix window will repeatedly pop up/down
1246
1247            error = di_calculate_matrix(aw_root, weighted_filter, bootstrap_flag, true, NULp);
1248            if (!error) {
1249                DI_MATRIX *matr = GLOBAL_MATRIX.get();
1250                if (!matr) {
1251                    error = "unexpected error in di_calculate_matrix_cb (data missing)";
1252                }
1253                else {
1254                    all_names = new StrArray;
1255                    all_names->reserve(matr->nentries+2);
1256
1257                    for (size_t i=0; i<matr->nentries; i++) {
1258                        all_names->put(ARB_strdup(matr->entries[i]->name));
1259                    }
1260                    ctree = new ConsensusTree(*all_names);
1261                }
1262            }
1263        }
1264    }
1265
1266    TreeNode *tree    = NULp;
1267    bool      aborted = false;
1268
1269    do {
1270        if (error) break;
1271
1272        aborted = false;
1273
1274        if (bootstrap_flag) {
1275            if (loop_count>0) { // in first loop we already have a valid matrix -> no need to recalculate
1276                GLOBAL_MATRIX.forget();
1277            }
1278        }
1279        else if (need_recalc.matrix) {
1280            GLOBAL_MATRIX.forget();
1281        }
1282
1283        if (tree_progress.isNull()) tree_progress = new arb_progress(WEIGHTED, NULp, phase1_fraction); // each tree calculation consists of 2 steps (matrix+neighbourjoining)
1284
1285        error = di_calculate_matrix(aw_root, weighted_filter, bootstrap_flag, !bootstrap_flag, &aborted); // [implicit progress increment only if matrix is recalculated]
1286        if (error && aborted) {
1287            error = NULp;       // clear error (otherwise no tree will be read below)
1288            break;              // end of bootstrap
1289        }
1290        tree_progress->inc();
1291
1292        if (!GLOBAL_MATRIX.exists()) {
1293            error = "unexpected error in di_calculate_matrix_cb (data missing)";
1294            break;
1295        }
1296
1297        DI_MATRIX  *matr  = GLOBAL_MATRIX.get();
1298        char      **names = ARB_calloc<char*>(matr->nentries+2);
1299
1300        for (size_t i=0; i<matr->nentries; i++) {
1301            names[i] = matr->entries[i]->name;
1302        }
1303        di_assert(matr->nentries == matr->matrix->size());
1304        tree = neighbourjoining(names, *matr->matrix, new SimpleRoot, &aborted); // [implicit progress increment]
1305        tree_progress->inc();
1306
1307        if (bootstrap_flag && !aborted) {
1308            error               = ctree->insert_tree_weighted(tree, matr->nentries, 1, false);
1309            UNCOVERED();
1310            destroy(tree); tree = NULp;
1311            // Warning: nodes stored as origin in deconstructed PARTs now point to deleted nodes.
1312
1313            if (!error) {
1314                ++loop_count;
1315                progress->inc();
1316            }
1317            if (!bootstrap_count) { // when waiting for kill
1318                time_t        t     = time(NULp);
1319                static time_t tlast = 0;
1320
1321                double diff_seconds = difftime(t, tlast);
1322                if (diff_seconds>3) {
1323                    progress->force_update();
1324                    tlast = t;
1325                }
1326            }
1327        }
1328        free(names);
1329
1330        if (aborted || error) {
1331            break;
1332        }
1333
1334        tree_progress.setNull();
1335
1336    } while (bootstrap_flag && loop_count != bootstrap_count);
1337
1338    if (aborted || error) tree_progress->done();
1339    tree_progress.setNull();
1340
1341    if (!error) {
1342        if (bootstrap_flag) {
1343            if (aborted) {
1344                // continuing to use progress would result in instant abort.
1345                // (abort state is already set in progress-implementation)
1346                // @@@L would need a mechanism to reset that 'abort' flag!
1347
1348                progress->done();
1349                progress.setNull();
1350                progress = new arb_progress("Collect consensus of bootstrapped trees", 1UL);
1351            }
1352
1353            tree = ctree->get_consensus_tree(error);
1354            progress->inc();
1355            if (!error) {
1356                error = GBT_is_invalid(tree);
1357                di_assert(!error);
1358            }
1359        }
1360        else {
1361            error = progress->error_if_aborted();
1362        }
1363
1364        if (!error) {
1365            char   *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1366            GB_begin_transaction(gb_main);
1367            error             = GBT_write_tree(gb_main, tree_name, tree);
1368
1369            if (!error) {
1370                char *filter_name = AWT_get_combined_filter_name(aw_root, "dist");
1371                int   transr      = aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1372
1373                {
1374                    const char *comment;
1375                    if (enum_trans_to_string[transr]) {
1376                        comment = GBS_global_string("PRG=dnadist CORR=%s FILTER=%s PKG=ARB", enum_trans_to_string[transr], filter_name);
1377                    }
1378                    else {
1379                        di_assert(transr == DI_TRANSFORMATION_FROM_TREE);
1380                        const char *treename = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_char_pntr();
1381                        comment = GBS_global_string("PRG=treedist (from '%s') PKG=ARB", treename);
1382                    }
1383
1384                    error = GBT_write_tree_remark(gb_main, tree_name, comment);
1385                }
1386
1387                if (!error) {
1388                    const char *log_message;
1389                    if (bootstrap_flag) {
1390                        log_message = GBS_global_string("generated consensus tree from %li bootstrap trees%s.",
1391                                                        loop_count,
1392                                                        aborted ? " (aborted by user)" : "");
1393                    }
1394                    else {
1395                        log_message = "calculation finished";
1396                    }
1397                    error = GBT_log_to_named_trees_remark(gb_main, tree_name, log_message, true);
1398                }
1399
1400                free(filter_name);
1401            }
1402            error = GB_end_transaction(gb_main, error);
1403            free(tree_name);
1404        }
1405    }
1406
1407    UNCOVERED();
1408    destroy(tree);
1409
1410    // aw_status(); // remove 'abort' flag (@@@ got no equiv for arb_progress yet. really needed?)
1411
1412    if (bootstrap_flag) {
1413        if (all_names) delete all_names;
1414        GLOBAL_MATRIX.forget();
1415    }
1416#if defined(DEBUG)
1417    else {
1418        di_assert(!all_names);
1419    }
1420#endif // DEBUG
1421
1422    if (progress.isSet()) {
1423        progress->done();
1424    }
1425
1426    if (error) {
1427        aw_message(error);
1428    }
1429    else {
1430        need_recalc.tree = false;
1431        aw_root->awar(AWAR_TREE_REFRESH)->touch();
1432    }
1433}
1434
1435
1436static void di_autodetect_callback(AW_window *aww, GBDATA *gb_main) {
1437    GB_push_transaction(gb_main);
1438
1439    GLOBAL_MATRIX.forget();
1440
1441    AW_root  *aw_root = aww->get_root();
1442    char     *use     = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1443    long      ali_len = GBT_get_alignment_len(gb_main, use);
1444    GB_ERROR  error   = NULp;
1445
1446    if (ali_len<=0) {
1447        GB_pop_transaction(gb_main);
1448        error = "Please select a valid alignment";
1449        GB_clear_error();
1450    }
1451    else {
1452        arb_progress progress("Analyzing data");
1453
1454        char *filter_str = aw_root->awar(AWAR_DIST_FILTER_FILTER)->read_string();
1455        // char *cancel     = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1456
1457        AliView *aliview = NULp;
1458        {
1459            AP_filter *ap_filter = NULp;
1460            long       flen  = strlen(filter_str);
1461
1462            if (flen == ali_len) {
1463                ap_filter = new AP_filter(filter_str, "0", ali_len);
1464            }
1465            else {
1466                if (flen) {
1467                    aw_message("Warning: your filter len is not equal to the alignment len\nfilter got truncated with zeros or cutted");
1468                    ap_filter = new AP_filter(filter_str, "0", ali_len);
1469                }
1470                else {
1471                    ap_filter = new AP_filter(ali_len); // unfiltered
1472                }
1473            }
1474
1475            error = ap_filter->is_invalid();
1476            if (!error) {
1477                AP_weights ap_weights(ap_filter);
1478                aliview = new AliView(gb_main, *ap_filter, ap_weights, use);
1479            }
1480            delete ap_filter;
1481        }
1482
1483        if (error) {
1484            GB_pop_transaction(gb_main);
1485        }
1486        else {
1487            DI_MATRIX phm(*aliview);
1488
1489            {
1490                const char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_char_pntr();
1491
1492                GLOBAL_MATRIX.forget();
1493
1494                MatrixOrder order(gb_main, sort_tree_name);
1495                error = phm.load(whatToLoad(), order, true, NULp);
1496            }
1497
1498            GB_pop_transaction(gb_main);
1499
1500            if (!error) {
1501                progress.subtitle("Search Correction");
1502
1503                string msg;
1504                DI_TRANSFORMATION detected = phm.detect_transformation(msg);
1505                aw_root->awar(AWAR_DIST_CORR_TRANS)->write_int(detected);
1506                aw_message(msg.c_str());
1507            }
1508        }
1509
1510        // free(cancel);
1511        delete aliview;
1512
1513        free(filter_str);
1514    }
1515
1516    if (error) aw_message(error);
1517
1518    free(use);
1519}
1520
1521__ATTR__NORETURN static void di_exit(AW_window *aww, GBDATA *gb_main) {
1522    if (gb_main) {
1523        AW_root *aw_root = aww->get_root();
1524        shutdown_macro_recording(aw_root);
1525        aw_root->unlink_awars_from_DB(gb_main);
1526        GB_close(gb_main);
1527    }
1528    GLOBAL_MATRIX.set_changed_cb(NULp);
1529    exit(EXIT_SUCCESS);
1530}
1531
1532static void di_calculate_full_matrix_cb(AW_window *aww, const WeightedFilter *weighted_filter) {
1533    recalculate_matrix_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1534
1535    GLOBAL_MATRIX.forget_if_not_has_type(DI_MATRIX_FULL);
1536    GB_ERROR error = di_calculate_matrix(aww->get_root(), weighted_filter, 0, true, NULp);
1537    aw_message_if(error);
1538    last_matrix_calculation_error = error;
1539    if (!error) tree_needs_recalc_cb();
1540}
1541
1542static void di_calculate_compressed_matrix_cb(AW_window *aww, WeightedFilter *weighted_filter) {
1543    recalculate_matrix_cb = new BoundWindowCallback(aww, makeWindowCallback(di_calculate_compressed_matrix_cb, weighted_filter));
1544
1545    GBDATA         *gb_main = weighted_filter->get_gb_main();
1546    GB_transaction  ta(gb_main);
1547
1548    AW_root  *aw_root  = aww->get_root();
1549    char     *treename = aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->read_string();
1550    GB_ERROR  error    = NULp;
1551    TreeNode  *tree    = GBT_read_tree(gb_main, treename, new SimpleRoot);
1552
1553    if (!tree) {
1554        error = GB_await_error();
1555    }
1556    else {
1557        {
1558            LocallyModify<MatrixDisplay*> skipRefresh(matrixDisplay, NULp); // skip refresh, until matrix has been compressed
1559
1560            GLOBAL_MATRIX.forget(); // always forget (as tree might have changed)
1561            error = di_calculate_matrix(aw_root, weighted_filter, 0, true, NULp);
1562            if (!error && !GLOBAL_MATRIX.exists()) {
1563                error = "Failed to calculate your matrix (bug?)";
1564            }
1565            if (!error) {
1566                error = GLOBAL_MATRIX.get()->compress(tree);
1567            }
1568        }
1569        UNCOVERED();
1570        destroy(tree);
1571
1572        // now force refresh
1573        if (matrixDisplay) {
1574            matrixDisplay->mark(MatrixDisplay::NEED_SETUP);
1575            matrixDisplay->update_display();
1576        }
1577    }
1578    free(treename);
1579    aw_message_if(error);
1580    last_matrix_calculation_error = error;
1581    if (!error) tree_needs_recalc_cb();
1582}
1583
1584static void di_define_sort_tree_name_cb(AW_window *aww) {
1585    AW_root *aw_root = aww->get_root();
1586    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1587    aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->write_string(tree_name);
1588    free(tree_name);
1589}
1590static void di_define_compression_tree_name_cb(AW_window *aww) {
1591    AW_root *aw_root = aww->get_root();
1592    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1593    aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->write_string(tree_name);
1594    free(tree_name);
1595}
1596
1597static void di_define_save_tree_name_cb(AW_window *aww) {
1598    AW_root *aw_root = aww->get_root();
1599    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1600    aw_root->awar(AWAR_DIST_TREE_STD_NAME)->write_string(tree_name);
1601    free(tree_name);
1602}
1603
1604
1605AW_window *DI_create_matrix_window(AW_root *aw_root, GBDATA *gb_main) {
1606    AW_window_simple_menu *aws = new AW_window_simple_menu;
1607    aws->init(aw_root, "NEIGHBOUR JOINING", "NEIGHBOUR JOINING [ARB_DIST]");
1608    aws->load_xfig("di_ge_ma.fig");
1609    aws->button_length(10);
1610
1611    WindowCallback close_cb = makeWindowCallback(di_exit, gb_main);
1612
1613    aws->at("close");
1614    aws->callback(close_cb);
1615    aws->create_button("CLOSE", "CLOSE", "C");
1616
1617    aws->at("help");
1618    aws->callback(makeHelpCallback("dist.hlp"));
1619    aws->create_button("HELP", "HELP", "H");
1620
1621    GB_push_transaction(gb_main);
1622
1623#if defined(DEBUG)
1624    AWT_create_debug_menu(aws);
1625#endif // DEBUG
1626
1627    aws->create_menu("File", "F", AWM_ALL);
1628    insert_macro_menu_entry(aws, false);
1629    aws->insert_menu_topic("quit", "Quit", "Q", "quit.hlp", AWM_ALL, close_cb);
1630
1631    aws->create_menu("Properties", "P", AWM_ALL);
1632    aws->insert_menu_topic("frame_props", "Frame settings ...", "F", "props_frame.hlp", AWM_ALL, AW_preset_window);
1633    aws->sep______________();
1634    AW_insert_common_property_menu_entries(aws);
1635    aws->sep______________();
1636    aws->insert_menu_topic("save_props", "Save Properties (dist.arb)", "S", "savedef.hlp", AWM_ALL, AW_save_properties);
1637
1638    aws->insert_help_topic("ARB_DIST help", "D", "dist.hlp", AWM_ALL, makeHelpCallback("dist.hlp"));
1639
1640    // ------------------
1641    //      left side
1642
1643    aws->at("which_species");
1644    aws->create_option_menu(AWAR_DIST_WHICH_SPECIES);
1645    aws->insert_option("all", "a", "all");
1646    aws->insert_default_option("marked",   "m", "marked");
1647    aws->update_option_menu();
1648
1649    aws->at("which_alignment");
1650    awt_create_ALI_selection_list(gb_main, (AW_window *)aws, AWAR_DIST_ALIGNMENT, "*=");
1651
1652    // filter & weights
1653
1654    AW_awar *awar_dist_alignment    = aws->get_root()->awar_string(AWAR_DIST_ALIGNMENT);
1655    WeightedFilter *weighted_filter = // do NOT free (bound to callbacks)
1656        new WeightedFilter(gb_main, aws->get_root(), AWAR_DIST_FILTER_NAME, AWAR_DIST_COLUMN_STAT_NAME, awar_dist_alignment);
1657
1658    aws->at("filter_select");
1659    aws->callback(makeCreateWindowCallback(awt_create_select_filter_win, weighted_filter->get_adfiltercbstruct()));
1660    aws->create_button("SELECT_FILTER", AWAR_DIST_FILTER_NAME);
1661
1662    aws->at("weights_select");
1663    aws->sens_mask(AWM_EXP);
1664    aws->callback(makeCreateWindowCallback(COLSTAT_create_selection_window, weighted_filter->get_column_stat()));
1665    aws->create_button("SELECT_COL_STAT", AWAR_DIST_COLUMN_STAT_NAME);
1666    aws->sens_mask(AWM_ALL);
1667
1668    aws->at("which_cancel");
1669    aws->create_input_field(AWAR_DIST_CANCEL_CHARS, 12);
1670
1671    aws->at("cancel_select");
1672    aws->callback(awt_create_select_cancel_window);
1673    aws->create_button("SELECT_CANCEL_CHARS", "Info", "C");
1674
1675    aws->at("change_matrix");
1676    aws->callback(create_dna_matrix_window);
1677    aws->create_button("EDIT_MATRIX", "Edit Matrix");
1678
1679    aws->at("enable");
1680    aws->create_toggle(AWAR_DIST_MATRIX_DNA_ENABLED);
1681
1682    aws->at("which_correction");
1683    aws->create_option_menu(AWAR_DIST_CORR_TRANS);
1684    aws->insert_option("none",                    "n", (int)DI_TRANSFORMATION_NONE);
1685    aws->insert_option("similarity",              "n", (int)DI_TRANSFORMATION_SIMILARITY);
1686    aws->insert_option("jukes-cantor (dna)",      "c", (int)DI_TRANSFORMATION_JUKES_CANTOR);
1687    aws->insert_option("felsenstein (dna)",       "f", (int)DI_TRANSFORMATION_FELSENSTEIN);
1688    aws->insert_option("olsen (dna)",             "o", (int)DI_TRANSFORMATION_OLSEN);
1689    aws->insert_option("felsenstein/voigt (exp)", "1", (int)DI_TRANSFORMATION_FELSENSTEIN_VOIGT);
1690    aws->insert_option("olsen/voigt (exp)",       "2", (int)DI_TRANSFORMATION_OLSEN_VOIGT);
1691    aws->insert_option("kimura (pro)",            "k", (int)DI_TRANSFORMATION_KIMURA);
1692    aws->insert_option("PAM (protein)",           "c", (int)DI_TRANSFORMATION_PAM);
1693    aws->insert_option("Cat. Hall(exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_HALL);
1694    aws->insert_option("Cat. Barker(exp)",        "c", (int)DI_TRANSFORMATION_CATEGORIES_BARKER);
1695    aws->insert_option("Cat.Chem (exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_CHEMICAL);
1696    aws->insert_option("from selected tree",      "t", (int)DI_TRANSFORMATION_FROM_TREE);
1697    aws->insert_default_option("unknown",         "u", (int)DI_TRANSFORMATION_NONE);
1698
1699    aws->update_option_menu();
1700
1701    aws->at("autodetect");   // auto
1702    aws->callback(makeWindowCallback(di_autodetect_callback, gb_main));
1703    aws->sens_mask(AWM_EXP);
1704    aws->create_button("AUTODETECT_CORRECTION", "AUTODETECT", "A");
1705    aws->sens_mask(AWM_ALL);
1706
1707    // -------------------
1708    //      right side
1709
1710
1711    aws->at("mark_distance");
1712    aws->callback(makeWindowCallback(di_mark_by_distance, weighted_filter));
1713    aws->create_autosize_button("MARK_BY_DIST", "Mark all species");
1714
1715    aws->at("mark_lower");
1716    aws->create_input_field(AWAR_DIST_MIN_DIST, 5);
1717
1718    aws->at("mark_upper");
1719    aws->create_input_field(AWAR_DIST_MAX_DIST, 5);
1720
1721    // -----------------
1722
1723    // tree selection
1724
1725    aws->at("tree_list");
1726    awt_create_TREE_selection_list(gb_main, aws, AWAR_DIST_TREE_CURR_NAME);
1727
1728    aws->at("detect_clusters");
1729    aws->callback(makeCreateWindowCallback(DI_create_cluster_detection_window, weighted_filter));
1730    aws->create_autosize_button("DETECT_CLUSTERS", "Detect homogenous clusters in tree", "D");
1731
1732    // matrix calculation
1733
1734    aws->button_length(18);
1735
1736    aws->at("calculate");
1737    aws->callback(makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1738    aws->create_button("CALC_FULL_MATRIX", "Calculate\nFull Matrix", "F");
1739
1740    aws->at("compress");
1741    aws->callback(makeWindowCallback(di_calculate_compressed_matrix_cb, weighted_filter));
1742    aws->create_button("CALC_COMPRESSED_MATRIX", "Calculate\nCompressed Matrix", "C");
1743
1744    recalculate_matrix_cb = new BoundWindowCallback(aws, makeWindowCallback(di_calculate_full_matrix_cb, weighted_filter));
1745
1746    aws->button_length(13);
1747
1748    {
1749        static save_matrix_params sparams;
1750
1751        sparams.awar_base       = AWAR_DIST_SAVE_MATRIX_BASE;
1752        sparams.weighted_filter = weighted_filter;
1753
1754        aws->at("save_matrix");
1755        aws->callback(makeCreateWindowCallback(DI_create_save_matrix_window, &sparams));
1756        aws->create_button("SAVE_MATRIX", "Save matrix", "M");
1757
1758        aws->at("view_matrix");
1759        aws->callback(makeWindowCallback(di_view_matrix_cb, &sparams));
1760        aws->create_button("VIEW_MATRIX", "View matrix", "V");
1761    }
1762
1763    aws->button_length(22);
1764    aws->at("use_compr_tree");
1765    aws->callback(di_define_compression_tree_name_cb);
1766    aws->create_button("USE_COMPRESSION_TREE", "Use to compress", "");
1767    aws->at("use_sort_tree");
1768    aws->callback(di_define_sort_tree_name_cb);
1769    aws->create_button("USE_SORT_TREE", "Use to sort", "");
1770
1771    aws->at("compr_tree_name"); aws->create_input_field(AWAR_DIST_TREE_COMP_NAME, 12);
1772    aws->at("sort_tree_name");  aws->create_input_field(AWAR_DIST_TREE_SORT_NAME, 12);
1773
1774    // tree calculation
1775
1776    aws->button_length(18);
1777
1778    aws->at("t_calculate");
1779    aws->callback(makeWindowCallback(di_calculate_tree_cb, weighted_filter, false));
1780    aws->create_button("CALC_TREE", "Calculate \ntree", "C");
1781
1782    aws->at("bootstrap");
1783    aws->callback(makeWindowCallback(di_calculate_tree_cb, weighted_filter, true));
1784    aws->create_button("CALC_BOOTSTRAP_TREE", "Calculate \nbootstrap tree");
1785
1786    recalculate_tree_cb = new BoundWindowCallback(aws, makeWindowCallback(di_calculate_tree_cb, weighted_filter, false));
1787
1788    aws->button_length(22);
1789    aws->at("use_existing");
1790    aws->callback(di_define_save_tree_name_cb);
1791    aws->create_button("USE_NAME", "Use as new tree name", "");
1792
1793    aws->at("calc_tree_name");
1794    aws->create_input_field(AWAR_DIST_TREE_STD_NAME, 12);
1795
1796    aws->at("bcount");
1797    aws->create_input_field(AWAR_DIST_BOOTSTRAP_COUNT, 7);
1798
1799    {
1800        aws->sens_mask(AWM_EXP);
1801
1802        aws->at("auto_calc_tree");
1803        aws->label("Auto calculate tree");
1804        aws->create_toggle(AWAR_DIST_MATRIX_AUTO_CALC_TREE);
1805
1806        aws->at("auto_recalc");
1807        aws->label("Auto recalculate");
1808        aws->create_toggle(AWAR_DIST_MATRIX_AUTO_RECALC);
1809
1810        aws->sens_mask(AWM_ALL);
1811    }
1812
1813    bool disable_autocalc = !ARB_in_expert_mode(aw_root);
1814    if (disable_autocalc) {
1815        aw_root->awar(AWAR_DIST_MATRIX_AUTO_RECALC)->write_int(0);
1816        aw_root->awar(AWAR_DIST_MATRIX_AUTO_CALC_TREE)->write_int(0);
1817    }
1818
1819    GB_pop_transaction(gb_main);
1820    return aws;
1821}
1822
1823// --------------------------------------------------------------------------------
1824
1825#ifdef UNIT_TESTS
1826#include <arb_diff.h>
1827#include <arb_file.h>
1828
1829#ifndef TEST_UNIT_H
1830#include <test_unit.h>
1831#endif
1832
1833class DIST_testenv : virtual Noncopyable {
1834    GB_shell  shell;
1835    GBDATA   *gb_main;
1836    AliView  *ali_view;
1837
1838public:
1839    DIST_testenv(const char *dbname, const char *aliName)
1840        : ali_view(NULp)
1841    {
1842        gb_main = GB_open(dbname, "r");
1843        TEST_REJECT_NULL(gb_main);
1844
1845        GB_transaction ta(gb_main);
1846
1847        int aliLength = GBT_get_alignment_len(gb_main, aliName);
1848        TEST_REJECT(aliLength<=0);
1849
1850        AP_filter filter(aliLength);
1851        if (!filter.is_invalid()) {
1852            AP_weights weights(&filter);
1853            ali_view = new AliView(gb_main, filter, weights, aliName);
1854        }
1855    }
1856    ~DIST_testenv() {
1857        delete ali_view;
1858        GB_close(gb_main);
1859    }
1860
1861    const AliView& aliview() const { return *ali_view; }
1862    GBDATA *gbmain() const { return gb_main; }
1863};
1864
1865void TEST_matrix() {
1866    for (int iat = GB_AT_RNA; iat<=GB_AT_AA; ++iat) {
1867        GB_alignment_type at = GB_alignment_type(iat);
1868        // ---------------
1869        //      setup
1870        //                               GB_AT_RNA         GB_AT_DNA           GB_AT_AA
1871        const char *db_name[]=  { NULp, "TEST_trees.arb", "TEST_realign.arb", "TEST_realign.arb", };
1872        const char *ali_name[]= { NULp, "ali_5s",         "ali_dna",          "ali_pro",          };
1873
1874        TEST_ANNOTATE(GBS_global_string("ali_name=%s", ali_name[at]));
1875        DIST_testenv env(db_name[at], ali_name[at]);
1876
1877        DI_MATRIX matrix(env.aliview());
1878        MatrixOrder order(env.gbmain(), "tree_abc"); // no such tree!
1879        TEST_EXPECT_NO_ERROR(matrix.load(DI_LOAD_MARKED, order, true, NULp));
1880
1881        // -------------------------------
1882        //      detect_transformation
1883        DI_TRANSFORMATION detected_trans;
1884        {
1885            string msg;
1886
1887            detected_trans = matrix.detect_transformation(msg);
1888            DI_TRANSFORMATION expected = DI_TRANSFORMATION_NONE_DETECTED;
1889            switch (at) {
1890                case GB_AT_RNA: expected = DI_TRANSFORMATION_NONE;         break;
1891                case GB_AT_DNA: expected = DI_TRANSFORMATION_JUKES_CANTOR; break;
1892                case GB_AT_AA:  expected = DI_TRANSFORMATION_PAM;          break;
1893                case GB_AT_UNKNOWN: di_assert(0); break;
1894            }
1895            TEST_EXPECT_EQUAL(detected_trans, expected);
1896        }
1897
1898        // ------------------------------
1899        //      calculate the matrix
1900
1901        // @@@ does not test user-defined transformation-matrix!
1902        if (at == GB_AT_AA) {
1903            matrix.calculate_pro(detected_trans, NULp);
1904        }
1905        else {
1906            if (at == GB_AT_RNA) detected_trans = DI_TRANSFORMATION_FELSENSTEIN; // force calculate_overall_freqs
1907            matrix.calculate("", detected_trans, NULp, NULp);
1908        }
1909
1910        // -----------------------------------
1911        //      save in available formats
1912
1913        for (DI_SAVE_TYPE saveType = DI_SAVE_PHYLIP_COMP; saveType<=DI_SAVE_TABBED; saveType = DI_SAVE_TYPE(saveType+1)) {
1914            const char *savename = "distance/matrix.out";
1915            matrix.save(savename, saveType);
1916
1917            const char *suffixAT[] = { NULp, "rna", "dna", "pro" };
1918            const char *suffixST[] = { "phylipComp", "readable", "tabbed" };
1919            char *expected = GBS_global_string_copy("distance/matrix.%s.%s.expected", suffixAT[at], suffixST[saveType]);
1920
1921// #define TEST_AUTO_UPDATE // uncomment to auto-update expected matrices
1922
1923#if defined(TEST_AUTO_UPDATE)
1924            TEST_COPY_FILE(savename, expected);
1925#else
1926            TEST_EXPECT_TEXTFILES_EQUAL(savename, expected);
1927#endif // TEST_AUTO_UPDATE
1928            TEST_EXPECT_ZERO_OR_SHOW_ERRNO(GB_unlink(savename));
1929
1930            free(expected);
1931        }
1932    }
1933}
1934
1935#endif // UNIT_TESTS
1936
1937// --------------------------------------------------------------------------------
1938
Note: See TracBrowser for help on using the repository browser.