source: branches/lib/NALIGNER/ali_aligner.cxx

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