ARB
ali_aligner.cxx
Go to the documentation of this file.
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 
16 float ali_aligner_dellist::update(unsigned long position) {
17  // increase all multi gaps and update the apropriate costs
18  float minimal_costs = 0.0;
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 
37 ALI_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;
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 
79 void 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;
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()) {
115  akt = list_of_dels.current();
116  }
117  else {
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  }
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 
141 inline float ALI_ALIGNER::minimum2(float v1, float v2) {
142  return (v1 < v2) ? v1 : v2;
143 }
144 
145 inline float ALI_ALIGNER::minimum3(float v1, float v2, float v3) {
146  return (v1 < v2) ? ((v1 < v3) ? v1 : v3) : ((v2 < v3) ? v2 : v3);
147 }
148 
149 inline 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 
160 inline 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 
202 void 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 
219 void 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 
251 void 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 
318 void 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 
340 void 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 
368 void 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 
427 void 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 
504 void 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 
526 void 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;
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 
699 void 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;
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 
799 void 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;
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 
821 void 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;
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 
844 void 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 
947 void ALI_ALIGNER::make_map_systematic(ALI_TSTACK<unsigned char> *stack) {
948  // find ALL paths and make the apropriate maps
949  float min;
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 
985 void ALI_ALIGNER::make_map() {
986  // make the list of result maps
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 
1022  unsigned long v1, v2, v3;
1023 };
1024 
1025 unsigned 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;
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 
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 }
ALI_TLIST< ali_aligner_dellist_elem * > list_of_dels
Definition: ali_aligner.hxx:94
void ali_out_of_memory_if(bool cond)
Definition: ali_misc.hxx:52
ALI_TLIST< ali_pathmap_up_pointer > up_starts
float update(unsigned long position)
Definition: ali_aligner.cxx:16
int has_next()
Definition: ali_tlist.hxx:168
long
Definition: AW_awar.cxx:154
void optimize(unsigned long position)
Definition: ali_aligner.cxx:79
void get(unsigned long x, unsigned long y, unsigned char *val, ALI_TARRAY< ali_pathmap_up_pointer > **up_pointer)
Definition: ali_pathmap.cxx:61
float w_ins(unsigned long, unsigned long)
float w_del(unsigned long begin, unsigned long end)
ali_aligner_cell ** cells
Definition: ali_aligner.hxx:51
unsigned long v3
#define ALI_DIAG
Definition: ali_pathmap.hxx:20
int cardinality()
Definition: ali_tlist.hxx:162
void delete_element()
Definition: ali_tlist.hxx:355
void set(unsigned long position, T value)
Definition: ali_tarray.hxx:49
float w_del_multi_cheap(unsigned long start, unsigned long end)
void set(unsigned long x, unsigned long y, unsigned char val, ALI_TARRAY< ali_pathmap_up_pointer > *up_pointer=NULp)
Definition: ali_pathmap.cxx:36
void push(T value, unsigned long count=1)
Definition: ali_tstack.hxx:36
#define ALI_ALIGNER_SUB
Definition: ali_aligner.hxx:25
void * CALLOC(long i, long j)
Definition: ali_misc.hxx:56
ALI_PROFILE * profile
Definition: ali_aligner.hxx:93
unsigned char operation
Definition: ali_pathmap.hxx:27
ALI_TARRAY< ali_pathmap_up_pointer > * starts(float costs, unsigned long y_offset)
Definition: ali_aligner.cxx:37
#define ALI_ALIGNER_DEL
Definition: ali_aligner.hxx:26
float w_ins_cheap(unsigned long, unsigned long)
#define ALI_ALIGNER_INS
Definition: ali_aligner.hxx:24
unsigned long akt_size()
Definition: ali_tstack.hxx:33
void ali_error(const char *message, const char *func)
Definition: ali_main.cxx:90
#define ALI_LEFT
Definition: ali_pathmap.hxx:19
void optimize(unsigned long x)
Definition: ali_pathmap.cxx:92
#define ALI_UP
Definition: ali_pathmap.hxx:21
void mark_element()
Definition: ali_tlist.hxx:146
int is_empty()
Definition: ali_tlist.hxx:165
T get(unsigned long position)
Definition: ali_tarray.hxx:54
void update_border(unsigned long start_x, unsigned long end_x, unsigned long start_y, unsigned long end_y)
void insert(unsigned long start, float costs, unsigned char operation)
GB_write_int const char GB_write_autoconv_string WRITE_SKELETON(write_pointer, GBDATA *,"%p", GB_write_pointer) char *AW_awa if)(!gb_var) return strdup("")
Definition: AW_awar.cxx:166
int GB_random(int range)
Definition: admath.cxx:88
void update_left(ali_aligner_cell *akt_cell, unsigned long akt_pos, unsigned long start_x, unsigned long end_x)
float w_del_cheap(unsigned long position)
ALI_ALIGNER(ALI_ALIGNER_CONTEXT *context, ALI_PROFILE *profile, unsigned long start_x, unsigned long end_x, unsigned long start_y, unsigned long end_y)
void set(unsigned long base, unsigned long pos, int insert=-1)
void marked()
Definition: ali_tlist.hxx:149
float gap_percent(unsigned long begin, unsigned long end)
void ali_fatal_error(const char *message, const char *func)
Definition: ali_main.cxx:85
unsigned long v2
float w_sub(unsigned long positiony, unsigned long positionx)
void ali_message(const char *message)
Definition: ali_misc.hxx:46
float w_ins_multi_cheap(unsigned long startx, unsigned long endx)
#define NULp
Definition: cxxforward.h:97
T pop(unsigned long count=1)
Definition: ali_tstack.hxx:42
unsigned long v1
void insert(ALI_MAP *in_map)
void update_up(ali_aligner_cell *akt_cell, unsigned long akt_pos, unsigned long start_y, unsigned long end_y)
#define ALI_LUP
Definition: ali_pathmap.hxx:22
unsigned long last_base()
ALI_TLIST< ali_pathmap_up_pointer > left_starts
unsigned long size()
Definition: ali_tarray.hxx:46
T get(unsigned long position)
Definition: ali_tstack.hxx:55
#define min(a, b)
Definition: f2c.h:153
unsigned long column_length
Definition: ali_aligner.hxx:50