source: trunk/PROBE/PT_buildtree.cxx @ 5828

Last change on this file since 5828 was 5828, checked in by westram, 16 years ago
  • changed ulong→ULONG and uint→UINT (ulong/uint are missing under OSX)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.3 KB
Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4// #include <malloc.h>
5#include <unistd.h>
6#include <stdint.h>
7
8#include <PT_server.h>
9#include <PT_server_prototypes.h>
10#include "probe.h"
11#include "probe_tree.hxx"
12#include "pt_prototypes.h"
13
14POS_TREE *build_pos_tree (POS_TREE *pt, int anfangs_pos, int apos, int RNS_nr, unsigned int end)
15{
16    static POS_TREE       *pthelp, *pt_next;
17    unsigned int i,j;
18    int          height = 0, anfangs_apos_ref, anfangs_rpos_ref, RNS_nr_ref;
19    pthelp = pt;
20    i = anfangs_pos;
21    while (PT_read_type(pthelp) == PT_NT_NODE) {    /* now we got an inner node */
22        if ((pt_next = PT_read_son_stage_1(psg.ptmain,pthelp, (PT_BASES)psg.data[RNS_nr].data[i])) == NULL) {
23                            // there is no son of that type -> simply add the new son to that path //
24            if (pthelp == pt) { // now we create a new root structure (size will change)
25                PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr].data[i], anfangs_pos, apos, RNS_nr);
26                return pthelp;  // return the new root
27            } else {
28                PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr].data[i], anfangs_pos, apos, RNS_nr);
29                return pt;  // return the old root
30            }
31        } else {            // go down the tree
32            pthelp = pt_next;
33            height++;
34            i++;
35            if (i >= end) {     // end of sequence reached -> change node to chain and add
36                        // should never be reached, because of the terminal symbol
37                        // of each sequence
38                if (PT_read_type(pthelp) == PT_NT_CHAIN)
39                    pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr,apos,anfangs_pos);
40                // if type == node then forget it
41                return pt;
42            }
43        }
44    }
45            // type == leaf or chain
46    if (PT_read_type(pthelp) == PT_NT_CHAIN) {      // old chain reached
47        pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr, apos, anfangs_pos);
48        return pt;
49    }
50    anfangs_rpos_ref = PT_read_rpos(psg.ptmain,pthelp); // change leave to node and create two sons
51    anfangs_apos_ref = PT_read_apos(psg.ptmain,pthelp);
52    RNS_nr_ref = PT_read_name(psg.ptmain,pthelp);
53    j = anfangs_rpos_ref + height;
54
55    while (psg.data[RNS_nr].data[i] == psg.data[RNS_nr_ref].data[j]) {  /* creates nodes until sequences are different */
56                                        // type != nt_node
57        if (PT_read_type(pthelp) == PT_NT_CHAIN) {          /* break */
58            pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr, apos, anfangs_pos);
59            return pt;
60        }
61        if (height >= PT_POS_TREE_HEIGHT) {
62            if (PT_read_type(pthelp) == PT_NT_LEAF) {
63                pthelp = PT_leaf_to_chain(psg.ptmain,pthelp);
64            }
65            pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr, apos, anfangs_pos);
66            return pt;
67        }
68        if ( ((i + 1)>= end) && (j + 1 >= (unsigned)(psg.data[RNS_nr_ref].size))) { /* end of both sequences */
69            return pt;
70        }
71        pthelp = PT_change_leaf_to_node(psg.ptmain,pthelp); /* change tip to node and append two new leafs */
72        if (i + 1 >= end) {                 /* end of source sequence reached   */
73            pthelp = PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr_ref].data[j],
74                    anfangs_rpos_ref, anfangs_apos_ref, RNS_nr_ref);
75            return pt;
76        }
77        if (j + 1 >= (unsigned)(psg.data[RNS_nr_ref].size)) {       /* end of reference sequence */
78            pthelp = PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr].data[i], anfangs_pos, apos, RNS_nr);
79            return pt;
80        }
81        pthelp = PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr].data[i], anfangs_rpos_ref, anfangs_apos_ref, RNS_nr_ref);
82                    // dummy leaf just to create a new node; may become a chain
83        i++;
84        j++;
85        height++;
86    }
87    if (height >= PT_POS_TREE_HEIGHT) {
88        if (PT_read_type(pthelp) == PT_NT_LEAF)
89            pthelp = PT_leaf_to_chain(psg.ptmain,pthelp);
90        pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr, apos, anfangs_pos);
91        return pt;
92    }
93    if (PT_read_type(pthelp) == PT_NT_CHAIN) {
94        pthelp = PT_add_to_chain(psg.ptmain,pthelp, RNS_nr, apos, anfangs_pos);
95    } else {
96        pthelp = PT_change_leaf_to_node(psg.ptmain,pthelp); /* Blatt loeschen */
97        PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr].data[i], anfangs_pos, apos, RNS_nr);  /* zwei neue Blaetter */
98        PT_create_leaf(psg.ptmain,&pthelp, (PT_BASES)psg.data[RNS_nr_ref].data[j], anfangs_rpos_ref, anfangs_apos_ref, RNS_nr_ref);
99    }
100    return pt;
101}
102
103/* get the absolute alignment position */
104static GB_INLINE void get_abs_align_pos(char *seq, int &pos)
105{
106    int q_exists = 0;
107    if(pos > 3) {
108        pos-=3;
109        while(pos > 0) {
110            uint32_t c= *((uint32_t*)(seq+pos));
111            if(c == 0x2E2E2E2E) {
112                q_exists = 1;
113                pos-=4;
114                continue;
115            }
116            if(c == 0x2D2D2D2D) {
117                pos-=4;
118                continue;
119            }
120            break;
121        }
122        pos+=3;
123    }
124    while(pos) {
125        unsigned char c= seq[pos];
126        if(c == '.') {
127            q_exists = 1;
128            pos--;
129            continue;
130        }
131        if(c == '-') {
132            pos--;
133            continue;
134        }
135        break;
136    }
137    pos+=q_exists;
138}
139
140long PTD_save_partial_tree(FILE *out,PTM2 *ptmain,POS_TREE * node,char *partstring,int partsize,long pos, long *ppos)
141{
142    if (partsize) {
143        POS_TREE *son = PT_read_son(ptmain,node,(enum PT_bases_enum)partstring[0]);
144        if (son) {
145            pos = PTD_save_partial_tree(out,ptmain,son,partstring+1,partsize-1,pos,ppos);
146        }
147    }else{
148        PTD_clear_fathers(ptmain,node);
149        long r_pos;
150        int blocked;
151        blocked = 1;
152        while (blocked) {
153            blocked = 0;
154            printf("*************** pass %li *************\n",pos);
155            fflush(stdout);
156            r_pos = PTD_write_leafs_to_disk(out,ptmain,node,pos,ppos,&blocked);
157            if (r_pos > pos) pos = r_pos;
158        }
159    }
160    return pos;
161}
162
163inline int ptd_string_shorter_than(const char *s, int len) {
164    int i;
165    for (i=0;i<len;i++){
166        if (!s[i]) {
167            return 1;
168        }
169    }
170    return 0;
171}
172
173/*initialize tree and call the build pos tree procedure*/
174void
175enter_stage_1_build_tree(PT_main * main,char *tname)
176{
177    main = main;
178    POS_TREE       *pt = NULL;
179    int             i, j, abs_align_pos;
180    static  int psize;
181    char        *align_abs;
182    pt = 0;
183    char *t2name;
184
185    FILE *out;                      /* write params */
186    long    pos;
187    long    last_obj;
188
189    if (unlink(tname) ) {
190        if (GB_size_of_file(tname)>=0) {
191            fprintf(stderr,"Cannot remove %s\n",tname);
192            exit(EXIT_FAILURE);
193        }
194    }
195
196    t2name = (char *)calloc(sizeof(char),strlen(tname) + 2);
197    sprintf(t2name,"%s%%",tname);
198    out = fopen(t2name,"w");                /* open file for output */
199    if (!out) {
200        fprintf(stderr,"Cannot open %s\n",t2name);
201        exit(EXIT_FAILURE);
202    }
203    GB_set_mode_of_file(t2name,0666);
204    putc(0,out);        /* disable zero father */
205    pos = 1;
206
207    // now temp file exists -> trigger ptserver-selectionlist-update in all
208    // ARB applications by writing to log
209    GBS_add_ptserver_logentry(GBS_global_string("Calculating probe tree (%s)", tname));
210
211    psg.ptmain = PT_init(PT_B_MAX);
212    psg.ptmain->stage1 = 1;             /* enter stage 1 */
213
214    pt = PT_create_leaf(psg.ptmain,NULL,PT_N,0,0,0);       /* create main node */
215    pt = PT_change_leaf_to_node(psg.ptmain,pt);
216    psg.stat.cut_offs = 0;                  /* statistic information */
217    GB_begin_transaction(psg.gb_main);
218
219    char partstring[256];
220    int  partsize = 0;
221    int  pass0    = 0;
222    int  passes   = 1;
223
224    {
225        ULONG total_size = psg.char_count;
226
227        printf("Overall bases: %li\n", total_size);
228
229        while (1) {
230#ifdef DEVEL_JB
231            ULONG estimated_kb = (total_size/1024)*55;  // value by try and error (for gene server)
232                                                        // TODO: estimated_kb depends on 32/64 bit...
233#else
234            ULONG estimated_kb = (total_size/1024)*35; /* value by try and error; 35 bytes per base*/
235#endif           
236            printf("Estimated memory usage for %i passes: %lu k\n", passes, estimated_kb);
237
238            if (estimated_kb <= physical_memory) break;
239
240            total_size /= 4;
241            partsize ++;
242            passes     *= 5;
243        }
244
245        //     while ( ULONG(total_size*35/1024) > physical_memory) {  /* value by try and error; 35 bytes per base*/
246        //         total_size /= 4;
247        //         partsize ++;
248        //         passes *=5;
249        //     }
250    }
251
252    printf ("Tree Build: %li bases in %i passes\n",psg.char_count, passes);
253
254    PT_init_base_string_counter(partstring,PT_N,partsize);
255    while (!PT_base_string_counter_eof(partstring)) {
256        ++pass0;
257        printf("\n Starting pass %i(%i)\n", pass0, passes);
258
259        for (i = 0; i < psg.data_count; i++) {
260            align_abs = probe_read_alignment(i, &psize);
261
262            if ((i % 1000) == 0) {
263                printf("%i(%i)\t cutoffs:%i\n", i, psg.data_count, psg.stat.cut_offs);
264                fflush(stdout);
265            }
266            abs_align_pos = psize-1;
267            for (j = psg.data[i].size - 1; j >= 0; j--, abs_align_pos--) {
268                get_abs_align_pos(align_abs, abs_align_pos); // may result in neg. abs_align_pos (seems to happen if sequences are short < 214bp )
269                if (abs_align_pos < 0) break; // -> in this case abort
270
271                if ( partsize && (*partstring!= psg.data[i].data[j] || strncmp(partstring,psg.data[i].data+j,partsize)) ) continue;
272                if (ptd_string_shorter_than(psg.data[i].data+j,9)) continue;
273
274                pt = build_pos_tree(pt, j, abs_align_pos, i, psg.data[i].size);
275            }
276            free(align_abs);
277        }
278        pos = PTD_save_partial_tree(out, psg.ptmain, pt, partstring, partsize, pos, &last_obj);
279#ifdef PTM_DEBUG
280        PTM_debug_mem();
281        PTD_debug_nodes();
282#endif
283        PT_inc_base_string_count(partstring,PT_N,PT_B_MAX,partsize);
284    }
285    if (partsize){
286        pos = PTD_save_partial_tree(out, psg.ptmain, pt, NULL, 0, pos, &last_obj);
287#ifdef PTM_DEBUG
288        PTM_debug_mem();
289        PTD_debug_nodes();
290#endif
291    }
292#ifdef DEVEL_JB
293    if (last_obj >= 0xffffffff) {           // last_obj is bigger than int
294        PTD_put_longlong(out, last_obj);    // write last_obj as long long (64 bit)
295        PTD_put_int(out, 0xffffffff);       // write 0xffffffff at the end to signalize that
296                                            // last_obj ist stored in the 8 bytes before
297                                            // 0xffffffff is used as a signal to be compatible
298                                            // with older pt-servers
299    } else {                               
300        PTD_put_int(out,last_obj);          // lost_obj fits into an int -> store it as usual
301    }
302#else
303    PTD_put_int(out,last_obj);
304#endif   
305
306    GB_commit_transaction(psg.gb_main);
307    fclose(out);
308    if (GB_rename_file(t2name,tname) ) {
309        GB_print_error();
310    }
311
312    if (GB_set_mode_of_file(tname,00666)) {
313        GB_print_error();
314    }
315    psg.pt = pt;
316}
317
318/* load tree from disk */
319void enter_stage_3_load_tree(PT_main * main,char *tname)
320{
321    main = main;
322    psg.ptmain = PT_init(PT_B_MAX);
323    psg.ptmain->stage3 = 1;             /* enter stage 3 */
324
325    FILE *in;
326    long size = GB_size_of_file(tname);
327    if (size<0) {
328        fprintf(stderr,"PT_SERVER: Error while opening file %s\n",tname);
329        exit(EXIT_FAILURE);
330    }
331
332    printf ("PT_SERVER: Reading Tree %s of size %li from disk\n",tname,size);
333    in = fopen(tname,"r");
334    if (!in) {
335        fprintf(stderr,"PT_SERVER: Error while opening file %s\n",tname);
336        exit(EXIT_FAILURE);
337    }
338    PTD_read_leafs_from_disk(tname,psg.ptmain,&psg.pt);
339    fclose(in);
340}
341
342#define DEBUG_MAX_CHAIN_SIZE 1000
343#define DEBUG_TREE_DEEP PT_POS_TREE_HEIGHT+50
344
345struct PT_debug_struct {
346    int chainsizes[DEBUG_MAX_CHAIN_SIZE][DEBUG_TREE_DEEP];
347    int chainsizes2[DEBUG_MAX_CHAIN_SIZE];
348    int splits[DEBUG_TREE_DEEP][PT_B_MAX];
349    int nodes[DEBUG_TREE_DEEP];
350    int tips[DEBUG_TREE_DEEP];
351    int chains[DEBUG_TREE_DEEP];
352    int chaincount;
353    } *ptds;
354
355struct PT_chain_debug {
356  int operator() (int name, int apos, int rpos)
357    {
358    psg.height++;
359    name = name;
360    apos = apos;
361    rpos = rpos;
362    return 0;
363}
364};
365void PT_analyse_tree(POS_TREE *pt,int height)
366    {
367    PT_NODE_TYPE type;
368    POS_TREE *pt_help;
369    int i;
370    int basecnt;
371    type = PT_read_type(pt);
372    switch (type) {
373        case PT_NT_NODE:
374            ptds->nodes[height]++;
375            basecnt = 0;
376            for (i=PT_QU; i<PT_B_MAX; i++) {
377                if ((pt_help = PT_read_son(psg.ptmain,pt, (PT_BASES)i)))
378                {
379                    basecnt++;
380                    PT_analyse_tree(pt_help,height+1);
381                };
382            }
383            ptds->splits[height][basecnt]++;
384            break;
385        case PT_NT_LEAF:
386            ptds->tips[height]++;
387            break;
388        default:
389            ptds->chains[height]++;
390            psg.height = 0;
391            PT_read_chain(psg.ptmain,pt, PT_chain_debug());
392            if (psg.height >= DEBUG_MAX_CHAIN_SIZE) psg.height = DEBUG_MAX_CHAIN_SIZE;
393            ptds->chainsizes[psg.height][height]++;
394            ptds->chainsizes2[psg.height]++;
395            ptds->chaincount++;
396            if (ptds->chaincount<20) {
397                printf("\n\n\n\n");
398                PT_read_chain(psg.ptmain,pt, PTD_chain_print());
399            }
400            break;
401    };
402}
403
404void PT_debug_tree(void)    // show various debug information about the tree
405    {
406    ptds = (struct PT_debug_struct *)calloc(sizeof(struct PT_debug_struct),1);
407    PT_analyse_tree(psg.pt,0);
408    int i,j,k;
409    int sum = 0;            // sum of chains
410    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
411        k = ptds->nodes[i];
412        if (k){
413            sum += k;
414            printf("nodes at deep %i: %i        sum %i\t\t",i,k,sum);
415            for (j=0;j<PT_B_MAX; j++) {
416                k =     ptds->splits[i][j];
417                printf("%i:%i\t",j,k);
418            }
419            printf("\n");
420        }
421    }
422    sum = 0;
423    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
424        k = ptds->tips[i];
425        sum += k;
426        if (k)  printf("tips at deep %i: %i     sum %i\n",i,k,sum);
427    }
428    sum = 0;
429    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
430        k = ptds->chains[i];
431        if (k)  {
432            sum += k;
433            printf("chains at deep %i: %i       sum %i\n",i,k,sum);
434        }
435    }
436    sum = 0;
437    int sume = 0;           // sum of chain entries
438    for (i=0;i<DEBUG_MAX_CHAIN_SIZE;i++ ) {
439            k =     ptds->chainsizes2[i];
440            sum += k;
441            sume += i*k;
442            if (k) printf("chain of size %i occur %i    sc %i   sce %i\n",i,k,sum,sume);
443    }
444#if 0
445    sum = 0;
446    for (i=0;i<DEBUG_MAX_CHAIN_SIZE;i++ ) {
447        for (j=0;j<DEBUG_TREE_DEEP;j++ ) {
448            k =     ptds->chainsizes[i][j];
449            sum += k;
450            if (k) printf("chain of size %i at deep %i occur %i     sum:%i\n",i,j,k,sum);
451        }
452    }
453#endif
454}
Note: See TracBrowser for help on using the repository browser.