Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

source: tags/initial/NALIGNER/ali_aligner.cxx

Visit:

Last change on this file was 2, checked in by oldcode, 24 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 34.2 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: