source: tags/testbuild/DIST/DI_matr.cxx

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