source: branches/stable/NALIGNER/ali_aligner.cxx

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