source: branches/port5/DIST/DI_matr.cxx

Last change on this file was 6143, checked in by westram, 16 years ago
  • backport [6141] (parts changing code, but only strings and comments)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 59.0 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <memory.h>
4// #include <malloc.h>
5#include <string.h>
6
7#include <math.h>
8#include <time.h>
9#include <limits.h>
10
11#include <arbdb.h>
12#include <arbdbt.h>
13
14#include <aw_root.hxx>
15#include <aw_device.hxx>
16#include <aw_window.hxx>
17#include <aw_preset.hxx>
18#include <aw_global.hxx>
19#include <aw_awars.hxx>
20
21#include <awt.hxx>
22#include <awt_tree.hxx>
23#include <awt_sel_boxes.hxx>
24#include <awt_macro.hxx>
25#include <awt_csp.hxx>
26
27#include "dist.hxx"
28#include <BI_helix.hxx>
29#include <CT_ctree.hxx>
30#include "di_matr.hxx"
31#include "di_protdist.hxx"
32#include "di_view_matrix.hxx"
33
34#define di_assert(cond) arb_assert(cond)
35
36// --------------------------------------------------------------------------------
37
38#define AWAR_DIST_FILTER_PREFIX      AWAR_DIST_PREFIX "filter/"
39#define AWAR_DIST_COLUMN_STAT_PREFIX AWAR_DIST_PREFIX "col_stat/"
40#define AWAR_DIST_TREE_PREFIX        AWAR_DIST_PREFIX "tree/"
41
42#define AWAR_DIST_MATRIX_DNA_BASE    AWAR_DIST_PREFIX "dna/matrix"
43#define AWAR_DIST_MATRIX_DNA_ENABLED "tmp/" AWAR_DIST_MATRIX_DNA_BASE "/enable"
44
45#define AWAR_DIST_WHICH_SPECIES   AWAR_DIST_PREFIX "which_species"
46#define AWAR_DIST_ALIGNMENT       AWAR_DIST_PREFIX "alignment"
47#define AWAR_DIST_BOOTSTRAP_COUNT AWAR_DIST_PREFIX "bootstrap/count"
48#define AWAR_DIST_CANCEL_CHARS    AWAR_DIST_PREFIX "cancel/chars"
49
50#define AWAR_DIST_FILTER_ALIGNMENT AWAR_DIST_FILTER_PREFIX "alignment"
51#define AWAR_DIST_FILTER_NAME      AWAR_DIST_FILTER_PREFIX "name"
52#define AWAR_DIST_FILTER_FILTER    AWAR_DIST_FILTER_PREFIX "filter"
53#define AWAR_DIST_FILTER_SIMPLIFY  AWAR_DIST_FILTER_PREFIX "simplify"
54
55#define AWAR_DIST_COLUMN_STAT_ALIGNMENT AWAR_DIST_COLUMN_STAT_PREFIX "alignment"
56#define AWAR_DIST_COLUMN_STAT_NAME      AWAR_DIST_COLUMN_STAT_PREFIX "name"
57
58#define AWAR_DIST_TREE_STD_NAME  AWAR_DIST_TREE_PREFIX "tree_name"
59#define AWAR_DIST_TREE_CURR_NAME "tmp/" AWAR_DIST_TREE_PREFIX "curr_tree_name"
60#define AWAR_DIST_TREE_SORT_NAME "tmp/" AWAR_DIST_TREE_PREFIX "sort_tree_name"
61#define AWAR_DIST_TREE_COMP_NAME "tmp/" AWAR_DIST_TREE_PREFIX "compr_tree_name"
62
63#define AWAR_DIST_SAVE_MATRIX_TYPE     AWAR_DIST_SAVE_MATRIX_BASE "/type"
64#define AWAR_DIST_SAVE_MATRIX_FILENAME AWAR_DIST_SAVE_MATRIX_BASE "/file_name"
65
66// --------------------------------------------------------------------------------
67
68DI_MATRIX   *DI_MATRIX::ROOT = 0;
69DI_dmatrix *di_dmatrix     = 0;
70AWT_csp    *global_csp     = 0;
71
72static adfiltercbstruct *dist_filter_cl = 0;
73
74AP_matrix DI_dna_matrix(AP_MAX);
75
76static void delete_matrix_cb(AW_root *dummy)
77{
78    AWUSE(dummy);
79    delete DI_MATRIX::ROOT;
80    DI_MATRIX::ROOT = 0;
81    if (di_dmatrix){
82        di_dmatrix->resized();
83        di_dmatrix->display(false);
84    }
85}
86
87static AW_window *create_dna_matrix_window(AW_root *aw_root){
88    AW_window_simple *aws = 0;
89    aws = new AW_window_simple();
90    aws->init( aw_root, "SET_DNA_MATRIX","SET MATRIX");
91    aws->auto_increment(50,50);
92    aws->button_length(10);
93    aws->callback(AW_POPDOWN);
94    aws->create_button("CLOSE", "CLOSE");
95
96    aws->callback(AW_POPUP_HELP,(AW_CL)"user_matrix.hlp");
97    aws->create_button("HELP", "HELP");
98
99    aws->at_newline();
100
101    DI_dna_matrix.create_input_fields(aws,AWAR_DIST_MATRIX_DNA_BASE);
102    aws->window_fit();
103    return aws;
104}
105
106void DI_create_matrix_variables(AW_root *aw_root, AW_default def)
107{
108    GB_transaction dummy(GLOBAL_gb_main);
109    DI_dna_matrix.set_description("","");
110    DI_dna_matrix.x_description[AP_A] = strdup("A");
111    DI_dna_matrix.x_description[AP_C] = strdup("C");
112    DI_dna_matrix.x_description[AP_G] = strdup("G");
113    DI_dna_matrix.x_description[AP_T] = strdup("TU");
114    DI_dna_matrix.x_description[AP_S] = strdup("GAP");
115
116    DI_dna_matrix.y_description[AP_A] = strdup("A");
117    DI_dna_matrix.y_description[AP_C] = strdup("C");
118    DI_dna_matrix.y_description[AP_G] = strdup("G");
119    DI_dna_matrix.y_description[AP_T] = strdup("TU");
120    DI_dna_matrix.y_description[AP_S] = strdup("GAP");
121
122    DI_dna_matrix.create_awars(aw_root,AWAR_DIST_MATRIX_DNA_BASE);
123
124    aw_root->awar_int(AWAR_DIST_MATRIX_DNA_ENABLED, 0)->add_callback(delete_matrix_cb); // user matrix disabled by default
125    {
126        GBDATA *gbd = GB_search((GBDATA *)aw_root->application_database, AWAR_DIST_MATRIX_DNA_BASE, GB_FIND);
127        GB_add_callback(gbd, GB_CB_CHANGED, (GB_CB)delete_matrix_cb, 0);
128    }
129
130    aw_root->awar_string("tmp/dummy_string", "0", def);
131
132    aw_root->awar_string(AWAR_DIST_WHICH_SPECIES, "marked", def)->add_callback(delete_matrix_cb);
133    {
134        char *default_ali = GBT_get_default_alignment(GLOBAL_gb_main);
135        aw_root->awar_string( AWAR_DIST_ALIGNMENT,"",def)      ->add_callback(delete_matrix_cb)->write_string(default_ali);
136        free(default_ali);
137    }
138    aw_root->awar_string(AWAR_DIST_FILTER_ALIGNMENT, "none", def);
139    aw_root->awar_string(AWAR_DIST_FILTER_NAME,      "none", def);
140    aw_root->awar_string(AWAR_DIST_FILTER_FILTER,    "",     def)->add_callback(delete_matrix_cb);
141    aw_root->awar_int   (AWAR_DIST_FILTER_SIMPLIFY,  0,      def)->add_callback(delete_matrix_cb);
142
143    aw_root->awar_string(AWAR_DIST_COLUMN_STAT_NAME,      "none", def);
144    aw_root->awar_string(AWAR_DIST_COLUMN_STAT_ALIGNMENT, "none", def);
145
146    aw_root->awar_string(AWAR_DIST_CANCEL_CHARS, ".", def)->add_callback(delete_matrix_cb);   
147    aw_root->awar_int(AWAR_DIST_CORR_TRANS, (int)DI_TRANSFORMATION_SIMILARITY, def)->add_callback(delete_matrix_cb);
148
149    aw_root->awar(AWAR_DIST_FILTER_ALIGNMENT)     ->map(AWAR_DIST_ALIGNMENT);
150    aw_root->awar(AWAR_DIST_COLUMN_STAT_ALIGNMENT)->map(AWAR_DIST_ALIGNMENT);
151
152    aw_create_selection_box_awars(aw_root, AWAR_DIST_SAVE_MATRIX_BASE, ".", "", "infile", def);
153    aw_root->awar_int(AWAR_DIST_SAVE_MATRIX_TYPE, 0, def);
154
155    aw_root->awar_string(AWAR_DIST_TREE_STD_NAME,  "tree_nj", def);
156    {
157        char *currentTree = aw_root->awar_string(AWAR_TREE, "", GLOBAL_gb_main)->read_string();
158        aw_root->awar_string(AWAR_DIST_TREE_CURR_NAME, currentTree, def);
159        aw_root->awar_string(AWAR_DIST_TREE_SORT_NAME, currentTree, def)->add_callback(delete_matrix_cb);
160        free(currentTree);
161    }
162    aw_root->awar_string(AWAR_DIST_TREE_COMP_NAME, "?????", def)->add_callback(delete_matrix_cb);
163
164    aw_root->awar_int(AWAR_DIST_BOOTSTRAP_COUNT, 1000, def);
165
166    aw_root->awar_float(AWAR_DIST_MIN_DIST, 0.0);
167    aw_root->awar_float(AWAR_DIST_MAX_DIST, 0.0);
168   
169    aw_root->awar_string(AWAR_SPECIES_NAME, "", GLOBAL_gb_main);
170
171    GB_push_transaction(GLOBAL_gb_main);
172
173    GBDATA *gb_species_data = GB_search(GLOBAL_gb_main,"species_data",GB_CREATE_CONTAINER);
174    GB_add_callback(gb_species_data,GB_CB_CHANGED,(GB_CB)delete_matrix_cb,0);
175
176    GB_pop_transaction(GLOBAL_gb_main);
177
178}
179
180DI_ENTRY::DI_ENTRY(GBDATA *gbd,class DI_MATRIX *phmatri) {
181    memset((char *)this,0,sizeof(DI_ENTRY));
182    phmatrix = phmatri;
183
184    GBDATA *gb_ali = GB_entry(gbd,phmatrix->use);
185    if (gb_ali) {
186        GBDATA *gb_data = GB_entry(gb_ali,"data");
187        if (gb_data) {
188            const char *seq = GB_read_char_pntr(gb_data);
189           
190            if (phmatrix->is_AA) {
191                sequence = sequence_protein = new AP_sequence_simple_protein(phmatrix->tree_root);
192            }
193            else {
194                sequence = sequence_parsimony = new AP_sequence_parsimony(phmatrix->tree_root);
195            }
196            sequence->set(seq);
197
198            name      = GBT_read_string(gbd, "name");
199            full_name = GBT_read_string(gbd, "full_name");
200        }
201    }
202}
203
204DI_ENTRY::DI_ENTRY(char *namei,class DI_MATRIX *phmatri)
205{
206    memset((char *)this,0,sizeof(DI_ENTRY));
207    phmatrix = phmatri;
208    this->name = strdup(namei);
209}
210
211DI_ENTRY::~DI_ENTRY(void)
212{
213    delete sequence_protein;
214    delete sequence_parsimony;
215    free(name);
216    free(full_name);
217
218}
219
220DI_MATRIX::DI_MATRIX(GBDATA *gb_maini,AW_root *awr)
221{
222    memset((char *)this,0,sizeof(DI_MATRIX));
223    aw_root = awr;
224    gb_main = gb_maini;
225}
226
227char *DI_MATRIX::unload(void)
228{
229    freeset(use, 0);
230    long i;
231    for (i=0;i<nentries;i++){
232        delete entries[i];
233    }
234    delete tree_root;
235    freeset(entries, 0);
236    nentries = 0;
237    return 0;
238}
239
240DI_MATRIX::~DI_MATRIX(void)
241{
242    unload();
243    delete matrix;
244}
245
246extern "C" int qsort_strcmp(const void *str1ptr, const void *str2ptr) {
247    GB_CSTR s1 = *(GB_CSTR*)str1ptr;
248    GB_CSTR s2 = *(GB_CSTR*)str2ptr;
249
250    return strcmp(s1, s2);
251}
252
253char *DI_MATRIX::load(char *usei, AP_filter *filter, AP_weights *weights, LoadWhat what, GB_CSTR sort_tree_name, bool show_warnings, GBDATA **species_list) {
254    delete tree_root;
255    tree_root = new AP_tree_root(gb_main,0,0);
256    tree_root->filter = filter;
257    tree_root->weights = weights;
258    this->use = strdup(usei);
259    this->gb_main = GLOBAL_gb_main;
260    GB_push_transaction(gb_main);
261
262    seq_len          = GBT_get_alignment_len(gb_main,use);
263    is_AA            = GBT_is_alignment_protein(gb_main,use);
264    gb_species_data  = GB_search(gb_main,"species_data",GB_CREATE_CONTAINER);
265    entries_mem_size = 1000;
266
267    entries = (DI_ENTRY **)calloc(sizeof(DI_ENTRY*),(size_t)entries_mem_size);
268
269    nentries = 0;
270    GBDATA *gb_species;
271    DI_ENTRY *phentry;
272
273    int tree_size;
274    GBT_TREE *sort_tree = GBT_read_tree_and_size(gb_main, sort_tree_name, sizeof(GBT_TREE), &tree_size);
275    tree_size++;
276    GB_CSTR *species_in_sort_tree = 0;
277
278    int no_of_species;
279    switch (what) {
280        case DI_LOAD_ALL:
281            no_of_species = GBT_get_species_count(gb_main);
282            break;
283        case DI_LOAD_MARKED:
284            no_of_species = GBT_count_marked_species(gb_main);
285            break;
286        case DI_LOAD_LIST:
287            di_assert(species_list);
288            for (no_of_species = 0; species_list[no_of_species]; ++no_of_species) ;
289            break;
290    }
291    int in_tree_species         = 0;
292    int out_tree_species        = no_of_species;
293    int unknown_species_in_tree = 0;
294
295    if (!sort_tree) {
296        if (show_warnings) {
297            aw_message("No valid tree given to sort matrix (using default database order)");
298        }
299    }
300    else {
301        species_in_sort_tree = GBT_get_species_names_of_tree(sort_tree);
302
303        {
304            GB_HASHI *shallLoad = 0;
305            if (what == DI_LOAD_LIST) {
306                shallLoad = GBS_create_hashi(no_of_species);
307                for (int i = 0; species_list[i]; ++i) {
308                    GBS_write_hashi(shallLoad, (long)species_list[i], 1);
309                }
310            }
311
312            for (int i=0; i<tree_size; i++) { // read all marked species in tree
313                gb_species = GBT_find_species_rel_species_data(gb_species_data, species_in_sort_tree[i]);
314                if (!gb_species) {
315                    if (show_warnings) {
316                        aw_message(GB_export_errorf("Species '%s' found in tree '%s' but NOT in database.",
317                                                    species_in_sort_tree[i], sort_tree_name));
318                    }
319                    unknown_species_in_tree++;
320                    continue;
321                }
322
323                if ((what == DI_LOAD_ALL)                                ||
324                    (what == DI_LOAD_MARKED && GB_read_flag(gb_species)) ||
325                    (what == DI_LOAD_LIST && GBS_read_hashi(shallLoad, (long)gb_species)))
326                {
327                    in_tree_species++;
328                    phentry = new DI_ENTRY(gb_species,this);
329                    if (phentry->sequence) {    // a species found
330                        entries[nentries++] = phentry;
331                        if (nentries == entries_mem_size) {
332                            entries_mem_size +=1000;
333                            entries = (DI_ENTRY **)realloc((MALLOC_T)entries,(size_t)(sizeof(DI_ENTRY*)*entries_mem_size));
334                        }
335                    }else{
336                        delete phentry;
337                    }
338                }
339            }
340
341            if (shallLoad) GBS_free_hashi(shallLoad);
342        }
343        out_tree_species = no_of_species-in_tree_species; // # of species not loaded yet (because not in tree)
344        gb_assert(out_tree_species>=0);
345        if (out_tree_species) {
346            qsort(species_in_sort_tree, tree_size, sizeof(*species_in_sort_tree), qsort_strcmp); // sort species names
347#if defined(DEBUG) && 0
348            printf("Sorted tree:\n");
349            for (int x=0; x<tree_size; x++) {
350                printf("- %s\n", species_in_sort_tree[x]);
351            }
352#endif
353        }
354    }
355
356    if (unknown_species_in_tree && show_warnings) {
357        aw_message(GB_export_errorf("We found %i species in tree '%s' which are not in database.\n"
358                                    "This does not affect the current calculation, but you should think about it.",
359                                    unknown_species_in_tree, sort_tree_name));
360    }
361
362    gb_assert(out_tree_species>=0);
363    if (out_tree_species) { // there are species outside the tree (or no valid sort tree selected) => load all marked species not loaded yet
364        int count   = 0;
365        int listIdx = 0;
366
367        switch (what) {
368            case DI_LOAD_ALL:    gb_species = GBT_first_species_rel_species_data(gb_species_data); break;
369            case DI_LOAD_MARKED: gb_species = GBT_first_marked_species_rel_species_data(gb_species_data); break;
370            case DI_LOAD_LIST:   gb_species = species_list[listIdx++]; break;
371        }
372
373        while (gb_species) {
374            int load_species = 0;
375
376            if (in_tree_species) { // test if species already loaded
377                const char *species_name = GBT_read_name(gb_species);
378                if (bsearch(&species_name, species_in_sort_tree, tree_size, sizeof(*species_in_sort_tree), qsort_strcmp)==0) {
379                    load_species = 1;
380                }
381            }
382            else {
383                load_species = 1;
384            }
385
386            if (load_species) {
387                count++;
388                phentry = new DI_ENTRY(gb_species,this);
389                if (phentry->sequence) {    // a species found
390                    entries[nentries++] = phentry;
391                    if (nentries == entries_mem_size) {
392                        entries_mem_size +=1000;
393                        entries = (DI_ENTRY **)realloc((MALLOC_T)entries,(size_t)(sizeof(DI_ENTRY*)*entries_mem_size));
394                    }
395                }else{
396                    delete phentry;
397                }
398            }
399
400            switch (what) {
401                case DI_LOAD_ALL:    gb_species = GBT_next_species(gb_species); break;
402                case DI_LOAD_MARKED: gb_species = GBT_next_marked_species(gb_species); break;
403                case DI_LOAD_LIST:   gb_species = species_list[listIdx++]; break;
404            }
405        }
406    }
407
408    free(species_in_sort_tree);
409    if (sort_tree) GBT_delete_tree(sort_tree);
410
411    GB_pop_transaction(gb_main);
412
413    return 0;
414}
415
416
417/******************************************************************************************
418            Calculate the global rate_matrix
419******************************************************************************************/
420void DI_MATRIX::clear(DI_MUT_MATR &hits)
421{
422    int i,j;
423    for (i=0;i<AP_MAX;i++){
424        for (j=0;j<AP_MAX;j++){
425            hits[i][j] = 0;
426        }
427    }
428}
429
430void DI_MATRIX::make_sym(DI_MUT_MATR &hits)
431{
432    int i,j;
433    for (i=AP_A;i<AP_MAX;i*=2){
434        for (j=AP_A;j<=i;j*=2){
435            hits[i][j] = hits[j][i] = hits[i][j] + hits[j][i];
436        }
437    }
438}
439
440void DI_MATRIX::rate_write(DI_MUT_MATR &hits,FILE *out){
441    int i,j;
442    for (i=AP_A;i<AP_MAX;i*=2){
443        for (j=AP_A;j<AP_MAX;j*=2){
444            fprintf(out,"%5li ",hits[i][j]);
445        }
446        fprintf(out,"\n");
447    }
448}
449
450long *DI_MATRIX::create_helix_filter(BI_helix *helix,AP_filter *filter){
451    long *result = (long *)calloc(sizeof(long),(size_t)filter->real_len);
452    long *temp = (long *)calloc(sizeof(long),(size_t)filter->real_len);
453    int i,count;
454    count = 0;
455    for (i=0;i<filter->filter_len;i++) {
456        if (filter->filter_mask[i]) {
457            temp[i] = count;
458            if (helix->pairtype(i)== HELIX_PAIR) {
459                result[count] = helix->opposite_position(i);
460            }else{
461                result[count] = -1;
462            }
463            count++;
464        }
465    }
466    while (--count >=0){
467        if (result[count] >= 0) {
468            result[count] = temp[result[count]];
469        }
470    }
471    free((char *)temp);
472    return result;
473}
474
475GB_ERROR DI_MATRIX::calculate_rates(DI_MUT_MATR &hrates,DI_MUT_MATR &nrates,DI_MUT_MATR &pairs,long *filter){
476
477    if (nentries<=1) {
478        return "There are no species selected";
479    }
480    long   row,col,pos,s_len;
481    unsigned char *seq1,*seq2;
482    s_len = tree_root->filter->real_len;
483    this->clear(hrates);
484    this->clear(nrates);
485    this->clear(pairs);
486    for (row = 0; row<nentries;row++){
487        double gauge = (double)row/(double)nentries;
488        if (aw_status(gauge*gauge)) return "Aborted";
489        for (col=0; col<row; col++) {
490            seq1 = (unsigned char *)entries[row]->sequence_parsimony->sequence;
491            seq2 = (unsigned char *)entries[col]->sequence_parsimony->sequence;
492            for (pos = 0; pos < s_len; pos++){
493                if(filter[pos]>=0) {
494                    hrates[*seq1][*seq2]++;
495                }else{
496                    nrates[*seq1][*seq2]++;
497                }
498                seq1++; seq2++;
499            }
500        }
501    }
502    for (row = 0; row<nentries;row++){
503        seq1 = (unsigned char *)entries[row]->sequence_parsimony->sequence;
504        for (pos = 0; pos < s_len; pos++){
505            if(filter[pos]>=0) {
506                pairs[seq1[pos]][seq1[filter[pos]]]++;
507            }
508        }
509    }
510    return 0;
511}
512
513/******************************************************************************************
514    Some test functions to check correlated base correction
515******************************************************************************************/
516
517
518GB_ERROR DI_MATRIX::haeschoe(const char *path)
519{
520    DI_MUT_MATR temp,temp2,pairs;
521    static BI_helix *helix = 0;
522    if (!helix) helix = new BI_helix();
523    const char *error = helix->init(gb_main);
524    if (error) return error;
525
526    FILE *out = fopen(path,"w");
527    if (!out) return "Cannot open file";
528
529    aw_openstatus("Calculating Distance Matrix");
530    aw_status("Calculating global rate matrix");
531
532    fprintf(out,"Pairs in helical regions:\n");
533    long *filter = create_helix_filter(helix,tree_root->filter);
534    error = calculate_rates(temp,temp2,pairs,filter);
535    if (error){
536        aw_closestatus();
537        fclose(out);
538        return error;
539    }
540    rate_write(pairs,out);
541    make_sym(temp);make_sym(temp2);
542    fprintf(out,"\nRatematrix helical parts:\n");
543    rate_write(temp,out);
544    fprintf(out,"\nRatematrix non helical parts:\n");
545    rate_write(temp2,out);
546
547    long   row,col,pos,s_len;
548    char *seq1,*seq2;
549    s_len = tree_root->filter->real_len;
550    fprintf(out,"\nDistance matrix (Helixdist Helixlen Nonhelixdist Nonhelixlen):");
551    fprintf(out,"\n%li",nentries);
552    const int MAXDISTDEBUG = 1000;
553    double distdebug[MAXDISTDEBUG];
554    for (pos = 0;pos<MAXDISTDEBUG;pos++) distdebug[pos] = 0.0;
555    for (row = 0; row<nentries;row++){
556        if ( nentries > 100 || ((row&0xf) == 0)){
557            double gauge = (double)row/(double)nentries;
558            if (aw_status(gauge*gauge)) {
559                aw_closestatus();
560                return "Aborted";
561            }
562        }
563        fprintf (out,"\n%s  ",entries[row]->name);
564
565        for (col=0; col<row; col++) {
566            seq1 = entries[row]->sequence_parsimony->sequence;
567            seq2= entries[col]->sequence_parsimony->sequence;
568            this->clear(temp);
569            this->clear(temp2);
570            for (pos = 0; pos < s_len; pos++){
571                if(filter[pos]>=0) {
572                    temp[(unsigned char)*seq1][(unsigned char)*seq2]++;
573                }else{
574                    temp2[(unsigned char)*seq1][(unsigned char)*seq2]++;
575                }
576                seq1++; seq2++;
577            }
578            long hdist= 0, hall2 = 0;
579            int i,j;
580            for (i=AP_A;i<AP_MAX;i*=2){
581                for (j=AP_A;j<AP_MAX;j*=2){
582                    hall2 += temp[i][j];
583                    if (i!=j) hdist += temp[i][j];
584                }
585            }
586
587            long dist= 0, all = 0;
588            for (i=AP_A;i<=AP_T;i*=2){
589                for (j=AP_A;j<=AP_T;j*=2){
590                    all += temp2[i][j];
591                    if (i!=j) dist += temp2[i][j];
592                }
593            }
594            fprintf(out,"(%4li:%4li %4li:%4li) ",hdist,hall2, dist,all);
595            if (all>100){
596                distdebug[dist*MAXDISTDEBUG/all] = (double)hdist/(double)hall2;
597            }
598        }
599    }
600
601    fprintf (out,"\n");
602    fprintf (out,"\nhdist/dist:\n");
603
604    for (pos = 1;pos<MAXDISTDEBUG;pos++) {
605        if (distdebug[pos]) {
606            fprintf(out,"%4f %5f\n",(double)pos/(double)MAXDISTDEBUG,distdebug[pos]/(double)pos*(double)MAXDISTDEBUG);
607        }
608    }
609    aw_closestatus();
610    fclose(out);
611
612    return 0;
613}
614
615char *DI_MATRIX::calculate_overall_freqs(double rel_frequencies[AP_MAX], char *cancel)
616{
617    long  hits2[AP_MAX];
618    long  sum   = 0;
619    int   i,row;
620    char *seq1;
621    int   pos;
622    int   b;
623    long  s_len = tree_root->filter->real_len;
624
625    memset((char *) &hits2[0], 0, sizeof(hits2));
626    for (row = 0; row < nentries; row++) {
627        seq1 = entries[row]->sequence_parsimony->sequence;
628        for (pos = 0; pos < s_len; pos++) {
629            b = *(seq1++);
630            if (cancel[b]) continue;
631            hits2[b]++;
632        }
633    }
634    for (i = 0; i < AP_MAX; i++)    sum += hits2[i];
635    for (i = 0; i < AP_MAX; i++)    rel_frequencies[i] = hits2[i] / (double) sum;
636    return 0;
637}
638
639double DI_MATRIX::corr(double dist, double b, double & sigma){
640    const double eps = 0.01;
641    double ar = 1.0 - dist/b;
642    sigma = 1000.0;
643    if (ar< eps) return 3.0;
644    sigma = b/ar;
645    return - b * log(1-dist/b);
646}
647
648GB_ERROR DI_MATRIX::calculate(AW_root *awr, char *cancel, double /*alpha*/, DI_TRANSFORMATION transformation)
649{
650
651    if (transformation == DI_TRANSFORMATION_HAESCH) {
652        GB_ERROR error = haeschoe("outfile");
653        if (error) return error;
654        return "Your matrices have been written to 'outfile'\nSorry I can not make a tree";
655    }
656    int user_matrix_enabled = awr->awar(AWAR_DIST_MATRIX_DNA_ENABLED)->read_int();
657    if (user_matrix_enabled){   // set transformation Matrix
658        switch (transformation){
659            case DI_TRANSFORMATION_NONE:
660            case DI_TRANSFORMATION_SIMILARITY:
661            case DI_TRANSFORMATION_JUKES_CANTOR:
662                break;
663            default:
664                return GB_export_error("Sorry not implemented:\n"
665                                       "    This kind of distance correction does not support\n"
666                                       "    a user defined matrix");
667        }
668        DI_dna_matrix.read_awars(awr,AWAR_DIST_MATRIX_DNA_BASE);
669        DI_dna_matrix.normize();
670    }
671
672
673
674    matrix = new AP_smatrix(nentries);
675
676    long   s_len = tree_root->filter->real_len;
677    long   hits[AP_MAX][AP_MAX];
678    int    row;
679    int    col;
680    size_t i;
681
682    if (nentries<=1) {
683        return "There are no species selected";
684    }
685    memset(&cancel_columns[0],0,256);
686
687    for (i=0;i<strlen(cancel);i++){
688        int j = cancel[i];
689        j = AP_sequence_parsimony::table[j];
690        cancel_columns[j] = 1;
691    }
692    long    columns;
693    double b;
694    long frequencies[AP_MAX];
695    double rel_frequencies[AP_MAX];
696    double S_square=0;
697
698    switch(transformation){
699        case DI_TRANSFORMATION_FELSENSTEIN:
700            this->calculate_overall_freqs(rel_frequencies,cancel_columns);
701            S_square = 0.0;
702            for (i=0;i<AP_MAX; i++) S_square+= rel_frequencies[i]*rel_frequencies[i];
703            break;
704        default:    break;
705    };
706
707    for (row = 0; row<nentries;row++){
708        if ( nentries > 100 || ((row&0xf) == 0)){
709            double gauge = (double)row/(double)nentries;
710            if (aw_status(gauge*gauge)) {
711                aw_closestatus();
712                return "Aborted";
713            }
714        }
715        for (col=0; col<=row; col++) {
716            columns    = 0;
717            char *seq1 = entries[row]->sequence_parsimony->sequence;
718            char *seq2 = entries[col]->sequence_parsimony->sequence;
719            b          = 0.0;
720            switch(transformation){
721                case  DI_TRANSFORMATION_JUKES_CANTOR:
722                    b = 0.75;
723                case  DI_TRANSFORMATION_NONE:
724                case  DI_TRANSFORMATION_SIMILARITY:{
725                    double  dist = 0.0;
726                    if (user_matrix_enabled){
727                        memset((char *)hits,0,sizeof(long) * AP_MAX * AP_MAX);
728                        int pos;
729                        for (pos = s_len; pos >=0; pos--){
730                            hits[(unsigned char)*(seq1++)][(unsigned char)*(seq2++)]++;
731                        }
732                        int x,y;
733                        double diffsum = 0.0;
734                        double all_sum = 0.001;
735                        for (x = AP_A; x < AP_MAX; x*=2){
736                            for (y = AP_A; y < AP_MAX; y*=2){
737                                if (x==y){
738                                    all_sum += hits[x][y];
739                                }else{
740                                    diffsum += hits[x][y] * DI_dna_matrix.get(x,y);
741                                    all_sum += hits[x][y] * DI_dna_matrix.get(x,y);
742                                }
743                            }
744                        }
745                        dist = diffsum / all_sum;
746                    }else{
747                        int pos;
748                        int b1,b2;
749                        for (pos = s_len; pos >=0; pos--){
750                            b1 =*(seq1++);
751                            b2 =*(seq2++);
752                            if (cancel_columns[b1]) continue;
753                            if (cancel_columns[b2]) continue;
754                            columns++;
755                            if (b1&b2) continue;
756                            dist+=1.0;
757                        }
758                        if (columns== 0) columns = 1;
759                        dist /= columns;
760                    }
761                    if (transformation==DI_TRANSFORMATION_SIMILARITY){
762                        dist =  (1.0-dist);
763                    }else   if (b){
764                        double sigma;
765                        dist = this->corr(dist,b,sigma);
766                    }
767                    matrix->set(row,col,dist);
768                    break;
769                }
770                case DI_TRANSFORMATION_KIMURA:
771                case DI_TRANSFORMATION_OLSEN:
772                case DI_TRANSFORMATION_FELSENSTEIN:{
773                    int    pos;
774                    double dist = 0.0;
775                    long   N,P,Q,M;
776                    double p,q;
777                   
778                    memset((char *)hits,0,sizeof(long) * AP_MAX * AP_MAX);
779                    for (pos = s_len; pos >=0; pos--){
780                        hits[(unsigned char)*(seq1++)][(unsigned char)*(seq2++)]++;
781                    }
782                    switch(transformation){
783                        case DI_TRANSFORMATION_KIMURA:
784                            P = hits[AP_A][AP_G] +
785                                hits[AP_G][AP_A] +
786                                hits[AP_C][AP_T] +
787                                hits[AP_T][AP_C];
788                            Q = hits[AP_A][AP_C] +
789                                hits[AP_A][AP_T] +
790                                hits[AP_C][AP_A] +
791                                hits[AP_T][AP_A] +
792                                hits[AP_G][AP_C] +
793                                hits[AP_G][AP_T] +
794                                hits[AP_C][AP_G] +
795                                hits[AP_T][AP_G];
796                            M = hits[AP_A][AP_A] +
797                                hits[AP_C][AP_C] +
798                                hits[AP_G][AP_G] +
799                                hits[AP_T][AP_T];
800                            N = P+Q+M;
801                            if (N==0) N=1;
802                            p = (double)P/(double)N;
803                            q = (double)Q/(double)N;
804                            dist = - .5 * log(
805                                              (1.0-2.0*p-q)*sqrt(1.0-2.0*q)
806                                              );
807                            break;
808
809                        case DI_TRANSFORMATION_OLSEN:
810                        case DI_TRANSFORMATION_FELSENSTEIN:
811
812                            memset((char *)frequencies,0,
813                                   sizeof(long) * AP_MAX);
814
815                            N = 0;
816                            M = 0;
817
818                            for (i=0;i<AP_MAX;i++) {
819                                if (cancel_columns[i]) continue;
820                                unsigned int j;
821                                for (j=0;j<i;j++) {
822                                    if (cancel_columns[j]) continue;
823                                    frequencies[i] +=
824                                        hits[i][j]+
825                                        hits[j][i];
826                                }
827                                frequencies[i] += hits[i][i];
828                                N += frequencies[i];
829                                M += hits[i][i];;
830                            }
831                            if (N==0) N=1;
832                            if (transformation == DI_TRANSFORMATION_OLSEN){ // Calc sum square freq individually for each line
833                                S_square = 0.0;
834                                for (i=0;i<AP_MAX; i++) S_square+= frequencies[i]*frequencies[i];
835                                b = 1.0 - S_square/((double)N*(double)N);
836                            }else{
837                                b = 1.0 - S_square;
838                            }
839
840                            dist = ((double)(N-M)) / (double) N;
841                            double sigma;
842                            dist = this->corr(dist,b,sigma);
843                            break;
844
845                        default: return "Sorry: Transformation not implemented";
846                    }
847                    matrix->set(row,col,dist);
848                    break;
849                }
850                default:;
851            }   /* switch */
852        }
853    }
854    return 0;
855}
856
857GB_ERROR DI_MATRIX::calculate_pro(DI_TRANSFORMATION transformation){
858    di_cattype whichcat;
859    di_codetype whichcode = universal;
860
861    switch (transformation){
862        case DI_TRANSFORMATION_NONE:
863            whichcat = none;
864            break;
865        case DI_TRANSFORMATION_SIMILARITY:
866            whichcat = similarity;
867            break;
868        case DI_TRANSFORMATION_KIMURA:
869            whichcat = kimura;
870            break;
871        case DI_TRANSFORMATION_PAM:
872            whichcat = pam;
873            break;
874        case DI_TRANSFORMATION_CATEGORIES_HALL:
875            whichcat = hall;
876            break;
877        case DI_TRANSFORMATION_CATEGORIES_BARKER:
878            whichcat = george;
879            break;
880        case DI_TRANSFORMATION_CATEGORIES_CHEMICAL:
881            whichcat = chemical;
882            break;
883        default:
884            return "This correction is not available for protein data";
885    }
886    matrix = new AP_smatrix(nentries);
887
888    di_protdist prodist(whichcode, whichcat, nentries,entries, tree_root->filter->real_len,matrix);
889    return prodist.makedists();
890}
891
892static GB_ERROR di_calculate_matrix_cb(AW_window *aww, AW_CL bootstrap_flag, AW_CL cl_show_warnings) // returns "Aborted" if stopped by user
893{
894    bool show_warnings = (bool)cl_show_warnings;
895    GB_push_transaction(GLOBAL_gb_main);
896    if (DI_MATRIX::ROOT){
897        GB_pop_transaction(GLOBAL_gb_main);
898        return 0;
899    }
900    AW_root *aw_root = aww->get_root();
901    char    *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
902    long ali_len = GBT_get_alignment_len(GLOBAL_gb_main,use);
903    if (ali_len<=0) {
904        GB_pop_transaction(GLOBAL_gb_main);
905        delete use;
906        aw_message("Please select a valid alignment");
907        return "Error";
908    }
909    if (!bootstrap_flag){
910        aw_openstatus("Calculating Matrix");
911        aw_status("Read the database");
912    }
913
914    char *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
915
916    AP_filter *ap_filter = awt_get_filter(aw_root, dist_filter_cl);
917    if (bootstrap_flag){
918        ap_filter->enable_bootstrap();
919    }
920
921    AP_weights *ap_weights = new AP_weights;
922    ap_weights->init(ap_filter);
923
924    char *load_what      = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
925    char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
926
927    LoadWhat all_flag = (strcmp(load_what,"all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
928    GB_ERROR error    = 0;
929    {
930        DI_MATRIX *phm = new DI_MATRIX(GLOBAL_gb_main,aw_root);
931        phm->matrix_type = DI_MATRIX_FULL;
932        phm->load(use, ap_filter, ap_weights, all_flag, sort_tree_name, show_warnings, NULL);
933        free(sort_tree_name);
934        GB_pop_transaction(GLOBAL_gb_main);
935        if (aw_status()) {
936            phm->unload();
937            return "Aborted";
938        }
939
940        DI_TRANSFORMATION trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
941
942        if (phm->is_AA){
943            error = phm->calculate_pro(trans);
944        }else{
945            error = phm->calculate(aw_root,cancel,0.0,trans);
946        }
947        if (!bootstrap_flag){
948            aw_closestatus();
949        }
950        delete phm->ROOT;
951        if (error && strcmp(error,"Aborted")){
952            aw_message(error);
953            delete phm;
954            phm->ROOT = 0;
955        }
956        else {
957            phm->ROOT = phm;
958
959        }
960    }
961    free(load_what);
962
963    free(cancel);
964    delete ap_filter;
965    delete ap_weights;
966
967    free(use);
968    return error;
969}
970
971static void di_mark_by_distance(AW_window *aww) {
972    AW_root *aw_root    = aww->get_root();
973    double   lowerBound = aw_root->awar(AWAR_DIST_MIN_DIST)->read_float();
974    double   upperBound = aw_root->awar(AWAR_DIST_MAX_DIST)->read_float();
975
976    GB_ERROR error = 0;
977    if (lowerBound >= upperBound) {
978        error = GBS_global_string("Lower bound (%f) has to be smaller than upper bound (%f)", lowerBound, upperBound);
979    }
980    else if (lowerBound<0.0 || lowerBound > 1.0) {
981        error = GBS_global_string("Lower bound (%f) is not in allowed range [0.0 .. 1.0]", lowerBound);
982    }
983    else if (upperBound<0.0 || upperBound > 1.0) {
984        error = GBS_global_string("Upper bound (%f) is not in allowed range [0.0 .. 1.0]", upperBound);
985    }
986    else {
987        GB_transaction ta(GLOBAL_gb_main);
988
989        char   *selected    = aw_root->awar(AWAR_SPECIES_NAME)->read_string();
990        GBDATA *gb_selected = 0;
991
992        if (!selected[0]) {
993            error = "Please select a species";
994        }
995        else {
996            gb_selected = GBT_find_species(GLOBAL_gb_main, selected);
997            if (!gb_selected) {
998                error = GBS_global_string("Couldn't find species '%s'", selected);
999            }
1000            else {
1001                char              *use           = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1002                char              *cancel        = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1003                AP_filter         *ap_filter     = awt_get_filter(aw_root,dist_filter_cl);
1004                AP_weights        *ap_weights    = new AP_weights;
1005                AP_smatrix        *ap_ratematrix = 0;
1006                DI_TRANSFORMATION  trans         = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1007
1008                ap_weights->init(ap_filter);
1009
1010                if (!error) {
1011                    aw_openstatus("Mark species by distance");
1012                    aw_status("Calculating distances");
1013                    aw_status(0.0);
1014
1015                    DI_MATRIX *old_root = DI_MATRIX::ROOT;
1016                    DI_MATRIX::ROOT = 0;
1017
1018                    size_t speciesCount = GBT_get_species_count(GLOBAL_gb_main);
1019                    size_t count        = 0;
1020                    bool markedSelected = false;
1021
1022                    for (GBDATA *gb_species = GBT_first_species(GLOBAL_gb_main);
1023                         gb_species;
1024                         gb_species = GBT_next_species(gb_species))
1025                    {
1026                        DI_MATRIX *phm         = new DI_MATRIX(GLOBAL_gb_main, aw_root);
1027                        phm->matrix_type       = DI_MATRIX_FULL;
1028                        GBDATA *species_pair[] = { gb_selected, gb_species, NULL };
1029
1030                        phm->load(use, ap_filter, ap_weights, DI_LOAD_LIST, NULL, false, species_pair);
1031
1032                        if (phm->is_AA){
1033                            error = phm->calculate_pro(trans);
1034                        }
1035                        else {
1036                            error = phm->calculate(aw_root,cancel,0.0,trans);
1037                        }
1038
1039                        double dist_value = phm->matrix->get(0, 1);                         // distance or conformance
1040                        bool   mark       = (lowerBound <= dist_value && dist_value <= upperBound);
1041                        GB_write_flag(gb_species, mark);
1042
1043                        if (!markedSelected) {
1044                            dist_value = phm->matrix->get(0, 0);                                     // distance or conformance to self
1045                            mark       = (lowerBound <= dist_value && dist_value <= upperBound);
1046                            GB_write_flag(gb_selected, mark);
1047                           
1048                            markedSelected = true;
1049                        }
1050
1051                        delete phm;
1052                        aw_status(++count/(double)speciesCount);
1053                    }
1054
1055                    di_assert(DI_MATRIX::ROOT == 0);
1056                    DI_MATRIX::ROOT = old_root;
1057
1058                    aw_closestatus();
1059                }
1060
1061                delete ap_ratematrix;
1062                delete ap_weights;
1063                delete ap_filter;
1064                free(cancel);
1065                free(use);
1066            }
1067        }
1068
1069        free(selected);
1070        error = ta.close(error);
1071    }
1072
1073    if (error) {
1074        aw_message(error);
1075    }
1076}
1077
1078static void DI_timer(AW_root *aw_root, AW_CL cl_gbmain, AW_CL cd2) {
1079    GBDATA *gb_main = reinterpret_cast<GBDATA*>(cl_gbmain);
1080    {
1081        GB_transaction ta(gb_main);
1082        GB_tell_server_dont_wait(gb_main); // trigger database callbacks
1083    }
1084    aw_root->add_timed_callback(500, DI_timer, cl_gbmain, cd2);
1085}
1086
1087static void selected_species_changed_cb(AW_root* , AW_CL cl_viewer) {
1088    if (di_dmatrix) {
1089        AW_window *viewer = reinterpret_cast<AW_window*>(cl_viewer);
1090        void       redisplay_needed(AW_window *,DI_dmatrix *dis);
1091        redisplay_needed(viewer, di_dmatrix);
1092    }
1093}
1094
1095static void di_view_matrix_cb(AW_window *aww){
1096    GB_ERROR error = di_calculate_matrix_cb(aww,0,true);
1097    if (error) return;
1098
1099    if (!di_dmatrix) di_dmatrix = new DI_dmatrix();
1100
1101    static AW_window *viewer = 0;
1102    if (!viewer) {
1103        AW_root *aw_root = aww->get_root();
1104        viewer           = DI_create_view_matrix_window(aw_root, di_dmatrix);
1105       
1106        AW_awar *awar_sel = aw_root->awar(AWAR_SPECIES_NAME);
1107        awar_sel->add_callback(selected_species_changed_cb, AW_CL(viewer));
1108
1109        aw_root->add_timed_callback(2000, DI_timer, AW_CL(GLOBAL_gb_main), 0);
1110    }
1111
1112    di_dmatrix->init();
1113    di_dmatrix->display(false);
1114   
1115    viewer->activate();
1116}
1117
1118static void di_save_matrix_cb(AW_window *aww)
1119{
1120    // save the matrix
1121    GB_ERROR error = di_calculate_matrix_cb(aww,0, true);
1122    if (!error) {
1123        char              *filename = aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_FILENAME)->read_string();
1124        enum DI_SAVE_TYPE  type     = (enum DI_SAVE_TYPE)aww->get_root()->awar(AWAR_DIST_SAVE_MATRIX_TYPE)->read_int();
1125
1126        DI_MATRIX::ROOT->save(filename,type);
1127        free(filename);
1128    }
1129    awt_refresh_selection_box(aww->get_root(), AWAR_DIST_SAVE_MATRIX_BASE);
1130    aww->hide_or_notify(error);
1131}
1132
1133AW_window *DI_create_save_matrix_window(AW_root *aw_root, char *base_name)
1134{
1135    static AW_window_simple *aws = 0;
1136    if (aws) return aws;
1137    aws = new AW_window_simple;
1138    aws->init( aw_root, "SAVE_MATRIX", "Save Matrix");
1139    aws->load_xfig("sel_box.fig");
1140
1141    aws->at("close");aws->callback((AW_CB0)AW_POPDOWN);
1142    aws->create_button("CLOSE", "CANCEL","C");
1143
1144
1145    aws->at("help");aws->callback(AW_POPUP_HELP,(AW_CL)"save_matrix.hlp");
1146    aws->create_button("HELP", "HELP","H");
1147
1148    aws->at("user");
1149    aws->create_option_menu(AWAR_DIST_SAVE_MATRIX_TYPE);
1150    aws->insert_default_option("Phylip Format (Lower Triangular Matrix)","P",DI_SAVE_PHYLIP_COMP);
1151    aws->insert_option("Readable (using NDS)","R",DI_SAVE_READABLE);
1152    aws->insert_option("Tabbed (using NDS)","R",DI_SAVE_TABBED);
1153    aws->update_option_menu();
1154
1155    awt_create_selection_box((AW_window *)aws,base_name);
1156
1157    aws->at("save2");aws->callback(di_save_matrix_cb);
1158    aws->create_button("SAVE", "SAVE","S");
1159
1160    aws->callback( (AW_CB0)AW_POPDOWN);
1161    aws->at("cancel2");
1162    aws->create_button("CLOSE", "CANCEL","C");
1163
1164    return aws;
1165}
1166
1167static AW_window *awt_create_select_cancel_window(AW_root *aw_root)
1168{
1169    AW_window_simple *aws = new AW_window_simple;
1170    aws->init( aw_root, "SELECT_CHARS_TO_CANCEL_COLUMN", "CANCEL SELECT");
1171    aws->load_xfig("di_cancel.fig");
1172
1173    aws->at("close");aws->callback((AW_CB0)AW_POPDOWN);
1174    aws->create_button("CLOSE", "CLOSE","C");
1175
1176    aws->at("cancel");
1177    aws->create_input_field(AWAR_DIST_CANCEL_CHARS,12);
1178
1179    return (AW_window *)aws;
1180}
1181
1182static const char *enum_trans_to_string[] = {
1183    "none",
1184    "similarity",
1185    "jukes_cantor",
1186    "felsenstein",
1187
1188    "pam",
1189    "hall",
1190    "barker",
1191    "chemical",
1192
1193    "haesch",
1194    "kimura",
1195    "olsen",
1196    "felsenstein voigt",
1197    "olsen voigt",
1198    "max ml"
1199};
1200
1201static int di_calculate_tree_show_status(int loop_count, int bootstrap_count) {
1202    // returns value of aw_status()
1203    int result = 0;
1204    if (bootstrap_count) {
1205        result = aw_status(GBS_global_string("Tree #%i of %i", loop_count, bootstrap_count));
1206        if (!result) result = aw_status(loop_count/double(bootstrap_count));
1207    }
1208    else {
1209        result = aw_status(GBS_global_string("Tree #%i of %u", loop_count, UINT_MAX));
1210        if (!result) result = aw_status(loop_count/double(UINT_MAX));
1211    }
1212    return result;
1213}
1214
1215static void di_calculate_tree_cb(AW_window *aww, AW_CL bootstrap_flag)
1216{
1217    AW_root   *aw_root         = aww->get_root();
1218    GB_ERROR   error           = 0;
1219    GBT_TREE  *tree            = 0;
1220    char     **all_names       = 0;
1221    int        loop_count      = 0;   
1222    bool       close_stat      = false;
1223    int        bootstrap_count = aw_root->awar(AWAR_DIST_BOOTSTRAP_COUNT)->read_int();
1224
1225    {
1226        char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1227        error           = GBT_check_tree_name(tree_name);
1228        free(tree_name);
1229    }
1230
1231    int starttime        = time(0);
1232    int trees_per_second = -1;
1233    int last_status_upd  = 0;
1234
1235    if (!error) {
1236        {
1237            const char *stat_msg = 0;
1238            if (bootstrap_flag) {
1239                if (bootstrap_count) stat_msg = "Calculating bootstrap trees";
1240                else stat_msg = "Calculating bootstrap trees (KILL to stop)";
1241            }
1242            else stat_msg = "Calculating tree";
1243
1244            aw_openstatus(stat_msg);
1245            close_stat = true;
1246        }
1247
1248        if (bootstrap_flag){
1249            loop_count++;
1250            di_calculate_tree_show_status(loop_count, bootstrap_count);
1251
1252            gb_assert(DI_MATRIX::ROOT == 0);
1253            di_calculate_matrix_cb(aww,bootstrap_flag, true);
1254
1255            DI_MATRIX *matr = DI_MATRIX::ROOT;
1256            if (!matr) {
1257                error = "unexpected error in di_calculate_matrix_cb (data missing)";
1258            }
1259            else {
1260                all_names = (char **)calloc(sizeof(char *),(size_t)matr->nentries+2);
1261                for (long i=0;i<matr->nentries;i++){
1262                    all_names[i] = strdup(matr->entries[i]->name);
1263                }
1264                ctree_init(matr->nentries,all_names);
1265            }
1266        }
1267    }
1268
1269    do {
1270        if (error) break;
1271
1272        if (bootstrap_flag) {
1273            bool show_update = true;
1274            if (trees_per_second == -1) {
1275                int timediff = time(0)-starttime;
1276                if (timediff>10) {
1277                    trees_per_second = loop_count/timediff;
1278                }
1279            }
1280            else {
1281                if ((loop_count-last_status_upd)*4 < trees_per_second) {
1282                    show_update = false; // do not update more often than 4 times a second
1283                }
1284            }
1285
1286            loop_count++;
1287            if (show_update) { // max. once per second
1288                if (di_calculate_tree_show_status(loop_count, bootstrap_count)) {
1289                    break; // user abort
1290                }
1291                last_status_upd = loop_count;
1292            }
1293
1294        }
1295
1296        delete DI_MATRIX::ROOT;
1297        DI_MATRIX::ROOT = 0;
1298
1299        error = di_calculate_matrix_cb(aww, bootstrap_flag, bootstrap_flag ? false : true);
1300        if (error && strcmp(error,"Aborted") == 0) {
1301            error = 0;          // clear error (otherwise no tree will be read below)
1302            break;              // end of bootstrap
1303        }
1304
1305        if (!DI_MATRIX::ROOT) {
1306            error = "unexpected error in di_calculate_matrix_cb (data missing)";
1307            break;
1308        }
1309
1310        DI_MATRIX  *matr  = DI_MATRIX::ROOT;
1311        char     **names = (char **)calloc(sizeof(char *),(size_t)matr->nentries+2);
1312        for (long i=0;i<matr->nentries;i++){
1313            names[i] = matr->entries[i]->name;
1314        }
1315        tree = neighbourjoining(names,matr->matrix->m,matr->nentries,sizeof(GBT_TREE));
1316
1317#if 0
1318        // save all generated trees for debugging
1319        FILE *out = fopen(GBS_global_string("tree_%i",loop_count),"w");
1320        TREE_export_tree(gb_main, out, tree, false, true, false);
1321        printf(";\n"); // aka empty tree
1322        fclose(out);
1323#endif
1324
1325        if (bootstrap_flag){
1326            insert_ctree(tree,1);
1327            GBT_delete_tree(tree); tree = 0;
1328        }
1329        free(names);
1330
1331    } while (bootstrap_flag && loop_count != bootstrap_count);
1332
1333    if (!error) {
1334        if (bootstrap_flag) tree = get_ctree();
1335
1336        char *tree_name = aw_root->awar(AWAR_DIST_TREE_STD_NAME)->read_string();
1337        GB_begin_transaction(GLOBAL_gb_main);
1338        error = GBT_write_tree(GLOBAL_gb_main,0,tree_name,tree);
1339
1340        if (!error) {
1341            char       *filter_name = AWT_get_combined_filter_name(aw_root, "dist");
1342            int         transr      = aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1343            const char *comment     = GBS_global_string("PRG=dnadist CORR=%s FILTER=%s PKG=ARB", enum_trans_to_string[transr], filter_name);
1344
1345            error = GBT_write_tree_rem(GLOBAL_gb_main,tree_name,comment);
1346            free(filter_name);
1347        }
1348
1349        error = GB_end_transaction(GLOBAL_gb_main, error);
1350        free(tree_name);
1351    }
1352
1353    if (tree) GBT_delete_tree(tree);
1354
1355    if (close_stat) aw_closestatus();
1356    aw_status();                      // remove 'abort' flag
1357
1358    if (bootstrap_flag) {
1359        if (all_names) GBT_free_names(all_names);
1360    }
1361#if defined(DEBUG)
1362    else {
1363        gb_assert(all_names == 0);
1364    }
1365#endif // DEBUG
1366
1367    delete DI_MATRIX::ROOT; DI_MATRIX::ROOT = 0;
1368
1369    if (error) aw_message(error);
1370}
1371
1372
1373static void di_autodetect_callback(AW_window *aww)
1374{
1375    GB_push_transaction(GLOBAL_gb_main);
1376    if (DI_MATRIX::ROOT){
1377        delete DI_MATRIX::ROOT;
1378        DI_MATRIX::ROOT = 0;
1379    }
1380    AW_root *aw_root = aww->get_root();
1381    char    *use = aw_root->awar(AWAR_DIST_ALIGNMENT)->read_string();
1382    long ali_len = GBT_get_alignment_len(GLOBAL_gb_main,use);
1383    if (ali_len<=0) {
1384        GB_pop_transaction(GLOBAL_gb_main);
1385        delete use;
1386        aw_message("Please select a valid alignment");
1387        return;
1388    }
1389
1390    aw_openstatus("Analyzing data...");
1391    aw_status("Read the database");
1392
1393    DI_MATRIX *phm = new DI_MATRIX(GLOBAL_gb_main,aw_root);
1394    char    *filter_str = aw_root->awar(AWAR_DIST_FILTER_FILTER)->read_string();
1395    char    *cancel = aw_root->awar(AWAR_DIST_CANCEL_CHARS)->read_string();
1396
1397    AP_filter *ap_filter = new AP_filter;
1398    long flen = strlen(filter_str);
1399    if (flen == ali_len){
1400        ap_filter->init(filter_str,"0",ali_len);
1401    }else{
1402        if (flen){
1403            aw_message("WARNING: YOUR FILTER LEN IS NOT THE ALIGNMENT LEN\nFILTER IS TRUNCATED WITH ZEROS OR CUTTED");
1404            ap_filter->init(filter_str,"0",ali_len);
1405        }else{
1406            ap_filter->init(ali_len);
1407        }
1408    }
1409
1410    AP_weights *ap_weights = new AP_weights;
1411    ap_weights->init(ap_filter);
1412
1413    char *load_what      = aw_root->awar(AWAR_DIST_WHICH_SPECIES)->read_string();
1414    char *sort_tree_name = aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->read_string();
1415   
1416    LoadWhat all_flag = (strcmp(load_what,"all") == 0) ? DI_LOAD_ALL : DI_LOAD_MARKED;
1417
1418    delete phm->ROOT;
1419    phm->ROOT = 0;
1420    phm->load(use, ap_filter, ap_weights, all_flag, sort_tree_name, true, NULL);
1421    free(sort_tree_name);
1422    GB_pop_transaction(GLOBAL_gb_main);
1423
1424    aw_status("Search Correction");
1425
1426    DI_TRANSFORMATION trans;
1427    trans = (DI_TRANSFORMATION)aw_root->awar(AWAR_DIST_CORR_TRANS)->read_int();
1428
1429
1430    char *error = phm->analyse(aw_root);
1431    aw_closestatus();
1432
1433    delete phm;
1434    if (error) aw_message(error);
1435
1436    free(load_what);
1437
1438    free(cancel);
1439    delete ap_filter;
1440    delete ap_weights;
1441
1442    free(use);
1443    free(filter_str);
1444}
1445
1446ATTRIBUTED(__ATTR__NORETURN, static void di_exit(AW_window *aww)) {
1447    if (GLOBAL_gb_main) {
1448        aww->get_root()->unlink_awars_from_DB(GLOBAL_gb_main);
1449        GB_close(GLOBAL_gb_main);
1450    }
1451    exit(0);
1452}
1453
1454static void di_calculate_full_matrix_cb(AW_window *aww){
1455    if (DI_MATRIX::ROOT && DI_MATRIX::ROOT->matrix_type == DI_MATRIX_FULL) return;
1456    delete_matrix_cb(0);
1457    di_calculate_matrix_cb(aww,0, true);
1458    if (di_dmatrix){
1459        di_dmatrix->resized();
1460        di_dmatrix->display(false);
1461    }
1462}
1463
1464static void di_compress_tree_cb(AW_window *aww){
1465    GB_transaction dummy(GLOBAL_gb_main);
1466    char *treename = aww->get_root()->awar(AWAR_DIST_TREE_COMP_NAME)->read_string();
1467    GB_ERROR error = 0;
1468    GBT_TREE *tree = GBT_read_tree(GLOBAL_gb_main,treename,sizeof(GBT_TREE));
1469    if (!tree) error = "Please select a valid tree first";
1470
1471    if (DI_MATRIX::ROOT && DI_MATRIX::ROOT->matrix_type!=DI_MATRIX_COMPRESSED){
1472        delete_matrix_cb(0);        // delete wrong matrix !!!
1473    }
1474
1475    di_calculate_matrix_cb(aww,0, true);
1476    if (!DI_MATRIX::ROOT) error = "Cannot calculate your matrix";
1477
1478    if (!error) {
1479        error = DI_MATRIX::ROOT->compress(tree);
1480    }
1481
1482    if (tree) GBT_delete_tree(tree);
1483
1484    if (error) {
1485        aw_message(error);
1486    }
1487    else {
1488        //  aww->hide();
1489        if (di_dmatrix){
1490            di_dmatrix->resized();
1491            di_dmatrix->display(false);
1492        }
1493    }
1494    free(treename);
1495}
1496
1497/*
1498  AW_window *create_select_compress_tree(AW_root *awr){
1499  static AW_window_simple *aws = 0;
1500  if (aws) return (AW_window *)aws;
1501
1502  aws = new AW_window_simple;
1503  aws->init( awr, "SELECT_TREE_TO_COMPRESS_MATRIX", "SELECT A TREE TO COMPRESS MATRIX", 400, 200 );
1504  aws->at(10,10);
1505  aws->auto_space(0,0);
1506  awt_create_selection_list_on_trees(gb_main,(AW_window *)aws,"dist/compress/tree_name");
1507  aws->at_newline();
1508
1509  aws->callback(AW_POPDOWN);
1510  aws->create_button("CLOSE", "CLOSE","C");
1511
1512  aws->callback(di_compress_tree_cb);
1513  aws->help_text("compress_tree.hlp");
1514  aws->create_button("GO", "DO IT","D");
1515
1516  aws->callback(AW_POPUP_HELP,(AW_CL)"compress_tree.hlp");
1517  aws->create_button("HELP", "HELP","H");
1518
1519  aws->window_fit();
1520  return (AW_window *)aws;
1521  }
1522*/
1523
1524static void di_define_sort_tree_name_cb(AW_window *aww) {
1525    AW_root *aw_root = aww->get_root();
1526    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1527    aw_root->awar(AWAR_DIST_TREE_SORT_NAME)->write_string(tree_name);
1528    free(tree_name);
1529}
1530static void di_define_compression_tree_name_cb(AW_window *aww) {
1531    AW_root *aw_root = aww->get_root();
1532    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1533    aw_root->awar(AWAR_DIST_TREE_COMP_NAME)->write_string(tree_name);
1534    free(tree_name);
1535}
1536
1537static void di_define_save_tree_name_cb(AW_window *aww) {
1538    AW_root *aw_root = aww->get_root();
1539    char *tree_name = aw_root->awar(AWAR_DIST_TREE_CURR_NAME)->read_string();
1540    aw_root->awar(AWAR_DIST_TREE_STD_NAME)->write_string(tree_name);
1541    free(tree_name);
1542}
1543
1544
1545AW_window *DI_create_matrix_window(AW_root *aw_root) {
1546    AW_window_simple_menu *aws = new AW_window_simple_menu;
1547    aws->init( aw_root, "NEIGHBOUR JOINING", "NEIGHBOUR JOINING");
1548    aws->load_xfig("di_ge_ma.fig");
1549    aws->button_length( 10 );
1550
1551    aws->at("close");
1552    aws->callback(di_exit);
1553    aws->create_button("CLOSE", "CLOSE","C");
1554
1555    aws->at("help");
1556    aws->callback(AW_POPUP_HELP,(AW_CL)"dist.hlp");
1557    aws->create_button("HELP", "HELP","H");
1558
1559
1560    GB_push_transaction(GLOBAL_gb_main);
1561
1562    //  aws->at("close");aws->callback((AW_CB0)AW_POPDOWN);
1563    //  aws->create_button("CLOSE","C");
1564
1565    aws->create_menu("FILE", "F", "HELP for Window", AWM_ALL);
1566    aws->insert_menu_topic("macros", "Macros  ...", "M", "macro.hlp", AWM_ALL, (AW_CB)AW_POPUP, (AW_CL)awt_open_macro_window,(AW_CL)"NEIGHBOUR_JOINING");
1567    aws->insert_menu_topic("quit",   "Quit",        "Q", "quit.hlp",  AWM_ALL, (AW_CB)di_exit,  0,  0                                                 );
1568
1569    aws->create_menu("Properties", "P", "properties.hlp", AWM_ALL);
1570    aws->insert_menu_topic("frame_props", "Frame ...",                                 "F", "props_frame.hlp", AWM_ALL, AW_POPUP, (AW_CL)AW_preset_window, 0);
1571    aws->insert_menu_topic("save_props",  "Save Properties (in ~/.arb_prop/dist.arb)", "S", "savedef.hlp",     AWM_ALL, (AW_CB)AW_save_defaults, 0, 0       );
1572
1573    aws->insert_help_topic("help ...", "h", "dist.hlp", AWM_ALL, (AW_CB)AW_POPUP_HELP, (AW_CL)"dist.hlp", 0);
1574
1575    aws->at( "which_species" );
1576    aws->create_option_menu( AWAR_DIST_WHICH_SPECIES );
1577    aws->insert_option( "all", "a", "all" );
1578    aws->insert_default_option( "marked",  "m", "marked" );
1579    aws->update_option_menu();
1580
1581    aws->at("which_alignment");
1582    awt_create_selection_list_on_ad(GLOBAL_gb_main, (AW_window *)aws,AWAR_DIST_ALIGNMENT,"*=");
1583
1584    aws->at("filter_select");
1585    dist_filter_cl = awt_create_select_filter(aws->get_root(), GLOBAL_gb_main, AWAR_DIST_FILTER_NAME);
1586    aws->callback(AW_POPUP, (AW_CL)awt_create_select_filter_win, (AW_CL)dist_filter_cl);
1587    aws->create_button("SELECT_FILTER", AWAR_DIST_FILTER_NAME);
1588
1589    aws->at("weights_select");
1590    aws->sens_mask(AWM_EXP);
1591    global_csp = new AWT_csp(GLOBAL_gb_main,aws->get_root(),AWAR_DIST_COLUMN_STAT_NAME);
1592    aws->callback(AW_POPUP,(AW_CL)create_csp_window, (AW_CL)global_csp);
1593    aws->create_button("SELECT_COL_STAT",AWAR_DIST_COLUMN_STAT_NAME);
1594    aws->sens_mask(AWM_ALL);
1595
1596    aws->at("which_cancel");
1597    aws->create_input_field(AWAR_DIST_CANCEL_CHARS,12);
1598
1599    aws->at("cancel_select");
1600    aws->callback((AW_CB1)AW_POPUP,(AW_CL)awt_create_select_cancel_window);
1601    aws->create_button("SELECT_CANCEL_CHARS", "Info","C");
1602
1603    aws->at("change_matrix");
1604    aws->callback(AW_POPUP,(AW_CL)create_dna_matrix_window,0);
1605    aws->create_button("EDIT_MATRIX", "Edit Matrix");
1606
1607    aws->at("enable");
1608    aws->create_toggle(AWAR_DIST_MATRIX_DNA_ENABLED);
1609
1610    aws->at("which_correction");
1611    aws->create_option_menu(AWAR_DIST_CORR_TRANS, NULL ,"");
1612    aws->insert_option("none",                    "n", (int)DI_TRANSFORMATION_NONE               );
1613    aws->insert_option("similarity",              "n", (int)DI_TRANSFORMATION_SIMILARITY         );
1614    aws->insert_option("jukes-cantor (dna)",      "c", (int)DI_TRANSFORMATION_JUKES_CANTOR       );
1615    aws->insert_option("felsenstein (dna)",       "f", (int)DI_TRANSFORMATION_FELSENSTEIN        );
1616    aws->insert_option("olsen (dna)",             "o", (int)DI_TRANSFORMATION_OLSEN              );
1617    aws->insert_option("felsenstein/voigt (exp)", "1", (int)DI_TRANSFORMATION_FELSENSTEIN_VOIGT  );
1618    aws->insert_option("olsen/voigt (exp)",       "2", (int)DI_TRANSFORMATION_OLSEN_VOIGT        );
1619    aws->insert_option("haesch (exp)",            "f", (int)DI_TRANSFORMATION_HAESCH             );
1620    aws->insert_option("kimura (pro)",            "k", (int)DI_TRANSFORMATION_KIMURA             );
1621    aws->insert_option("PAM (protein)",           "c", (int)DI_TRANSFORMATION_PAM                );
1622    aws->insert_option("Cat. Hall(exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_HALL    );
1623    aws->insert_option("Cat. Barker(exp)",        "c", (int)DI_TRANSFORMATION_CATEGORIES_BARKER  );
1624    aws->insert_option("Cat.Chem (exp)",          "c", (int)DI_TRANSFORMATION_CATEGORIES_CHEMICAL);
1625    // aws->insert_option("Maximum Likelihood (exp)", "l", (int)DI_TRANSFORMATION_NONE);
1626    aws->insert_default_option("unknown", "u", (int)DI_TRANSFORMATION_NONE);
1627
1628    aws->update_option_menu();
1629
1630    // -------------------
1631    //      right side
1632    // -------------------
1633
1634    aws->at("mark_distance");
1635    aws->callback(di_mark_by_distance);
1636    aws->create_autosize_button("MARK_BY_DIST", "Mark all species");
1637
1638    aws->at("mark_lower");
1639    aws->create_input_field(AWAR_DIST_MIN_DIST, 5);
1640   
1641    aws->at("mark_upper");
1642    aws->create_input_field(AWAR_DIST_MAX_DIST, 5);
1643
1644    // -----------------
1645
1646    aws->button_length(18);
1647
1648    //     aws->at("compress");
1649    //     aws->callback(AW_POPUP,  (AW_CL)create_select_compress_tree,0);
1650    //     aws->help_text("compress_tree.hlp");
1651    //     aws->create_button("CALC_COMPRESSED_MATRIX", "Calculate\nCompressed M. ...","C");
1652
1653    aws->at("compress");
1654    aws->callback(di_compress_tree_cb);
1655    aws->create_button("CALC_COMPRESSED_MATRIX", "Calculate\nCompressed Matrix","C");
1656
1657    aws->at("calculate");
1658    aws->callback(di_calculate_full_matrix_cb);
1659    aws->create_button("CALC_FULL_MATRIX", "Calculate\nFull Matrix","F");
1660
1661    aws->button_length(13);
1662
1663    aws->at("save_matrix");
1664    aws->callback(AW_POPUP, (AW_CL)DI_create_save_matrix_window,(AW_CL)AWAR_DIST_SAVE_MATRIX_BASE);
1665    aws->create_button("SAVE_MATRIX", "Save matrix","M");
1666
1667    aws->at("view_matrix");
1668    aws->callback(di_view_matrix_cb);
1669    aws->create_button("VIEW_MATRIX", "View matrix","V");
1670
1671    aws->at("tree_list");
1672    awt_create_selection_list_on_trees(GLOBAL_gb_main,(AW_window *)aws,AWAR_DIST_TREE_CURR_NAME);
1673
1674    aws->button_length(18);
1675
1676    aws->at("t_calculate");
1677    aws->callback(di_calculate_tree_cb,0);
1678    aws->create_button("CALC_TREE", "Calculate \ntree","C");
1679
1680    aws->at("bootstrap");
1681    aws->callback(di_calculate_tree_cb,1);
1682    aws->create_button("CALC_BOOTSTRAP_TREE", "Calculate \nbootstrap tree");
1683
1684    aws->at("bcount");
1685    aws->create_input_field(AWAR_DIST_BOOTSTRAP_COUNT, 7);
1686
1687    aws->at("autodetect");   // auto
1688    aws->callback(di_autodetect_callback);
1689    aws->sens_mask(AWM_EXP);
1690    aws->create_button("AUTODETECT_CORRECTION", "AUTODETECT","A");
1691    aws->sens_mask(AWM_ALL);
1692
1693    aws->at("sort_tree_name");
1694    aws->create_input_field(AWAR_DIST_TREE_SORT_NAME,12);
1695    aws->at("compr_tree_name");
1696    aws->create_input_field(AWAR_DIST_TREE_COMP_NAME,12);
1697    aws->at("calc_tree_name");
1698    aws->create_input_field(AWAR_DIST_TREE_STD_NAME,12);
1699
1700    aws->button_length(22);
1701
1702    aws->at("use_compr_tree");
1703    aws->callback(di_define_compression_tree_name_cb);
1704    aws->create_button("USE_COMPRESSION_TREE", "Use to compress", "");
1705
1706    aws->at("use_sort_tree");
1707    aws->callback(di_define_sort_tree_name_cb);
1708    aws->create_button("USE_SORT_TREE", "Use to sort", "");
1709
1710    aws->at("use_existing");
1711    aws->callback(di_define_save_tree_name_cb);
1712    aws->create_button("USE_NAME", "Use as new tree name", "");
1713
1714    GB_pop_transaction(GLOBAL_gb_main);
1715    return (AW_window *)aws;
1716}
Note: See TracBrowser for help on using the repository browser.