Context Navigation

source: tags/arb-6.0/NALIGNER/ali_aligner.cxx

Visit:

Last change on this file was 6659, checked in by westram, 14 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
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;
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
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;
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
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;
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
141	inline float ALI_ALIGNER::minimum2(float v1, float v2)
142	{
143	return (v1 < v2) ? v1 : v2;
144	}
145
146	inline 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
151	inline 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
163	inline 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
205	void 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
223	void 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
255	void 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
322	void 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
344	void 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
373	void 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
433	void 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
510	void 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
532	void 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
705	void 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
805	void 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
827	void 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
850	void 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
953	void 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
991	void 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
1027	struct ali_aligner_tripel {
1028	unsigned long v1, v2, v3;
1029	};
1030
1031	unsigned 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
1288	ALI_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.

Download in other formats: