source: branches/stable/DIST/DI_matr.cxx

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