source: branches/port5/NALIGNER/ali_aligner.cxx

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