source: branches/nameserver/DIST/DI_matr.cxx

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