Context Navigation

source: branches/port5/NALIGNER/ali_prealigner.cxx

Visit:

Last change on this file was 5630, checked in by westram, 16 years ago
fixed a few gcc 4.3.2 warnings
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 44.4 KB

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

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

Download in other formats: