source: trunk/NALIGNER/ali_prealigner.cxx @ 5440

Last change on this file since 5440 was 5440, checked in by westram, 17 years ago
  • fix operator precedence warnings
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 44.3 KB
Line 
1#include <stdlib.h>
2#include "ali_prealigner.hxx"
3#include "ali_aligner.hxx"
4
5unsigned long   random_stat[6] = {0, 0, 0, 0, 0, 0};
6
7/*****************************************************************************
8 *
9 * structure ali_prealigner_mask
10 *
11 *****************************************************************************/
12
13/*
14 * Insert a new map
15 */
16void            ali_prealigner_mask::
17insert(ALI_MAP * in_map, float costs)
18{
19    unsigned long   i;
20
21    calls++;
22    if (map == 0) {
23        map = in_map;
24        cost_of_binding = costs;
25    } else {
26        if (costs > cost_of_binding) {
27            delete          in_map;
28            return;
29        }
30        if (map->first_base() != in_map->first_base() ||
31            map->last_base() != in_map->last_base() ||
32            map->first_reference_base() != in_map->first_reference_base() ||
33            map->last_reference_base() != in_map->last_reference_base()) {
34            ali_fatal_error("Incompatible maps",
35                            "ali_prealigner_mask::insert()");
36        }
37        if (costs < cost_of_binding) {
38            delete          map;
39            map = in_map;
40            cost_of_binding = costs;
41            last_new = calls;
42            last_joins = 0;
43            return;
44        }
45        joins++;
46        last_joins++;
47        for (i = map->first_base(); i <= map->last_base(); i++)
48            if ((!map->is_undefined(i)) && map->position(i) != in_map->position(i))
49                map->undefine(i);
50        delete          in_map;
51    }
52}
53
54/*
55 * Delete expensive parts of solution
56 */
57void            ali_prealigner_mask::
58delete_expensive(ALI_PREALIGNER_CONTEXT * context,
59                 ALI_PROFILE * profile)
60{
61    ALI_MAP        *inverse_map;
62    unsigned long   start_hel, end_hel;
63    unsigned long   start_seq, end_seq;
64    unsigned long   start_mapped, end_mapped;
65    unsigned long   start_ok=0, end_ok=0;
66    int             start_ok_flag;
67    unsigned long   found_helix;
68    unsigned long   error_counter;
69
70    unsigned long   map_pos, i, j;
71    float           max_cost, helix_cost;
72    unsigned long   helix_counter;
73    long            compl_pos;
74    unsigned char   base1, base2;
75
76    printf("MASK : calls = %ld  joins = %ld last_new = %ld last_joins = %ld\n",
77           calls, joins, last_new, last_joins);
78
79    max_cost = profile->w_sub_maximum() * context->max_cost_of_sub_percent;
80
81    /*
82     * Delete expensive Bases
83     */
84    error_counter = 0;
85    for (i = map->first_base(); i <= map->last_base(); i++) {
86        if (!(map->is_inserted(i)) &&
87            profile->w_sub(map->position(i), i) > max_cost) {
88            error_counter++;
89            if (error_counter > context->error_count)
90                map->undefine(i);
91            else {
92                if (error_counter == context->error_count) {
93                    for (j = i - error_counter + 1; j <= i; j++)
94                        map->undefine(j);
95                }
96            }
97        } else {
98            /*
99             * If error was in helix => delete helix total
100             */
101            if (error_counter > 0 &&
102                profile->is_in_helix(map->position(i - 1), &start_hel, &end_hel)) {
103                for (j = i - 1; map->position(j) >= long(start_hel); j--);
104                for (; map->position(j) <= long(end_hel); j++)
105                    map->undefine(j);
106            }
107            error_counter = 0;
108        }
109    }
110
111    /*
112     * Delete expensive Helizes
113     */
114    inverse_map = map->inverse_without_inserts();
115    for (i = inverse_map->first_base(); i <= inverse_map->last_base(); i++) {
116        /*
117         * found a helix
118         */
119        if (profile->is_in_helix(i, &start_hel, &end_hel)) {
120            if (i != start_hel)
121                ali_fatal_error("Inconsistent positions",
122                                "ali_prealigner_mask::delete_expensive()");
123            compl_pos = profile->complement_position(start_hel);
124            /*
125             * only forward bindings
126             */
127            if (compl_pos > long(end_hel)) {
128                helix_cost = 0.0;
129                helix_counter = 0;
130                while (i <= end_hel) {
131                    /*
132                     * is binding ?
133                     */
134                    if (compl_pos > 0) {
135                        if (!inverse_map->is_undefined(i)) {
136                            base1 = (profile->sequence())->base(
137                                                                inverse_map->first_reference_base() +
138                                                                inverse_map->position(i));
139                        } else {
140                            base1 = ALI_GAP_CODE;
141                        }
142                        if (!inverse_map->is_undefined(compl_pos)) {
143                            base2 = (profile->sequence())->base(
144                                                                inverse_map->first_reference_base() +
145                                                                inverse_map->position(compl_pos));
146                        } else {
147                            base2 = ALI_GAP_CODE;
148                        }
149                        if (base1 != ALI_GAP_CODE || base2 != ALI_GAP_CODE) {
150                            helix_cost += profile->w_bind(i, base1, compl_pos, base2);
151                            helix_counter++;
152                        }
153                    }
154                    i++;
155                    compl_pos = profile->complement_position(i);
156                }
157                if (helix_counter > 0) helix_cost /= helix_counter;
158                if (helix_cost > context->max_cost_of_helix) {
159                    for (j = start_hel; j <= end_hel; j++) {
160                        if (!inverse_map->is_undefined(j)) {
161                            map->undefine(inverse_map->first_reference_base() + inverse_map->position(j));
162                        }
163                    }
164                    for (j = profile->complement_position(end_hel); long(j) <= profile->complement_position(start_hel); j++) {
165                        if (!inverse_map->is_undefined(j)) {
166                            map->undefine(inverse_map->first_reference_base() + inverse_map->position(j));
167                        }
168                    }
169                }
170            }
171            i = end_hel;
172        }
173    }
174    delete          inverse_map;
175
176    /*
177     * Check for good parts
178     */
179    for (map_pos = map->first_base(); map_pos <= map->last_base(); map_pos++) {
180        /*
181         * search next defined segment
182         */
183        if (!map->is_undefined(map_pos)) {
184            /*
185             * find start and end of segment
186             */
187            start_seq = map_pos;
188            start_mapped = map->position(map_pos);
189            for (map_pos++; map_pos <= map->last_base() &&
190                     (!map->is_undefined(map_pos)); map_pos++);
191            end_seq = map_pos - 1;
192            end_mapped = map->position(end_seq);
193
194            /*
195             * Check segment for helizes
196             */
197            found_helix = 0;
198            start_ok_flag = 0;
199            for (i = start_seq; i <= end_seq; i++) {
200                if (profile->is_in_helix(map->position(i), &start_hel, &end_hel)) {
201                    found_helix++;
202                    /*
203                     * Helix is inside the segment
204                     */
205                    if (start_hel >= start_mapped && end_hel <= end_mapped) {
206                        if (start_ok_flag == 0) {
207                            start_ok = start_hel;
208                            start_ok_flag = 1;
209                        }
210                        end_ok = end_hel;
211                    }
212                }
213            }
214
215            /*
216             * Found good helizes
217             */
218            if (start_ok_flag == 1) {
219                for (i = start_seq; map->position(i) < long(start_ok); i++) map->undefine(i);
220                for (i = end_seq; map->position(i) > long(end_ok); i--) map->undefine(i);
221            } else {
222                /*
223                 * Found bad helizes
224                 */
225                if (found_helix > 0) {
226                    for (i = start_seq; i <= end_seq; i++)
227                        map->undefine(i);
228                }
229                /*
230                 * Segment without helix
231                 */
232                else {
233                    if (end_seq - start_seq + 1 >= (unsigned long)((2 * context->intervall_border) + context->intervall_center)) {
234                        for (i = start_seq; i < start_seq + context->intervall_border; i++) map->undefine(i);
235                        for (i = end_seq; i > end_seq - context->intervall_border; i--) map->undefine(i);
236                    } else {
237                        for (i = start_seq; i <= end_seq; i++) map->undefine(i);
238                    }
239                }
240            }
241        }
242    }
243}
244
245
246/*****************************************************************************
247 *
248 * class ALI_PREALIGNER  (privat)
249 *
250 *****************************************************************************/
251
252GB_INLINE float    ALI_PREALIGNER::
253minimum2(float a, float b)
254{
255    return ((a < b) ? a : b);
256}
257
258GB_INLINE float    ALI_PREALIGNER::
259minimum3(float a, float b, float c)
260{
261    return ((a < b) ? ((a < c) ? a : c) : ((b < c) ? b : c));
262}
263
264
265GB_INLINE void     ALI_PREALIGNER::
266calculate_first_column_first_cell(
267                                  ali_prealigner_cell * akt_cell)
268{
269    float           v1, v2, v3;
270
271    /***
272        v1 = profile->w_ins(start_x,start_y) + profile->w_del(start_y,start_y);
273    ***/
274    v1 = profile->w_ins_multi_cheap(start_x, start_y) +
275        profile->w_sub_multi_gap_cheap(start_y, start_y);
276    v2 = profile->w_sub(start_y, start_x);
277    v3 = v1;
278
279    akt_cell->d = minimum2(v1, v2);
280
281    if (akt_cell->d == v1)
282        path_map->set(0, 0, ALI_UP | ALI_LEFT);
283    if (akt_cell->d == v2)
284        path_map->set(0, 0, ALI_DIAG);
285}
286
287GB_INLINE void     ALI_PREALIGNER::
288calculate_first_column_cell(
289                            ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell,
290                            unsigned long pos_y)
291{
292    float           v1, v2, v3;
293    unsigned long   positiony;
294
295    positiony = start_y + pos_y;
296
297    /***
298        v1 = up_cell->d + profile->w_del(positiony,positiony);
299        v2 = profile->w_del_multi_unweighted(start_y, positiony - 1) +
300        profile->w_sub(positiony, start_x);
301        v3 = profile->w_del_multi_unweighted(start_y, positiony) +
302        profile->w_ins(start_x,positiony);
303    ***/
304    v1 = up_cell->d + profile->w_sub_gap(positiony);
305    v2 = profile->w_sub_multi_gap_cheap(start_y, positiony - 1) +
306        profile->w_sub(positiony, start_x);
307    v3 = profile->w_sub_multi_gap_cheap(start_y, positiony) +
308        profile->w_ins(start_x, positiony);
309
310    akt_cell->d = minimum3(v1, v2, v3);
311
312    if (v1 == akt_cell->d)
313        path_map->set(0, pos_y, ALI_UP);
314    if (v2 == akt_cell->d)
315        path_map->set(0, pos_y, ALI_DIAG);
316    if (v3 == akt_cell->d)
317        path_map->set(0, pos_y, ALI_LEFT);
318}
319
320void            ALI_PREALIGNER::
321calculate_first_column(ali_prealigner_column * akt_column)
322{
323    unsigned long   pos_y;
324
325    calculate_first_column_first_cell(&(*akt_column->cells)[0]);
326
327    for (pos_y = 1; pos_y < akt_column->column_length; pos_y++)
328        calculate_first_column_cell(&(*akt_column->cells)[pos_y - 1],
329                                    &(*akt_column->cells)[pos_y],
330                                    pos_y);
331}
332
333
334GB_INLINE void     ALI_PREALIGNER::
335calculate_first_cell(
336                     ali_prealigner_cell * left_cell, ali_prealigner_cell * akt_cell,
337                     unsigned long pos_x)
338{
339    float           v1, v2, v3;
340    unsigned long   positionx;
341
342    positionx = start_x + pos_x;
343
344    /***
345        v1 = profile->w_ins_multi_unweighted(start_x, positionx) +
346        profile->w_del(start_y,start_y);
347        v2 = profile->w_ins_multi_unweighted(start_x, positionx - 1) +
348        profile->w_sub(start_y, positionx);
349    ***/
350    v1 = profile->w_ins_multi_cheap(start_x, positionx) +
351        profile->w_sub_gap(start_y);
352    v2 = profile->w_ins_multi_cheap(start_x, positionx - 1) +
353        profile->w_sub(start_y, positionx);
354    v3 = left_cell->d + profile->w_ins(positionx, start_y);
355
356    akt_cell->d = minimum3(v1, v2, v3);
357
358    if (v1 == akt_cell->d)
359        path_map->set(pos_x, 0, ALI_UP);
360    if (v2 == akt_cell->d)
361        path_map->set(pos_x, 0, ALI_DIAG);
362    if (v3 == akt_cell->d)
363        path_map->set(pos_x, 0, ALI_LEFT);
364}
365
366GB_INLINE void     ALI_PREALIGNER::
367calculate_cell(
368               ali_prealigner_cell * diag_cell, ali_prealigner_cell * left_cell,
369               ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell,
370               unsigned long pos_x, unsigned long pos_y)
371{
372    float           v1, v2, v3;
373    unsigned long   positionx, positiony;
374
375    positionx = start_x + pos_x;
376    positiony = start_y + pos_y;
377
378    /***
379        v1 = up_cell->d + profile->w_del(positiony,positiony);
380    ***/
381    v1 = up_cell->d + profile->w_sub_gap(positiony);
382    v2 = diag_cell->d + profile->w_sub(positiony, positionx);
383    v3 = left_cell->d + profile->w_ins(positionx, positiony);
384
385    akt_cell->d = minimum3(v1, v2, v3);
386
387    if (v1 == akt_cell->d)
388        path_map->set(pos_x, pos_y, ALI_UP);
389    if (v2 == akt_cell->d)
390        path_map->set(pos_x, pos_y, ALI_DIAG);
391    if (v3 == akt_cell->d)
392        path_map->set(pos_x, pos_y, ALI_LEFT);
393}
394
395void            ALI_PREALIGNER::
396calculate_column(
397                 ali_prealigner_column * prev_col, ali_prealigner_column * akt_col,
398                 unsigned long pos_x)
399{
400    unsigned long   pos_y;
401
402    calculate_first_cell(&(*prev_col->cells)[0], &(*akt_col->cells)[0], pos_x);
403
404    for (pos_y = 1; pos_y < akt_col->column_length; pos_y++)
405        calculate_cell(&(*prev_col->cells)[pos_y - 1], &(*prev_col->cells)[pos_y],
406                       &(*akt_col->cells)[pos_y - 1], &(*akt_col->cells)[pos_y],
407                       pos_x, pos_y);
408}
409
410GB_INLINE void     ALI_PREALIGNER::
411calculate_last_column_first_cell(
412                                 ali_prealigner_cell * left_cell, ali_prealigner_cell * akt_cell,
413                                 unsigned long pos_x)
414{
415    float           v1, v2, v3;
416    unsigned long   positionx;
417
418    positionx = start_x + pos_x;
419
420    /***
421        v1 = profile->w_ins_multi_unweighted(start_x, positionx) +
422        profile->w_del(start_y,start_y);
423        v2 = profile->w_ins_multi_unweighted(start_x, positionx - 1) +
424        profile->w_sub(start_y, positionx);
425    ***/
426    v1 = profile->w_ins_multi_cheap(start_x, positionx) +
427        profile->w_sub_gap_cheap(start_y);
428    v2 = profile->w_ins_multi_cheap(start_x, positionx - 1) +
429        profile->w_sub(start_y, positionx);
430    v3 = left_cell->d + profile->w_ins(positionx, start_y);
431
432    akt_cell->d = minimum3(v1, v2, v3);
433
434    if (v1 == akt_cell->d)
435        path_map->set(pos_x, 0, ALI_UP);
436    if (v2 == akt_cell->d)
437        path_map->set(pos_x, 0, ALI_DIAG);
438    if (v3 == akt_cell->d)
439        path_map->set(pos_x, 0, ALI_LEFT);
440}
441
442GB_INLINE void     ALI_PREALIGNER::
443calculate_last_column_cell(
444                           ali_prealigner_cell * diag_cell, ali_prealigner_cell * left_cell,
445                           ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell,
446                           unsigned long pos_x, unsigned long pos_y)
447{
448    float           v1, v2, v3;
449    unsigned long   positionx, positiony;
450
451    positionx = start_x + pos_x;
452    positiony = start_y + pos_y;
453
454    /***
455        v1 = up_cell->d + profile->w_del(positiony,positiony);
456    ***/
457    v1 = up_cell->d + profile->w_sub_gap_cheap(positiony);
458    v2 = diag_cell->d + profile->w_sub(positiony, positionx);
459    v3 = left_cell->d + profile->w_ins(positionx, positiony);
460
461    akt_cell->d = minimum3(v1, v2, v3);
462
463    if (v1 == akt_cell->d)
464        path_map->set(pos_x, pos_y, ALI_UP);
465    if (v2 == akt_cell->d)
466        path_map->set(pos_x, pos_y, ALI_DIAG);
467    if (v3 == akt_cell->d)
468        path_map->set(pos_x, pos_y, ALI_LEFT);
469}
470
471void            ALI_PREALIGNER::
472calculate_last_column(
473                      ali_prealigner_column * prev_col, ali_prealigner_column * akt_col,
474                      unsigned long pos_x)
475{
476    unsigned long   pos_y;
477
478    calculate_last_column_first_cell(&(*prev_col->cells)[0],
479                                     &(*akt_col->cells)[0], pos_x);
480
481    for (pos_y = 1; pos_y < akt_col->column_length; pos_y++)
482        calculate_last_column_cell(&(*prev_col->cells)[pos_y - 1],
483                                   &(*prev_col->cells)[pos_y],
484                                   &(*akt_col->cells)[pos_y - 1], &(*akt_col->cells)[pos_y],
485                                   pos_x, pos_y);
486}
487
488
489void            ALI_PREALIGNER::
490calculate_matrix(void)
491{
492    unsigned long   pos_x;
493    ali_prealigner_column *akt_col, *prev_col, *h_col;
494
495    akt_col = new ali_prealigner_column(end_y - start_y + 1);
496    prev_col = new ali_prealigner_column(end_y - start_y + 1);
497
498    calculate_first_column(prev_col);
499
500    for (pos_x = 1; pos_x < end_x - start_x + 1; pos_x++) {
501        if (pos_x == end_x - start_x || profile->is_internal_last(pos_x))
502            calculate_last_column(prev_col, akt_col, pos_x);
503        else
504            calculate_column(prev_col, akt_col, pos_x);
505        h_col = akt_col;
506        akt_col = prev_col;
507        prev_col = h_col;
508    }
509
510    delete          akt_col;
511    delete          prev_col;
512}
513
514
515#if 0
516
517int             ALI_PREALIGNER::
518find_multiple_path(
519                   unsigned long start_x, unsigned long start_y,
520                   unsigned long *end_x, unsigned long *end_y)
521{
522    unsigned char   value;
523    long            x1, y1, x2, y2;
524
525    value = path_map->get_value(start_x, start_y);
526
527    /*
528     * (x1,y1) is moving left, diag, up
529     */
530    if (value & ALI_LEFT) {
531        x1 = start_x - 1;
532        y1 = start_y;
533    } else {
534        if (value & ALI_DIAG) {
535            x1 = start_x - 1;
536            y1 = start_y - 1;
537        } else {
538            *end_x = start_x;
539            *end_y = start_y - 1;
540            return 1;
541        }
542    }
543
544    /*
545     * (x2,y2) is moving up, diag, left
546     */
547    if (value & ALI_UP) {
548        x2 = start_x;
549        y2 = start_y - 1;
550    } else {
551        if (value & ALI_DIAG) {
552            x2 = start_x - 1;
553            y2 = start_y - 1;
554        } else {
555            *end_x = start_x - 1;
556            *end_y = start_y;
557            return 1;
558        }
559    }
560
561    while (x1 > 0 && x2 > 0 && y1 > 0 && y2 > 0 &&
562           (x1 != x2 || y1 != y2)) {
563        /*
564         * move (x2,y2)
565         */
566        if (y2 < y1) {
567            value = path_map->get_value((unsigned long) x1, (unsigned long) y1);
568            if (value & ALI_LEFT) {
569                x1--;
570            } else {
571                if (value & ALI_DIAG) {
572                    x1--;
573                    y1--;
574                } else {
575                    y1--;
576                }
577            }
578        }
579        /*
580         * move (x1,y2)
581         */
582        else {
583            value = path_map->get_value((unsigned long) x2, (unsigned long) y2);
584            if (value & ALI_UP) {
585                y2--;
586            } else {
587                if (value & ALI_DIAG) {
588                    x2--;
589                    y2--;
590                } else {
591                    x2--;
592                }
593            }
594        }
595    }
596
597    if (x1 == x2 && y1 == y2) {
598        *end_x = x1;
599        *end_y = y1;
600        return 1;
601    }
602    return 0;
603}
604#endif
605
606
607/*
608 * Generate a sub_solution by deleting all undefined segments
609 */
610void            ALI_PREALIGNER::
611generate_solution(ALI_MAP * map)
612{
613    ALI_MAP        *seg_map;
614    unsigned long   map_pos, map_len;
615    unsigned long   start_seg, end_seg, pos_seg;
616
617    sub_solution = new ALI_SUB_SOLUTION(profile);
618
619    map_len = map->last_base() - map->first_base() + 1;
620    for (map_pos = map->first_base(); map_pos <= map->last_base(); map_pos++) {
621        /*
622         * search for segment
623         */
624        for (start_seg = map_pos; start_seg <= map->last_base() &&
625                 map->is_undefined(start_seg); start_seg++);
626        if (start_seg <= map->last_base()) {
627            for (end_seg = start_seg; end_seg <= map->last_base() &&
628                     (!map->is_undefined(end_seg)); end_seg++);
629            end_seg--;
630
631            seg_map = new ALI_MAP(start_seg,
632                                  end_seg,
633                                  map->position(start_seg),
634                                  map->position(end_seg));
635            /*
636             * Copy segment
637             */
638            for (pos_seg = start_seg; pos_seg <= end_seg; pos_seg++) {
639                if (map->is_inserted(pos_seg))
640                    seg_map->set(pos_seg,
641                                 map->position(pos_seg) - map->position(start_seg), 1);
642                else
643                    seg_map->set(pos_seg,
644                                 map->position(pos_seg) - map->position(start_seg), 0);
645            }
646
647            if (sub_solution->insert(seg_map) != 1)
648                ali_fatal_error("Inconsistent solution?",
649                                "ALI_PREALIGNER::generate_solution()");
650
651            map_pos = end_seg;
652        }
653    }
654}
655
656/*
657 * generate the result mask from an stack of operations
658 */
659void            ALI_PREALIGNER::
660generate_result_mask(ALI_TSTACK < unsigned char >*stack)
661{
662    ALI_SEQUENCE   *seq;
663    float           cost_of_bindings;
664    ALI_MAP        *map;
665    unsigned long   seq_pos, dest_pos;
666    long            i;
667
668    map = new ALI_MAP(start_x, end_x, start_y, end_y);
669
670    seq_pos = start_x;
671    dest_pos = 0;
672    for (i = (long) stack->akt_size() - 1; i >= 0; i--) {
673        switch (stack->get(i)) {
674            case ALI_PREALIGNER_INS:
675                map->set(seq_pos++, dest_pos, 1);
676                break;
677            case ALI_PREALIGNER_SUB:
678                map->set(seq_pos++, dest_pos++, 0);
679                break;
680            case ALI_PREALIGNER_DEL:
681                dest_pos++;
682                break;
683            case ALI_PREALIGNER_INS | ALI_PREALIGNER_MULTI_FLAG:
684                map->set(seq_pos, dest_pos, 1);
685                map->undefine(seq_pos++);
686                break;
687            case ALI_PREALIGNER_SUB | ALI_PREALIGNER_MULTI_FLAG:
688                map->set(seq_pos, dest_pos++, 0);
689                map->undefine(seq_pos++);
690                break;
691            case ALI_PREALIGNER_DEL | ALI_PREALIGNER_MULTI_FLAG:
692                dest_pos++;
693                break;
694            default:
695                ali_fatal_error("Unexpected value",
696                                "ALI_PREALIGNER::generate_result_mask()");
697        }
698    }
699
700    if (result_mask_counter > 0)
701        result_mask_counter--;
702
703    seq = map->sequence_without_inserts(profile->sequence());
704    cost_of_bindings = profile->w_binding(map->first_reference_base(), seq);
705    delete          seq;
706
707    /*
708     * make the intersection
709     */
710    result_mask.insert(map, cost_of_bindings);
711}
712
713/*
714 * Fill the stack with rest DELs or INSs
715 */
716void            ALI_PREALIGNER::
717mapper_post(ALI_TSTACK < unsigned char >*stack,
718            unsigned long ins_nu, unsigned long del_nu)
719{
720    if (ins_nu > 0 && del_nu > 0)
721        ali_fatal_error("Unexpected values",
722                        "ALI_PREALIGNER::mapper_post()");
723
724    if (ins_nu > 0) {
725        stack->push(ALI_PREALIGNER_INS, ins_nu);
726        generate_result_mask(stack);
727        stack->pop(ins_nu);
728    } else {
729        if (del_nu > 0) {
730            stack->push(ALI_PREALIGNER_DEL, del_nu);
731            generate_result_mask(stack);
732            stack->pop(del_nu);
733        } else
734            generate_result_mask(stack);
735    }
736}
737
738/*
739 * Fill the stack with rest DELs or INSs (with MULTI_FLAG)
740 */
741void            ALI_PREALIGNER::
742mapper_post_multi(ALI_TSTACK < unsigned char >*stack,
743                  unsigned long ins_nu, unsigned long del_nu)
744{
745    if (ins_nu > 0 && del_nu > 0)
746        ali_fatal_error("Unexpected values",
747                        "ALI_PREALIGNER::mapper_post_multi()");
748
749    if (ins_nu > 0) {
750        stack->push(ALI_PREALIGNER_INS | ALI_PREALIGNER_MULTI_FLAG, ins_nu);
751        generate_result_mask(stack);
752        stack->pop(ins_nu);
753    } else {
754        if (del_nu > 0) {
755            stack->push(ALI_PREALIGNER_DEL | ALI_PREALIGNER_MULTI_FLAG, del_nu);
756            generate_result_mask(stack);
757            stack->pop(del_nu);
758        } else
759            generate_result_mask(stack);
760    }
761}
762
763/*
764 * generate a stack of operations by taking a random path of the pathmap
765 */
766void            ALI_PREALIGNER::
767mapper_random(ALI_TSTACK < unsigned char >*stack,
768              unsigned long pos_x, unsigned long pos_y)
769{
770    unsigned long   next_x, next_y;
771    unsigned long   random;
772    unsigned char   value;
773    unsigned long   stack_counter = 0;
774
775
776
777    next_x = pos_x;
778    next_y = pos_y;
779    while (next_x <= pos_x && next_y <= pos_y) {
780        stack_counter++;
781
782        random = GB_random(6);
783
784        value = path_map->get_value(next_x, next_y);
785        if (value == 0)
786            ali_fatal_error("Unexpected value (1)",
787                            "ALI_PREALIGNER::mapper_random()");
788
789        switch (random) {
790            case 0:
791                if (value & ALI_UP) {
792                    stack->push(ALI_PREALIGNER_DEL);
793                    next_y--;
794                } else {
795                    if (value & ALI_DIAG) {
796                        stack->push(ALI_PREALIGNER_SUB);
797                        next_x--;
798                        next_y--;
799                    } else {
800                        stack->push(ALI_PREALIGNER_INS);
801                        next_x--;
802                    }
803                }
804                break;
805            case 1:
806                if (value & ALI_UP) {
807                    stack->push(ALI_PREALIGNER_DEL);
808                    next_y--;
809                } else {
810                    if (value & ALI_LEFT) {
811                        stack->push(ALI_PREALIGNER_INS);
812                        next_x--;
813                    } else {
814                        stack->push(ALI_PREALIGNER_SUB);
815                        next_x--;
816                        next_y--;
817                    }
818                }
819                break;
820            case 2:
821                if (value & ALI_DIAG) {
822                    stack->push(ALI_PREALIGNER_SUB);
823                    next_x--;
824                    next_y--;
825                } else {
826                    if (value & ALI_UP) {
827                        stack->push(ALI_PREALIGNER_DEL);
828                        next_y--;
829                    } else {
830                        stack->push(ALI_PREALIGNER_INS);
831                        next_x--;
832                    }
833                }
834                break;
835            case 3:
836                if (value & ALI_DIAG) {
837                    stack->push(ALI_PREALIGNER_SUB);
838                    next_x--;
839                    next_y--;
840                } else {
841                    if (value & ALI_LEFT) {
842                        stack->push(ALI_PREALIGNER_INS);
843                        next_x--;
844                    } else {
845                        stack->push(ALI_PREALIGNER_DEL);
846                        next_y--;
847                    }
848                }
849                break;
850            case 4:
851                if (value & ALI_LEFT) {
852                    stack->push(ALI_PREALIGNER_INS);
853                    next_x--;
854                } else {
855                    if (value & ALI_UP) {
856                        stack->push(ALI_PREALIGNER_DEL);
857                        next_y--;
858                    } else {
859                        stack->push(ALI_PREALIGNER_SUB);
860                        next_x--;
861                        next_y--;
862                    }
863                }
864                break;
865            case 5:
866                if (value & ALI_LEFT) {
867                    stack->push(ALI_PREALIGNER_INS);
868                    next_x--;
869                } else {
870                    if (value & ALI_DIAG) {
871                        stack->push(ALI_PREALIGNER_SUB);
872                        next_x--;
873                        next_y--;
874                    } else {
875                        stack->push(ALI_PREALIGNER_DEL);
876                        next_y--;
877                    }
878                }
879                break;
880            default:
881                ali_fatal_error("Unexpected random value",
882                                "ALI_PREALIGNER::mapper_random()");
883        }
884    }
885
886    if (next_x <= pos_x) {
887        mapper_post(stack, next_x + 1, 0);
888    } else {
889        if (next_y <= pos_y) {
890            mapper_post(stack, 0, next_y + 1);
891        } else {
892            mapper_post(stack, 0, 0);
893        }
894    }
895
896    if (stack_counter > 0)
897        stack->pop(stack_counter);
898}
899
900/*
901 * generate a stack of operations by taking every path
902 */
903void            ALI_PREALIGNER::
904mapper(ALI_TSTACK < unsigned char >*stack,
905       unsigned long pos_x, unsigned long pos_y)
906{
907    unsigned char   value;
908    unsigned long   stack_counter = 0;
909
910    value = path_map->get_value(pos_x, pos_y);
911
912    if (pos_x == 0 || pos_y == 0) {
913        if (value & ALI_UP) {
914            stack->push(ALI_PREALIGNER_DEL);
915            if (pos_y == 0)
916                mapper_post(stack, pos_x + 1, 0);
917            else
918                mapper(stack, pos_x, pos_y - 1);
919            stack->pop();
920        }
921        if (value & ALI_DIAG) {
922            stack->push(ALI_PREALIGNER_SUB);
923            if (pos_y > 0) {
924                mapper_post(stack, 0, pos_y);
925            } else {
926                if (pos_x > 0) {
927                    mapper_post(stack, pos_x, 0);
928                } else {
929                    mapper_post(stack, 0, 0);
930                }
931            }
932            stack->pop();
933        }
934        if (value & ALI_LEFT) {
935            stack->push(ALI_PREALIGNER_INS);
936            if (pos_x == 0)
937                mapper_post(stack, 0, pos_y + 1);
938            else
939                mapper(stack, pos_x - 1, pos_y);
940            stack->pop();
941        }
942        return;
943    }
944    /*
945     * follow an unique path
946     */
947    while (value == ALI_UP || value == ALI_DIAG || value == ALI_LEFT) {
948        stack_counter++;
949        switch (value) {
950            case ALI_UP:
951                stack->push(ALI_PREALIGNER_DEL);
952                pos_y--;
953                break;
954            case ALI_DIAG:
955                stack->push(ALI_PREALIGNER_SUB);
956                pos_x--;
957                pos_y--;
958                break;
959            case ALI_LEFT:
960                stack->push(ALI_PREALIGNER_INS);
961                pos_x--;
962                break;
963        }
964
965        value = path_map->get_value(pos_x, pos_y);
966
967        if (pos_x == 0 || pos_y == 0) {
968            if (value & ALI_UP) {
969                stack->push(ALI_PREALIGNER_DEL);
970                if (pos_y == 0)
971                    mapper_post(stack, pos_x + 1, 0);
972                else
973                    mapper(stack, pos_x, pos_y - 1);
974                stack->pop();
975            }
976            if (value & ALI_DIAG) {
977                stack->push(ALI_PREALIGNER_SUB);
978                if (pos_y > 0) {
979                    mapper_post(stack, 0, pos_y);
980                } else {
981                    if (pos_x > 0) {
982                        mapper_post(stack, pos_x, 0);
983                    } else {
984                        mapper_post(stack, 0, 0);
985                    }
986                }
987                stack->pop();
988            }
989            if (value & ALI_LEFT) {
990                stack->push(ALI_PREALIGNER_INS);
991                if (pos_x == 0)
992                    mapper_post(stack, 0, pos_y + 1);
993                else
994                    mapper(stack, pos_x - 1, pos_y);
995                stack->pop();
996            }
997            if (stack_counter > 0)
998                stack->pop(stack_counter);
999
1000            return;
1001        }
1002    }
1003
1004    if (value & ALI_UP) {
1005        stack->push(ALI_PREALIGNER_DEL);
1006        mapper(stack, pos_x, pos_y - 1);
1007        stack->pop();
1008    }
1009    if (value & ALI_DIAG) {
1010        stack->push(ALI_PREALIGNER_SUB);
1011        mapper(stack, pos_x - 1, pos_y - 1);
1012        stack->pop();
1013    }
1014    if (value & ALI_LEFT) {
1015        stack->push(ALI_PREALIGNER_INS);
1016        mapper(stack, pos_x - 1, pos_y);
1017        stack->pop();
1018    }
1019    if (stack_counter > 0)
1020        stack->pop(stack_counter);
1021}
1022
1023#if 0
1024
1025void            ALI_PREALIGNER::
1026quick_mapper(ALI_TSTACK < unsigned char >*stack,
1027             unsigned long pos_x, unsigned long pos_y)
1028{
1029    unsigned long   end_x, end_y;
1030    unsigned char   value;
1031
1032
1033    value = path_map->get_value(pos_x, pos_y);
1034
1035    while (pos_x > 0 && pos_y > 0) {
1036        if (value == ALI_UP || value == ALI_DIAG || value == ALI_LEFT) {
1037            switch (value) {
1038                case ALI_UP:
1039                    stack->push(ALI_PREALIGNER_DEL);
1040                    pos_y--;
1041                    break;
1042                case ALI_DIAG:
1043                    stack->push(ALI_PREALIGNER_SUB);
1044                    pos_x--;
1045                    pos_y--;
1046                    break;
1047                case ALI_LEFT:
1048                    stack->push(ALI_PREALIGNER_INS);
1049                    pos_x--;
1050                    break;
1051            }
1052            value = path_map->get_value(pos_x, pos_y);
1053        } else {
1054            if (find_multiple_path(pos_x, pos_y, &end_x, &end_y) == 0)
1055                end_x = end_y = 0;
1056            while (pos_x != end_x || pos_y != end_y) {
1057                if (value & ALI_UP) {
1058                    stack->push(ALI_PREALIGNER_DEL | ALI_PREALIGNER_MULTI_FLAG);
1059                    pos_y--;
1060                } else {
1061                    if (value & ALI_DIAG) {
1062                        stack->push(ALI_PREALIGNER_SUB | ALI_PREALIGNER_MULTI_FLAG);
1063                        pos_x--;
1064                        pos_y--;
1065                    } else {
1066                        stack->push(ALI_PREALIGNER_INS | ALI_PREALIGNER_MULTI_FLAG);
1067                        pos_x--;
1068                    }
1069                }
1070                value = path_map->get_value(pos_x, pos_y);
1071            }
1072        }
1073    }
1074
1075    if (pos_x == 0) {
1076        while (pos_y > 0 && value == ALI_UP) {
1077            stack->push(ALI_PREALIGNER_DEL);
1078            pos_y--;
1079            value = path_map->get_value(pos_x, pos_y);
1080        }
1081        switch (value) {
1082            case ALI_UP:
1083                stack->push(ALI_PREALIGNER_DEL);
1084                mapper_post(stack, pos_x, pos_y);
1085                break;
1086            case ALI_DIAG:
1087                stack->push(ALI_PREALIGNER_SUB);
1088                mapper_post(stack, pos_x, pos_y);
1089                break;
1090            case ALI_LEFT:
1091                stack->push(ALI_PREALIGNER_INS);
1092                mapper_post(stack, pos_x, pos_y);
1093                break;
1094            default:
1095                stack->push(ALI_PREALIGNER_SUB | ALI_PREALIGNER_MULTI_FLAG);
1096                mapper_post_multi(stack, pos_x, pos_y);
1097        }
1098    } else {
1099        while (pos_x > 0 && value == ALI_LEFT) {
1100            stack->push(ALI_PREALIGNER_INS);
1101            pos_x--;
1102            value = path_map->get_value(pos_x, pos_y);
1103        }
1104        switch (value) {
1105            case ALI_UP:
1106                stack->push(ALI_PREALIGNER_DEL);
1107                mapper_post(stack, pos_x, pos_y);
1108                break;
1109            case ALI_DIAG:
1110                stack->push(ALI_PREALIGNER_SUB);
1111                mapper_post(stack, pos_x, pos_y);
1112                break;
1113            case ALI_LEFT:
1114                stack->push(ALI_PREALIGNER_INS);
1115                mapper_post(stack, pos_x, pos_y);
1116                break;
1117            default:
1118                stack->push(ALI_PREALIGNER_SUB | ALI_PREALIGNER_MULTI_FLAG);
1119                mapper_post_multi(stack, pos_x, pos_y);
1120        }
1121    }
1122}
1123
1124#endif
1125
1126/*
1127 * make the result map from the path matrix
1128 */
1129void            ALI_PREALIGNER::
1130make_map(void)
1131{
1132    unsigned long   number_of_sol;
1133    ALI_TSTACK < unsigned char >*stack;
1134
1135    stack = new ALI_TSTACK < unsigned char >(end_x - start_x + end_y - start_y + 3);
1136    if (stack == 0) ali_fatal_error("Out of memory");
1137
1138    number_of_sol = number_of_solutions();
1139    printf("%lu solutions generated\n", number_of_sol);
1140
1141    if (number_of_sol == 0 || number_of_sol > result_mask_counter) {
1142        ali_message("Starting random mapping");
1143        do {
1144            mapper_random(stack, end_x - start_x, end_y - start_y);
1145        } while (result_mask_counter > 0);
1146    } else {
1147        ali_message("Starting systematic mapping");
1148        mapper(stack, end_x - start_x, end_y - start_y);
1149    }
1150
1151    delete          stack;
1152
1153    ali_message("Mapping finished");
1154}
1155
1156#if 0
1157
1158void            ALI_PREALIGNER::
1159make_quick_map(void)
1160{
1161    ALI_TSTACK < unsigned char >*stack;
1162
1163    ali_message("Starting quick mapping");
1164
1165    stack = new ALI_TSTACK < unsigned char >(end_x - start_x + end_y - start_y + 3);
1166    if (stack == 0)
1167        ali_fatal_error("Out of memory");
1168
1169    quick_mapper(stack, end_x - start_x, end_y - start_y);
1170
1171    delete          stack;
1172
1173    ali_message("Quick mapping finished");
1174}
1175
1176#endif
1177
1178/*
1179 * generate an approximation of a complete solution
1180 */
1181void            ALI_PREALIGNER::
1182generate_approximation(ALI_SUB_SOLUTION * work_sol)
1183{
1184    ALI_MAP        *map;
1185    ALI_SEQUENCE   *seq;
1186    char           *ins_marker;
1187    float           binding_costs;
1188
1189    map = work_sol->make_one_map();
1190    if (map == 0)
1191        ali_fatal_error("Can't make one map",
1192                        "ALI_PREALIGNER::generate_approximation()");
1193
1194    seq = map->sequence_without_inserts(profile->sequence());
1195    binding_costs = profile->w_binding(map->first_base(), seq);
1196    delete          seq;
1197
1198    ins_marker = map->insert_marker();
1199
1200    result_approx.insert(map, ins_marker, binding_costs);
1201    // delete ins_marker;   @@@
1202    // delete map;          @@@
1203}
1204
1205/*
1206 * combine subsolutions for an approximation
1207 */
1208void            ALI_PREALIGNER::
1209mapper_approximation(unsigned long area_no,
1210                     ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>*map_lists,
1211                     ALI_SUB_SOLUTION * work_sol)
1212{
1213    ALI_TLIST < ALI_MAP * >*map_list;
1214    ALI_MAP        *map;
1215
1216    /*
1217     * stop mapping at last area
1218     */
1219    if (area_no > map_lists->size())
1220        return;
1221
1222    if (area_no == map_lists->size()) {
1223        generate_approximation(work_sol);
1224        return;
1225    }
1226    /*
1227     * map area number 'area_no'
1228     */
1229    map_list = map_lists->get(area_no);
1230    if (map_list->is_empty())
1231        ali_fatal_error("Found empty list",
1232                        "ALI_PREALIGNER::mapper_approximation()");
1233
1234    /*
1235     * combine all possibilities
1236     */
1237    map = map_list->first();
1238    do {
1239        if (!work_sol->insert(map))
1240            ali_fatal_error("Can't insert map",
1241                            "ALI_PREALIGNER::mapper_approximation()");
1242
1243        mapper_approximation(area_no + 1, map_lists, work_sol);
1244
1245        if (!work_sol->delete_map(map))
1246            ali_fatal_error("Can't delete map",
1247                            "ALI_PREALIGNER::mapper_approximation()");
1248
1249        if (map_list->is_next())
1250            map = map_list->next();
1251        else
1252            map = 0;
1253    } while (map != 0);
1254}
1255
1256/*
1257 * Make an approximation by aligning the undefined sections
1258 */
1259void            ALI_PREALIGNER::
1260make_approximation(ALI_PREALIGNER_CONTEXT * context)
1261{
1262    ALI_SUB_SOLUTION *work_solution;
1263    ALI_ALIGNER_CONTEXT aligner_context;
1264    ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>*map_lists;
1265    ALI_ALIGNER    *aligner;
1266    unsigned long   area_number;
1267    unsigned long   start_seq, end_seq, start_ref, end_ref;
1268
1269    ali_message("Align free areas");
1270
1271    work_solution = new ALI_SUB_SOLUTION(sub_solution);
1272
1273    aligner_context.max_number_of_maps = context->max_number_of_maps_aligner;
1274
1275    area_number = sub_solution->number_of_free_areas();
1276    printf("number of areas = %ld (Maximal %ld solutions)\n", area_number,
1277           aligner_context.max_number_of_maps);
1278
1279    map_lists = new ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>(area_number);
1280
1281    /*
1282     * generate Solutions for all free areas
1283     */
1284    area_number = 0;
1285    while (sub_solution->free_area(&start_seq, &end_seq, &start_ref, &end_ref,
1286                                   area_number)) {
1287        printf("aligning area %ld (%ld,%ld) - (%ld,%ld)\n",
1288               area_number, start_seq, end_seq, start_ref, end_ref);
1289
1290        aligner = new ALI_ALIGNER(&aligner_context, profile,
1291                                  start_seq, end_seq, start_ref, end_ref);
1292        map_lists->set(area_number, aligner->solutions());
1293
1294        printf("%d solutions generated\n",
1295               map_lists->get(area_number)->cardinality());
1296
1297        delete          aligner;
1298        area_number++;
1299    }
1300
1301    /*
1302     * combine and evaluate the solutions
1303     */
1304    mapper_approximation(0, map_lists, work_solution);
1305
1306    delete          work_solution;
1307    // delete               map_lists;      // @@@
1308
1309    ali_message("Free areas aligned");
1310}
1311
1312
1313/*
1314 * approximate the number of solutions in the pathmap
1315 */
1316unsigned long   ALI_PREALIGNER::
1317number_of_solutions(void)
1318{
1319#define INFINIT 1000000
1320#define ADD(a,b)  if (a>=INFINIT || b>=INFINIT) {a = INFINIT;} else {a += b;}
1321
1322    unsigned long   pos_x, pos_y, col_length;
1323    unsigned long   number;
1324    unsigned char   value;
1325    unsigned long  *column1, *column2, *elem_akt_col, *elem_left_col;
1326
1327    col_length = end_y - start_y + 1;
1328    column1 = (unsigned long *) CALLOC((unsigned int) col_length,
1329                                       sizeof(unsigned long));
1330    column2 = (unsigned long *) CALLOC((unsigned int) col_length,
1331                                       sizeof(unsigned long));
1332    if (column1 == 0 || column2 == 0)
1333        ali_fatal_error("Out of memory");
1334
1335    ali_message("Start: Checking number of solutions");
1336
1337    if (end_x - (start_x & 0x01)) {
1338        elem_akt_col = column1 + col_length - 1;
1339        elem_left_col = column2 + col_length - 1;
1340    } else {
1341        elem_akt_col = column2 + col_length - 1;
1342        elem_left_col = column1 + col_length - 1;
1343    }
1344
1345    number = 0;
1346    *elem_akt_col = 1;
1347    for (pos_x = end_x - start_x; pos_x > 0;) {
1348        *(elem_left_col) = 0;
1349        for (pos_y = end_y - start_y; pos_y > 0; pos_y--) {
1350            *(elem_left_col - 1) = 0;
1351            value = path_map->get_value(pos_x, pos_y);
1352            if (value & ALI_UP) {
1353                /* *(elem_akt_col - 1) += *elem_akt_col; */
1354                ADD(*(elem_akt_col - 1), *elem_akt_col);
1355            }
1356            if (value & ALI_DIAG) {
1357                /* *(elem_left_col - 1) += *elem_akt_col; */
1358                ADD(*(elem_left_col - 1), *elem_akt_col);
1359            }
1360            if (value & ALI_LEFT) {
1361                /* *(elem_left_col) += *elem_akt_col; */
1362                ADD(*(elem_left_col), *elem_akt_col);
1363            }
1364            elem_akt_col--;
1365            elem_left_col--;
1366        }
1367        value = path_map->get_value(pos_x, 0);
1368        if (value & ALI_UP) {
1369            /* number += *elem_akt_col; */
1370            ADD(number, *elem_akt_col);
1371        }
1372        if (value & ALI_DIAG) {
1373            /* number += *elem_akt_col; */
1374            ADD(number, *elem_akt_col);
1375        }
1376        if (value & ALI_LEFT) {
1377            /* *(elem_left_col) += *elem_akt_col; */
1378            ADD(*(elem_left_col), *elem_akt_col);
1379        }
1380        pos_x--;
1381        /*
1382         * toggle the columns
1383         */
1384        if (pos_x & 0x01) {
1385            elem_akt_col = column1 + col_length - 1;
1386            elem_left_col = column2 + col_length - 1;
1387        } else {
1388            elem_akt_col = column2 + col_length - 1;
1389            elem_left_col = column1 + col_length - 1;
1390        }
1391    }
1392
1393    for (pos_y = end_y - start_y; pos_y > 0; pos_y--) {
1394        value = path_map->get_value(0, pos_y);
1395        if (value & ALI_UP) {
1396            /* *(elem_akt_col - 1) += *elem_akt_col; */
1397            ADD(*(elem_akt_col - 1), *elem_akt_col);
1398        }
1399        if (value & ALI_DIAG) {
1400            /* number += *elem_akt_col; */
1401            ADD(number, *elem_akt_col);
1402        }
1403        if (value & ALI_LEFT) {
1404            /* number += *elem_akt_col; */
1405            ADD(number, *elem_akt_col);
1406        }
1407        elem_akt_col--;
1408    }
1409
1410    /* number += *elem_akt_col; */
1411    ADD(number, *elem_akt_col);
1412
1413    ali_message("End: Checking number of solutions");
1414
1415    free((char *) column1);
1416    free((char *) column2);
1417
1418    return number;
1419}
1420
1421
1422/*****************************************************************************
1423 *
1424 * class ALI_PREALIGNER (public)
1425 *
1426 *****************************************************************************/
1427
1428ALI_PREALIGNER::ALI_PREALIGNER(ALI_PREALIGNER_CONTEXT * context,
1429                               ALI_PROFILE * prof,
1430                               unsigned long sx, unsigned long ex,
1431                               unsigned long sy, unsigned long ey)
1432{
1433    profile = prof;
1434
1435    start_x = sx;
1436    end_x = ex;
1437    start_y = sy;
1438    end_y = ey;
1439
1440    result_mask_counter = context->max_number_of_maps;
1441
1442    ali_message("Prealigning");
1443
1444    path_map = new ALI_PATHMAP(end_x - start_x + 1, end_y - start_y + 1);
1445
1446    calculate_matrix();
1447
1448    make_map();
1449
1450    result_mask.delete_expensive(context, profile);
1451    delete          path_map;
1452
1453    generate_solution(result_mask.map);
1454
1455    make_approximation(context);
1456
1457    ali_message("Prealigning finished");
1458}
Note: See TracBrowser for help on using the repository browser.