// =============================================================== // // // // File : ali_prealigner.cxx // // Purpose : // // // // Institute of Microbiology (Technical University Munich) // // http://www.arb-home.de/ // // // // =============================================================== // #include "ali_prealigner.hxx" #include "ali_aligner.hxx" void ali_prealigner_mask::insert(ALI_MAP * in_map, float costs) { // Insert a new map unsigned long i; calls++; if (!map) { map = in_map; cost_of_binding = costs; } else { if (costs > cost_of_binding) { delete in_map; return; } if (map->first_base() != in_map->first_base() || map->last_base() != in_map->last_base() || map->first_reference_base() != in_map->first_reference_base() || map->last_reference_base() != in_map->last_reference_base()) { ali_fatal_error("Incompatible maps", "ali_prealigner_mask::insert()"); } if (costs < cost_of_binding) { delete map; map = in_map; cost_of_binding = costs; last_new = calls; last_joins = 0; return; } joins++; last_joins++; for (i = map->first_base(); i <= map->last_base(); i++) if ((!map->is_undefined(i)) && map->position(i) != in_map->position(i)) map->undefine(i); delete in_map; } } void ali_prealigner_mask::delete_expensive(ALI_PREALIGNER_CONTEXT * context, ALI_PROFILE * profile) { // Delete expensive parts of solution ALI_MAP *inverse_map; unsigned long start_hel, end_hel; unsigned long start_seq, end_seq; unsigned long start_mapped, end_mapped; unsigned long start_ok=0, end_ok=0; int start_ok_flag; unsigned long found_helix; unsigned long error_counter; unsigned long map_pos, i, j; float max_cost, helix_cost; unsigned long helix_counter; long compl_pos; unsigned char base1, base2; printf("MASK : calls = %ld joins = %ld last_new = %ld last_joins = %ld\n", calls, joins, last_new, last_joins); max_cost = profile->w_sub_maximum() * context->max_cost_of_sub_percent; // Delete expensive Bases error_counter = 0; for (i = map->first_base(); i <= map->last_base(); i++) { if (!(map->is_inserted(i)) && profile->w_sub(map->position(i), i) > max_cost) { error_counter++; if (error_counter > context->error_count) map->undefine(i); else { if (error_counter == context->error_count) { for (j = i - error_counter + 1; j <= i; j++) map->undefine(j); } } } else { // If error was in helix => delete helix total if (error_counter > 0 && profile->is_in_helix(map->position(i - 1), &start_hel, &end_hel)) { for (j = i - 1; map->position(j) >= long(start_hel); j--) ; for (; map->position(j) <= long(end_hel); j++) map->undefine(j); } error_counter = 0; } } // Delete expensive Helizes inverse_map = map->inverse_without_inserts(); for (i = inverse_map->first_base(); i <= inverse_map->last_base(); i++) { // found a helix if (profile->is_in_helix(i, &start_hel, &end_hel)) { if (i != start_hel) ali_fatal_error("Inconsistent positions", "ali_prealigner_mask::delete_expensive()"); compl_pos = profile->complement_position(start_hel); // only forward bindings if (compl_pos > long(end_hel)) { helix_cost = 0.0; helix_counter = 0; while (i <= end_hel) { // is binding ? if (compl_pos > 0) { if (!inverse_map->is_undefined(i)) { base1 = (profile->sequence())->base( inverse_map->first_reference_base() + inverse_map->position(i)); } else { base1 = ALI_GAP_CODE; } if (!inverse_map->is_undefined(compl_pos)) { base2 = (profile->sequence())->base( inverse_map->first_reference_base() + inverse_map->position(compl_pos)); } else { base2 = ALI_GAP_CODE; } if (base1 != ALI_GAP_CODE || base2 != ALI_GAP_CODE) { helix_cost += profile->w_bind(i, base1, compl_pos, base2); helix_counter++; } } i++; compl_pos = profile->complement_position(i); } if (helix_counter > 0) helix_cost /= helix_counter; if (helix_cost > context->max_cost_of_helix) { for (j = start_hel; j <= end_hel; j++) { if (!inverse_map->is_undefined(j)) { map->undefine(inverse_map->first_reference_base() + inverse_map->position(j)); } } for (j = profile->complement_position(end_hel); long(j) <= profile->complement_position(start_hel); j++) { if (!inverse_map->is_undefined(j)) { map->undefine(inverse_map->first_reference_base() + inverse_map->position(j)); } } } } i = end_hel; } } delete inverse_map; // Check for good parts for (map_pos = map->first_base(); map_pos <= map->last_base(); map_pos++) { // search next defined segment if (!map->is_undefined(map_pos)) { // find start and end of segment start_seq = map_pos; start_mapped = map->position(map_pos); for (map_pos++; map_pos <= map->last_base() && (!map->is_undefined(map_pos)); map_pos++) ; end_seq = map_pos - 1; end_mapped = map->position(end_seq); // Check segment for helizes found_helix = 0; start_ok_flag = 0; for (i = start_seq; i <= end_seq; i++) { if (profile->is_in_helix(map->position(i), &start_hel, &end_hel)) { found_helix++; // Helix is inside the segment if (start_hel >= start_mapped && end_hel <= end_mapped) { if (start_ok_flag == 0) { start_ok = start_hel; start_ok_flag = 1; } end_ok = end_hel; } } } // Found good helizes if (start_ok_flag == 1) { for (i = start_seq; map->position(i) < long(start_ok); i++) map->undefine(i); for (i = end_seq; map->position(i) > long(end_ok); i--) map->undefine(i); } else { // Found bad helizes if (found_helix > 0) { for (i = start_seq; i <= end_seq; i++) map->undefine(i); } // Segment without helix else { if (end_seq - start_seq + 1 >= (unsigned long)((2 * context->interval_border) + context->interval_center)) { for (i = start_seq; i < start_seq + context->interval_border; i++) map->undefine(i); for (i = end_seq; i > end_seq - context->interval_border; i--) map->undefine(i); } else { for (i = start_seq; i <= end_seq; i++) map->undefine(i); } } } } } } // ----------------------- // ALI_PREALIGNER inline float ALI_PREALIGNER::minimum2(float a, float b) { return (a < b) ? a : b; } inline float ALI_PREALIGNER::minimum3(float a, float b, float c) { return (a < b) ? ((a < c) ? a : c) : ((b < c) ? b : c); } inline void ALI_PREALIGNER::calculate_first_column_first_cell(ali_prealigner_cell * akt_cell) { float v1, v2; v1 = profile->w_ins_multi_cheap(start_x, start_y) + profile->w_sub_multi_gap_cheap(start_y, start_y); v2 = profile->w_sub(start_y, start_x); akt_cell->d = minimum2(v1, v2); if (akt_cell->d == v1) path_map->set(0, 0, ALI_UP | ALI_LEFT); if (akt_cell->d == v2) path_map->set(0, 0, ALI_DIAG); } inline void ALI_PREALIGNER::calculate_first_column_cell(ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell, unsigned long pos_y) { float v1, v2, v3; unsigned long positiony; positiony = start_y + pos_y; v1 = up_cell->d + profile->w_sub_gap(positiony); v2 = profile->w_sub_multi_gap_cheap(start_y, positiony - 1) + profile->w_sub(positiony, start_x); v3 = profile->w_sub_multi_gap_cheap(start_y, positiony) + profile->w_ins(start_x, positiony); akt_cell->d = minimum3(v1, v2, v3); if (v1 == akt_cell->d) path_map->set(0, pos_y, ALI_UP); if (v2 == akt_cell->d) path_map->set(0, pos_y, ALI_DIAG); if (v3 == akt_cell->d) path_map->set(0, pos_y, ALI_LEFT); } void ALI_PREALIGNER::calculate_first_column(ali_prealigner_column * akt_column) { unsigned long pos_y; calculate_first_column_first_cell(&(*akt_column->cells)[0]); for (pos_y = 1; pos_y < akt_column->column_length; pos_y++) calculate_first_column_cell(&(*akt_column->cells)[pos_y - 1], &(*akt_column->cells)[pos_y], pos_y); } inline void ALI_PREALIGNER::calculate_first_cell(ali_prealigner_cell * left_cell, ali_prealigner_cell * akt_cell, unsigned long pos_x) { float v1, v2, v3; unsigned long positionx; positionx = start_x + pos_x; v1 = profile->w_ins_multi_cheap(start_x, positionx) + profile->w_sub_gap(start_y); v2 = profile->w_ins_multi_cheap(start_x, positionx - 1) + profile->w_sub(start_y, positionx); v3 = left_cell->d + profile->w_ins(positionx, start_y); akt_cell->d = minimum3(v1, v2, v3); if (v1 == akt_cell->d) path_map->set(pos_x, 0, ALI_UP); if (v2 == akt_cell->d) path_map->set(pos_x, 0, ALI_DIAG); if (v3 == akt_cell->d) path_map->set(pos_x, 0, ALI_LEFT); } inline void ALI_PREALIGNER::calculate_cell(ali_prealigner_cell * diag_cell, ali_prealigner_cell * left_cell, ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell, unsigned long pos_x, unsigned long pos_y) { float v1, v2, v3; unsigned long positionx, positiony; positionx = start_x + pos_x; positiony = start_y + pos_y; v1 = up_cell->d + profile->w_sub_gap(positiony); v2 = diag_cell->d + profile->w_sub(positiony, positionx); v3 = left_cell->d + profile->w_ins(positionx, positiony); akt_cell->d = minimum3(v1, v2, v3); if (v1 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_UP); if (v2 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_DIAG); if (v3 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_LEFT); } void ALI_PREALIGNER::calculate_column(ali_prealigner_column * prev_col, ali_prealigner_column * akt_col, unsigned long pos_x) { unsigned long pos_y; calculate_first_cell(&(*prev_col->cells)[0], &(*akt_col->cells)[0], pos_x); for (pos_y = 1; pos_y < akt_col->column_length; pos_y++) calculate_cell(&(*prev_col->cells)[pos_y - 1], &(*prev_col->cells)[pos_y], &(*akt_col->cells)[pos_y - 1], &(*akt_col->cells)[pos_y], pos_x, pos_y); } inline void ALI_PREALIGNER::calculate_last_column_first_cell(ali_prealigner_cell * left_cell, ali_prealigner_cell * akt_cell, unsigned long pos_x) { float v1, v2, v3; unsigned long positionx; positionx = start_x + pos_x; v1 = profile->w_ins_multi_cheap(start_x, positionx) + profile->w_sub_gap_cheap(start_y); v2 = profile->w_ins_multi_cheap(start_x, positionx - 1) + profile->w_sub(start_y, positionx); v3 = left_cell->d + profile->w_ins(positionx, start_y); akt_cell->d = minimum3(v1, v2, v3); if (v1 == akt_cell->d) path_map->set(pos_x, 0, ALI_UP); if (v2 == akt_cell->d) path_map->set(pos_x, 0, ALI_DIAG); if (v3 == akt_cell->d) path_map->set(pos_x, 0, ALI_LEFT); } inline void ALI_PREALIGNER::calculate_last_column_cell(ali_prealigner_cell * diag_cell, ali_prealigner_cell * left_cell, ali_prealigner_cell * up_cell, ali_prealigner_cell * akt_cell, unsigned long pos_x, unsigned long pos_y) { float v1, v2, v3; unsigned long positionx, positiony; positionx = start_x + pos_x; positiony = start_y + pos_y; v1 = up_cell->d + profile->w_sub_gap_cheap(positiony); v2 = diag_cell->d + profile->w_sub(positiony, positionx); v3 = left_cell->d + profile->w_ins(positionx, positiony); akt_cell->d = minimum3(v1, v2, v3); if (v1 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_UP); if (v2 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_DIAG); if (v3 == akt_cell->d) path_map->set(pos_x, pos_y, ALI_LEFT); } void ALI_PREALIGNER::calculate_last_column(ali_prealigner_column * prev_col, ali_prealigner_column * akt_col, unsigned long pos_x) { unsigned long pos_y; calculate_last_column_first_cell(&(*prev_col->cells)[0], &(*akt_col->cells)[0], pos_x); for (pos_y = 1; pos_y < akt_col->column_length; pos_y++) calculate_last_column_cell(&(*prev_col->cells)[pos_y - 1], &(*prev_col->cells)[pos_y], &(*akt_col->cells)[pos_y - 1], &(*akt_col->cells)[pos_y], pos_x, pos_y); } void ALI_PREALIGNER::calculate_matrix() { unsigned long pos_x; ali_prealigner_column *akt_col, *prev_col, *h_col; akt_col = new ali_prealigner_column(end_y - start_y + 1); prev_col = new ali_prealigner_column(end_y - start_y + 1); calculate_first_column(prev_col); for (pos_x = 1; pos_x < end_x - start_x + 1; pos_x++) { if (pos_x == end_x - start_x || profile->is_internal_last(pos_x)) calculate_last_column(prev_col, akt_col, pos_x); else calculate_column(prev_col, akt_col, pos_x); h_col = akt_col; akt_col = prev_col; prev_col = h_col; } delete akt_col; delete prev_col; } void ALI_PREALIGNER::generate_solution(ALI_MAP * map) { // Generate a sub_solution by deleting all undefined segments ALI_MAP *seg_map; unsigned long map_pos; unsigned long start_seg, end_seg, pos_seg; sub_solution = new ALI_SUB_SOLUTION(profile); for (map_pos = map->first_base(); map_pos <= map->last_base(); map_pos++) { // search for segment for (start_seg = map_pos; start_seg <= map->last_base() && map->is_undefined(start_seg); start_seg++) ; if (start_seg <= map->last_base()) { for (end_seg = start_seg; end_seg <= map->last_base() && (!map->is_undefined(end_seg)); end_seg++) ; end_seg--; seg_map = new ALI_MAP(start_seg, end_seg, map->position(start_seg), map->position(end_seg)); // Copy segment for (pos_seg = start_seg; pos_seg <= end_seg; pos_seg++) { if (map->is_inserted(pos_seg)) seg_map->set(pos_seg, map->position(pos_seg) - map->position(start_seg), 1); else seg_map->set(pos_seg, map->position(pos_seg) - map->position(start_seg), 0); } if (sub_solution->insert(seg_map) != 1) ali_fatal_error("Inconsistent solution?", "ALI_PREALIGNER::generate_solution()"); map_pos = end_seg; } } } void ALI_PREALIGNER::generate_result_mask(ALI_TSTACK < unsigned char >*stack) { // generate the result mask from an stack of operations ALI_SEQUENCE *seq; float cost_of_bindings; ALI_MAP *map; unsigned long seq_pos, dest_pos; long i; map = new ALI_MAP(start_x, end_x, start_y, end_y); seq_pos = start_x; dest_pos = 0; for (i = (long) stack->akt_size() - 1; i >= 0; i--) { switch (stack->get(i)) { case ALI_PREALIGNER_INS: map->set(seq_pos++, dest_pos, 1); break; case ALI_PREALIGNER_SUB: map->set(seq_pos++, dest_pos++, 0); break; case ALI_PREALIGNER_DEL: dest_pos++; break; case ALI_PREALIGNER_INS | ALI_PREALIGNER_MULTI_FLAG : map->set(seq_pos, dest_pos, 1); map->undefine(seq_pos++); break; case ALI_PREALIGNER_SUB | ALI_PREALIGNER_MULTI_FLAG : map->set(seq_pos, dest_pos++, 0); map->undefine(seq_pos++); break; case ALI_PREALIGNER_DEL | ALI_PREALIGNER_MULTI_FLAG : dest_pos++; break; default: ali_fatal_error("Unexpected value", "ALI_PREALIGNER::generate_result_mask()"); } } if (result_mask_counter > 0) result_mask_counter--; seq = map->sequence_without_inserts(profile->sequence()); cost_of_bindings = profile->w_binding(map->first_reference_base(), seq); delete seq; // make the intersection result_mask.insert(map, cost_of_bindings); } void ALI_PREALIGNER::mapper_post(ALI_TSTACK < unsigned char >*stack, unsigned long ins_nu, unsigned long del_nu) { // Fill the stack with rest DELs or INSs if (ins_nu > 0 && del_nu > 0) ali_fatal_error("Unexpected values", "ALI_PREALIGNER::mapper_post()"); if (ins_nu > 0) { stack->push(ALI_PREALIGNER_INS, ins_nu); generate_result_mask(stack); stack->pop(ins_nu); } else { if (del_nu > 0) { stack->push(ALI_PREALIGNER_DEL, del_nu); generate_result_mask(stack); stack->pop(del_nu); } else generate_result_mask(stack); } } void ALI_PREALIGNER::mapper_post_multi(ALI_TSTACK < unsigned char >*stack, unsigned long ins_nu, unsigned long del_nu) { // Fill the stack with rest DELs or INSs (with MULTI_FLAG) if (ins_nu > 0 && del_nu > 0) ali_fatal_error("Unexpected values", "ALI_PREALIGNER::mapper_post_multi()"); if (ins_nu > 0) { stack->push(ALI_PREALIGNER_INS | ALI_PREALIGNER_MULTI_FLAG, ins_nu); generate_result_mask(stack); stack->pop(ins_nu); } else { if (del_nu > 0) { stack->push(ALI_PREALIGNER_DEL | ALI_PREALIGNER_MULTI_FLAG, del_nu); generate_result_mask(stack); stack->pop(del_nu); } else generate_result_mask(stack); } } void ALI_PREALIGNER::mapper_random(ALI_TSTACK < unsigned char >*stack, unsigned long pos_x, unsigned long pos_y) { // generate a stack of operations by taking a random path of the pathmap unsigned long next_x, next_y; unsigned long random; unsigned char value; unsigned long stack_counter = 0; next_x = pos_x; next_y = pos_y; while (next_x <= pos_x && next_y <= pos_y) { stack_counter++; random = GB_random(6); value = path_map->get_value(next_x, next_y); if (value == 0) ali_fatal_error("Unexpected value (1)", "ALI_PREALIGNER::mapper_random()"); switch (random) { case 0: if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); next_y--; } else { if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } else { stack->push(ALI_PREALIGNER_INS); next_x--; } } break; case 1: if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); next_y--; } else { if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); next_x--; } else { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } } break; case 2: if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } else { if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); next_y--; } else { stack->push(ALI_PREALIGNER_INS); next_x--; } } break; case 3: if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } else { if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); next_x--; } else { stack->push(ALI_PREALIGNER_DEL); next_y--; } } break; case 4: if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); next_x--; } else { if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); next_y--; } else { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } } break; case 5: if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); next_x--; } else { if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); next_x--; next_y--; } else { stack->push(ALI_PREALIGNER_DEL); next_y--; } } break; default: ali_fatal_error("Unexpected random value", "ALI_PREALIGNER::mapper_random()"); } } if (next_x <= pos_x) { mapper_post(stack, next_x + 1, 0); } else { if (next_y <= pos_y) { mapper_post(stack, 0, next_y + 1); } else { mapper_post(stack, 0, 0); } } if (stack_counter > 0) stack->pop(stack_counter); } void ALI_PREALIGNER::mapper(ALI_TSTACK < unsigned char >*stack, unsigned long pos_x, unsigned long pos_y) { // generate a stack of operations by taking every path unsigned char value; unsigned long stack_counter = 0; value = path_map->get_value(pos_x, pos_y); if (pos_x == 0 || pos_y == 0) { if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); if (pos_y == 0) mapper_post(stack, pos_x + 1, 0); else mapper(stack, pos_x, pos_y - 1); stack->pop(); } if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); if (pos_y > 0) { mapper_post(stack, 0, pos_y); } else { if (pos_x > 0) { mapper_post(stack, pos_x, 0); } else { mapper_post(stack, 0, 0); } } stack->pop(); } if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); if (pos_x == 0) mapper_post(stack, 0, pos_y + 1); else mapper(stack, pos_x - 1, pos_y); stack->pop(); } return; } // follow an unique path while (value == ALI_UP || value == ALI_DIAG || value == ALI_LEFT) { stack_counter++; switch (value) { case ALI_UP: stack->push(ALI_PREALIGNER_DEL); pos_y--; break; case ALI_DIAG: stack->push(ALI_PREALIGNER_SUB); pos_x--; pos_y--; break; case ALI_LEFT: stack->push(ALI_PREALIGNER_INS); pos_x--; break; } value = path_map->get_value(pos_x, pos_y); if (pos_x == 0 || pos_y == 0) { if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); if (pos_y == 0) mapper_post(stack, pos_x + 1, 0); else mapper(stack, pos_x, pos_y - 1); stack->pop(); } if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); if (pos_y > 0) { mapper_post(stack, 0, pos_y); } else { if (pos_x > 0) { mapper_post(stack, pos_x, 0); } else { mapper_post(stack, 0, 0); } } stack->pop(); } if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); if (pos_x == 0) mapper_post(stack, 0, pos_y + 1); else mapper(stack, pos_x - 1, pos_y); stack->pop(); } if (stack_counter > 0) stack->pop(stack_counter); return; } } if (value & ALI_UP) { stack->push(ALI_PREALIGNER_DEL); mapper(stack, pos_x, pos_y - 1); stack->pop(); } if (value & ALI_DIAG) { stack->push(ALI_PREALIGNER_SUB); mapper(stack, pos_x - 1, pos_y - 1); stack->pop(); } if (value & ALI_LEFT) { stack->push(ALI_PREALIGNER_INS); mapper(stack, pos_x - 1, pos_y); stack->pop(); } if (stack_counter > 0) stack->pop(stack_counter); } void ALI_PREALIGNER::make_map() { // make the result map from the path matrix unsigned long number_of_sol; ALI_TSTACK < unsigned char >*stack; stack = new ALI_TSTACK < unsigned char >(end_x - start_x + end_y - start_y + 3); ali_out_of_memory_if(!stack); number_of_sol = number_of_solutions(); printf("%lu solutions generated\n", number_of_sol); if (number_of_sol == 0 || number_of_sol > result_mask_counter) { ali_message("Starting random mapping"); do { mapper_random(stack, end_x - start_x, end_y - start_y); } while (result_mask_counter > 0); } else { ali_message("Starting systematic mapping"); mapper(stack, end_x - start_x, end_y - start_y); } delete stack; ali_message("Mapping finished"); } void ALI_PREALIGNER::generate_approximation(ALI_SUB_SOLUTION * work_sol) { // generate an approximation of a complete solution ALI_MAP *map; ALI_SEQUENCE *seq; char *ins_marker; float binding_costs; map = work_sol->make_one_map(); if (!map) ali_fatal_error("Can't make one map", "ALI_PREALIGNER::generate_approximation()"); seq = map->sequence_without_inserts(profile->sequence()); binding_costs = profile->w_binding(map->first_base(), seq); delete seq; ins_marker = map->insert_marker(); result_approx.insert(map, ins_marker, binding_costs); } void ALI_PREALIGNER::mapper_approximation(unsigned long area_no, ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>*map_lists, ALI_SUB_SOLUTION * work_sol) { // combine subsolutions for an approximation ALI_TLIST < ALI_MAP * >*map_list; ALI_MAP *map; // stop mapping at last area if (area_no > map_lists->size()) return; if (area_no == map_lists->size()) { generate_approximation(work_sol); return; } // map area number 'area_no' map_list = map_lists->get(area_no); if (map_list->is_empty()) ali_fatal_error("Found empty list", "ALI_PREALIGNER::mapper_approximation()"); // combine all possibilities map = map_list->first(); do { if (!work_sol->insert(map)) ali_fatal_error("Can't insert map", "ALI_PREALIGNER::mapper_approximation()"); mapper_approximation(area_no + 1, map_lists, work_sol); if (!work_sol->delete_map(map)) ali_fatal_error("Can't delete map", "ALI_PREALIGNER::mapper_approximation()"); if (map_list->has_next()) map = map_list->next(); else map = NULp; } while (map); } void ALI_PREALIGNER::make_approximation(ALI_PREALIGNER_CONTEXT * context) { // Make an approximation by aligning the undefined sections ALI_SUB_SOLUTION *work_solution; ALI_ALIGNER_CONTEXT aligner_context; ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>*map_lists; ALI_ALIGNER *aligner; unsigned long area_number; unsigned long start_seq, end_seq, start_ref, end_ref; ali_message("Align free areas"); work_solution = new ALI_SUB_SOLUTION(sub_solution); aligner_context.max_number_of_maps = context->max_number_of_maps_aligner; area_number = sub_solution->number_of_free_areas(); printf("number of areas = %ld (Maximal %ld solutions)\n", area_number, aligner_context.max_number_of_maps); map_lists = new ALI_TARRAY < ALI_TLIST < ALI_MAP * >*>(area_number); // generate Solutions for all free areas area_number = 0; while (sub_solution->free_area(&start_seq, &end_seq, &start_ref, &end_ref, area_number)) { printf("aligning area %ld (%ld,%ld) - (%ld,%ld)\n", area_number, start_seq, end_seq, start_ref, end_ref); aligner = new ALI_ALIGNER(&aligner_context, profile, start_seq, end_seq, start_ref, end_ref); map_lists->set(area_number, aligner->solutions()); printf("%d solutions generated\n", map_lists->get(area_number)->cardinality()); delete aligner; area_number++; } // combine and evaluate the solutions mapper_approximation(0, map_lists, work_solution); delete work_solution; ali_message("Free areas aligned"); } unsigned long ALI_PREALIGNER::number_of_solutions() { // approximate the number of solutions in the pathmap #define INFINIT 1000000 #define ADD(a, b) if (a>=INFINIT || b>=INFINIT) { a = INFINIT; } else { a += b; } unsigned long pos_x, pos_y, col_length; unsigned long number; unsigned char value; unsigned long *column1, *column2, *elem_akt_col, *elem_left_col; col_length = end_y - start_y + 1; column1 = (unsigned long *) CALLOC((unsigned int) col_length, sizeof(unsigned long)); column2 = (unsigned long *) CALLOC((unsigned int) col_length, sizeof(unsigned long)); ali_out_of_memory_if(!column1 || !column2); ali_message("Start: Checking number of solutions"); if (end_x - (start_x & 0x01)) { elem_akt_col = column1 + col_length - 1; elem_left_col = column2 + col_length - 1; } else { elem_akt_col = column2 + col_length - 1; elem_left_col = column1 + col_length - 1; } number = 0; *elem_akt_col = 1; for (pos_x = end_x - start_x; pos_x > 0;) { *(elem_left_col) = 0; for (pos_y = end_y - start_y; pos_y > 0; pos_y--) { *(elem_left_col - 1) = 0; value = path_map->get_value(pos_x, pos_y); if (value & ALI_UP) { ADD(*(elem_akt_col - 1), *elem_akt_col); } if (value & ALI_DIAG) { ADD(*(elem_left_col - 1), *elem_akt_col); } if (value & ALI_LEFT) { ADD(*(elem_left_col), *elem_akt_col); } elem_akt_col--; elem_left_col--; } value = path_map->get_value(pos_x, 0); if (value & ALI_UP) { ADD(number, *elem_akt_col); } if (value & ALI_DIAG) { ADD(number, *elem_akt_col); } if (value & ALI_LEFT) { ADD(*(elem_left_col), *elem_akt_col); } pos_x--; // toggle the columns if (pos_x & 0x01) { elem_akt_col = column1 + col_length - 1; elem_left_col = column2 + col_length - 1; } else { elem_akt_col = column2 + col_length - 1; elem_left_col = column1 + col_length - 1; } } for (pos_y = end_y - start_y; pos_y > 0; pos_y--) { value = path_map->get_value(0, pos_y); if (value & ALI_UP) { ADD(*(elem_akt_col - 1), *elem_akt_col); } if (value & ALI_DIAG) { ADD(number, *elem_akt_col); } if (value & ALI_LEFT) { ADD(number, *elem_akt_col); } elem_akt_col--; } ADD(number, *elem_akt_col); ali_message("End: Checking number of solutions"); free(column1); free(column2); return number; } ALI_PREALIGNER::ALI_PREALIGNER(ALI_PREALIGNER_CONTEXT * context, ALI_PROFILE * prof, unsigned long sx, unsigned long ex, unsigned long sy, unsigned long ey) { profile = prof; start_x = sx; end_x = ex; start_y = sy; end_y = ey; result_mask_counter = context->max_number_of_maps; ali_message("Prealigning"); path_map = new ALI_PATHMAP(end_x - start_x + 1, end_y - start_y + 1); calculate_matrix(); make_map(); result_mask.delete_expensive(context, profile); delete path_map; generate_solution(result_mask.map); make_approximation(context); ali_message("Prealigning finished"); }