source: tags/arb_5.3/PROBE/PT_buildtree.cxx

Last change on this file was 5919, checked in by westram, 15 years ago
  • PT-server
    • DEVEL_JB → ARB_64
    • save information block (containing magic and db-version)
    • fail on 32bit when trying to load big database (>4Gb) created with 64-bit version
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.8 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 enter_stage_1_build_tree(PT_main * main,char *tname) {
175    main = main;
176    POS_TREE       *pt = NULL;
177    int             i, j, abs_align_pos;
178    static  int psize;
179    char        *align_abs;
180    pt = 0;
181    char *t2name;
182
183    FILE *out;                      /* write params */
184    long    pos;
185    long    last_obj;
186
187    if (unlink(tname) ) {
188        if (GB_size_of_file(tname)>=0) {
189            fprintf(stderr,"Cannot remove %s\n",tname);
190            exit(EXIT_FAILURE);
191        }
192    }
193
194    t2name = (char *)calloc(sizeof(char),strlen(tname) + 2);
195    sprintf(t2name,"%s%%",tname);
196    out = fopen(t2name,"w");                /* open file for output */
197    if (!out) {
198        fprintf(stderr,"Cannot open %s\n",t2name);
199        exit(EXIT_FAILURE);
200    }
201    GB_set_mode_of_file(t2name,0666);
202
203    putc(0,out);        /* disable zero father */
204    pos = 1;
205
206    // now temp file exists -> trigger ptserver-selectionlist-update in all
207    // ARB applications by writing to log
208    GBS_add_ptserver_logentry(GBS_global_string("Calculating probe tree (%s)", tname));
209
210    psg.ptmain = PT_init(PT_B_MAX);
211    psg.ptmain->stage1 = 1;             /* enter stage 1 */
212
213    pt = PT_create_leaf(psg.ptmain,NULL,PT_N,0,0,0);       /* create main node */
214    pt = PT_change_leaf_to_node(psg.ptmain,pt);
215    psg.stat.cut_offs = 0;                  /* statistic information */
216    GB_begin_transaction(psg.gb_main);
217
218    char partstring[256];
219    int  partsize = 0;
220    int  pass0    = 0;
221    int  passes   = 1;
222
223    {
224        ULONG total_size = psg.char_count;
225
226        printf("Overall bases: %li\n", total_size);
227
228        while (1) {
229#ifdef ARB_64
230            ULONG estimated_kb = (total_size/1024)*55;  // value by try and error (for gene server)
231                                                        // TODO: estimated_kb depends on 32/64 bit...
232#else
233            ULONG estimated_kb = (total_size/1024)*35; /* value by try and error; 35 bytes per base*/
234#endif           
235            printf("Estimated memory usage for %i passes: %lu k\n", passes, estimated_kb);
236
237            if (estimated_kb <= physical_memory) break;
238
239            total_size /= 4;
240            partsize ++;
241            passes     *= 5;
242        }
243
244        //     while ( ULONG(total_size*35/1024) > physical_memory) {  /* value by try and error; 35 bytes per base*/
245        //         total_size /= 4;
246        //         partsize ++;
247        //         passes *=5;
248        //     }
249    }
250
251    printf ("Tree Build: %li bases in %i passes\n",psg.char_count, passes);
252
253    PT_init_base_string_counter(partstring,PT_N,partsize);
254    while (!PT_base_string_counter_eof(partstring)) {
255        ++pass0;
256        printf("\n Starting pass %i(%i)\n", pass0, passes);
257
258        for (i = 0; i < psg.data_count; i++) {
259            align_abs = probe_read_alignment(i, &psize);
260
261            if ((i % 1000) == 0) {
262                printf("%i(%i)\t cutoffs:%i\n", i, psg.data_count, psg.stat.cut_offs);
263                fflush(stdout);
264            }
265            abs_align_pos = psize-1;
266            for (j = psg.data[i].size - 1; j >= 0; j--, abs_align_pos--) {
267                get_abs_align_pos(align_abs, abs_align_pos); // may result in neg. abs_align_pos (seems to happen if sequences are short < 214bp )
268                if (abs_align_pos < 0) break; // -> in this case abort
269
270                if ( partsize && (*partstring!= psg.data[i].data[j] || strncmp(partstring,psg.data[i].data+j,partsize)) ) continue;
271                if (ptd_string_shorter_than(psg.data[i].data+j,9)) continue;
272
273                pt = build_pos_tree(pt, j, abs_align_pos, i, psg.data[i].size);
274            }
275            free(align_abs);
276        }
277        pos = PTD_save_partial_tree(out, psg.ptmain, pt, partstring, partsize, pos, &last_obj);
278#ifdef PTM_DEBUG
279        PTM_debug_mem();
280        PTD_debug_nodes();
281#endif
282        PT_inc_base_string_count(partstring,PT_N,PT_B_MAX,partsize);
283    }
284    if (partsize){
285        pos = PTD_save_partial_tree(out, psg.ptmain, pt, NULL, 0, pos, &last_obj);
286#ifdef PTM_DEBUG
287        PTM_debug_mem();
288        PTD_debug_nodes();
289#endif
290    }
291
292    bool need64bit                        = false;  // does created db need a 64bit ptserver ?
293#ifdef ARB_64
294    if (last_obj >= 0xffffffff) need64bit = true;   // last_obj is bigger than int
295#else
296    if (last_obj <= 0) {                            // overflow ?
297        GBK_terminate("Overflow - out of memory");
298    }
299#endif
300
301    // write information about database
302    long info_pos = pos;
303    PTD_put_int(out, PT_SERVER_MAGIC);              // marker to identify PT-Server file
304    fputc(PT_SERVER_VERSION, out);                  // version of PT-Server file
305    pos += 4+1;
306
307    // as last element of info block, write it's size (2byte)
308    long info_size = pos-info_pos;
309    PTD_put_short(out, info_size);
310    pos += 2;
311
312    // save DB footer (which is the entry point on load)
313
314    if (need64bit) {                                // last_obj is bigger than int
315#ifdef ARB_64
316        PTD_put_longlong(out, last_obj);            // write last_obj as long long (64 bit)
317        PTD_put_int(out, 0xffffffff);               // write 0xffffffff at the end to signal 64bit ptserver is needed
318#else
319        gb_assert(0);
320#endif
321    }
322    else {
323        PTD_put_int(out,last_obj);                  // last_obj fits into an int -> store it as usual (compatible to old unversioned format)
324    }
325
326    GB_commit_transaction(psg.gb_main);
327    fclose(out);
328    if (GB_rename_file(t2name,tname) ) {
329        GB_print_error();
330    }
331
332    if (GB_set_mode_of_file(tname,00666)) {
333        GB_print_error();
334    }
335    psg.pt = pt;
336}
337
338/* load tree from disk */
339void enter_stage_3_load_tree(PT_main * main,char *tname)
340{
341    main = main;
342    psg.ptmain = PT_init(PT_B_MAX);
343    psg.ptmain->stage3 = 1;             /* enter stage 3 */
344
345    FILE *in;
346    long size = GB_size_of_file(tname);
347    if (size<0) {
348        fprintf(stderr,"PT_SERVER: Error while opening file %s\n",tname);
349        exit(EXIT_FAILURE);
350    }
351
352    printf ("PT_SERVER: Reading Tree %s of size %li from disk\n",tname,size);
353    in = fopen(tname,"r");
354    if (!in) {
355        fprintf(stderr,"PT_SERVER: Error while opening file %s\n",tname);
356        exit(EXIT_FAILURE);
357    }
358    PTD_read_leafs_from_disk(tname,psg.ptmain,&psg.pt);
359    fclose(in);
360}
361
362#define DEBUG_MAX_CHAIN_SIZE 1000
363#define DEBUG_TREE_DEEP PT_POS_TREE_HEIGHT+50
364
365struct PT_debug_struct {
366    int chainsizes[DEBUG_MAX_CHAIN_SIZE][DEBUG_TREE_DEEP];
367    int chainsizes2[DEBUG_MAX_CHAIN_SIZE];
368    int splits[DEBUG_TREE_DEEP][PT_B_MAX];
369    int nodes[DEBUG_TREE_DEEP];
370    int tips[DEBUG_TREE_DEEP];
371    int chains[DEBUG_TREE_DEEP];
372    int chaincount;
373    } *ptds;
374
375struct PT_chain_debug {
376  int operator() (int name, int apos, int rpos)
377    {
378    psg.height++;
379    name = name;
380    apos = apos;
381    rpos = rpos;
382    return 0;
383}
384};
385void PT_analyse_tree(POS_TREE *pt,int height)
386    {
387    PT_NODE_TYPE type;
388    POS_TREE *pt_help;
389    int i;
390    int basecnt;
391    type = PT_read_type(pt);
392    switch (type) {
393        case PT_NT_NODE:
394            ptds->nodes[height]++;
395            basecnt = 0;
396            for (i=PT_QU; i<PT_B_MAX; i++) {
397                if ((pt_help = PT_read_son(psg.ptmain,pt, (PT_BASES)i)))
398                {
399                    basecnt++;
400                    PT_analyse_tree(pt_help,height+1);
401                };
402            }
403            ptds->splits[height][basecnt]++;
404            break;
405        case PT_NT_LEAF:
406            ptds->tips[height]++;
407            break;
408        default:
409            ptds->chains[height]++;
410            psg.height = 0;
411            PT_read_chain(psg.ptmain,pt, PT_chain_debug());
412            if (psg.height >= DEBUG_MAX_CHAIN_SIZE) psg.height = DEBUG_MAX_CHAIN_SIZE;
413            ptds->chainsizes[psg.height][height]++;
414            ptds->chainsizes2[psg.height]++;
415            ptds->chaincount++;
416            if (ptds->chaincount<20) {
417                printf("\n\n\n\n");
418                PT_read_chain(psg.ptmain,pt, PTD_chain_print());
419            }
420            break;
421    };
422}
423
424void PT_debug_tree(void)    // show various debug information about the tree
425    {
426    ptds = (struct PT_debug_struct *)calloc(sizeof(struct PT_debug_struct),1);
427    PT_analyse_tree(psg.pt,0);
428    int i,j,k;
429    int sum = 0;            // sum of chains
430    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
431        k = ptds->nodes[i];
432        if (k){
433            sum += k;
434            printf("nodes at deep %i: %i        sum %i\t\t",i,k,sum);
435            for (j=0;j<PT_B_MAX; j++) {
436                k =     ptds->splits[i][j];
437                printf("%i:%i\t",j,k);
438            }
439            printf("\n");
440        }
441    }
442    sum = 0;
443    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
444        k = ptds->tips[i];
445        sum += k;
446        if (k)  printf("tips at deep %i: %i     sum %i\n",i,k,sum);
447    }
448    sum = 0;
449    for (i=0;i<DEBUG_TREE_DEEP;i++ ) {
450        k = ptds->chains[i];
451        if (k)  {
452            sum += k;
453            printf("chains at deep %i: %i       sum %i\n",i,k,sum);
454        }
455    }
456    sum = 0;
457    int sume = 0;           // sum of chain entries
458    for (i=0;i<DEBUG_MAX_CHAIN_SIZE;i++ ) {
459            k =     ptds->chainsizes2[i];
460            sum += k;
461            sume += i*k;
462            if (k) printf("chain of size %i occur %i    sc %i   sce %i\n",i,k,sum,sume);
463    }
464#if 0
465    sum = 0;
466    for (i=0;i<DEBUG_MAX_CHAIN_SIZE;i++ ) {
467        for (j=0;j<DEBUG_TREE_DEEP;j++ ) {
468            k =     ptds->chainsizes[i][j];
469            sum += k;
470            if (k) printf("chain of size %i at deep %i occur %i     sum:%i\n",i,j,k,sum);
471        }
472    }
473#endif
474}
Note: See TracBrowser for help on using the repository browser.