diff --git a/ksw2_extd2.cpp b/ksw2_extd2.cpp index 499f77f..38a979a 100644 --- a/ksw2_extd2.cpp +++ b/ksw2_extd2.cpp @@ -264,11 +264,13 @@ int ksw_update_diag ( y2[t-1] = y2_new; H[t-1] = H_new; #ifdef DEBUG - // if (t >= t_st && t <= t_en+1) - // fprintf(align_debug_file, - // "t %d, s %d, u %d, v %d, x %d, y %d, x2 %d, y2 - // %d H %d\n", t - 1, sc[t-1], u_new, v_new, x_new, - // y_new, x2_new, y2_new, H_new); + if (r == 129){ + // if (t >= t_st && t <= t_en+1) + // fprintf(align_debug_file, + // "t %d, s %d, u %d, v %d, x %d, y %d, x2 %d, y2 %d H %d\n", + // t - 1, sc[t-1], u_new, v_new, x_new, + // y_new, x2_new, y2_new, H_new); + } #endif } @@ -333,7 +335,7 @@ int ksw_update_diag ( rmax[t] = r; } #ifdef DEBUG - // fprintf(align_debug_file, "%d|%d ", prev_H, H_new); + // fprintf(align_debug_file, "p %x ", p[t]); #endif } // update KZ matrix @@ -577,22 +579,29 @@ void ksw_extd2_cpp( // DEBUG: debug output #ifdef DEBUG fprintf(align_debug_file, "#%d (st=%d en=%d) ", r, st, en); + // for (int t = get_t(st0, r, n_col); t <= get_t(en0, r, n_col); ++t) { for (int t = t_st; t <= t_en; ++t) { // if (!align_debug_file) { // align_debug_file = fopen("debug/test_sample_debug.output", "w+"); // } - fprintf(align_debug_file, "%d ", H[t]); - // fprintf(align_debug_file, "(%d,%d,%d|%d,%d,%d,%d,%d,%d,%d,%d,%d,0x%x)\n", - // r, t, get_i(t, r, n_col), ((int8_t *)u)[t], ((int8_t *)v)[t], + // fprintf(align_debug_file, "%d ", H[t]); + fprintf(align_debug_file, "[%d %d v-1 %d s %d]%d ", v[t], u[t], v[t-1], sc[t], H[t]); + // fprintf(align_debug_file, + // "(%d,%d,%d|%d,%d,%d,%d,%d,%d,%d,%d,%d,0x%x)\n", + // r, t, get_i(t, r, n_col), ((int8_t *)u)[t], ((int8_t + // *)v)[t], // ((int8_t *)x)[t], ((int8_t *)y)[t], ((int8_t *)x2)[t], - // ((int8_t *)y2)[t], ((int32_t *)H)[t], ((int32_t *)Hmax)[t], - // ((int *)rmax)[t], p[r*min_n_col - t_st + t]); // for debugging + // ((int8_t *)y2)[t], ((int32_t *)H)[t], ((int32_t + // *)Hmax)[t], + // ((int *)rmax)[t], p[r*min_n_col - t_st + t]); // for + // debugging // if (!align_score_file) { - // align_score_file = fopen("debug/test_sample_score.output", "w+"); - // fprintf(align_score_file, "(r, t | u, v, x, y)\n"); + // align_score_file = fopen("debug/test_sample_score.output", + // "w+"); fprintf(align_score_file, "(r, t | u, v, x, y)\n"); // } - // fprintf(align_score_file, "(%d,%d|%d,%d,%d,%d)", - // r, t-t_st+st, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], + // fprintf(align_score_file, "(%d,%d|%d,%d,%d,%d)", + // r, t-t_st+st, ((int8_t*)u)[t], ((int8_t*)v)[t], + // ((int8_t*)x)[t], // ((int8_t*)y)[t]); // for debugging } // fprintf(align_score_file, "\n"); @@ -622,8 +631,8 @@ void ksw_extd2_cpp( // ez->score); for (int i = 0; i < qlen + tlen - 1; i++) { int len = i + 1 < qlen + tlen - i - 1 ? i + 1 : qlen + tlen - i - 1; - len = min_n_col < len ? min_n_col : len; - fprintf(align_score_file, "#%d ", i); + len = w < len ? w : len; + fprintf(align_score_file, "#%d (%d) ", i, len); for (int j = 0; j < len; j++) fprintf(align_score_file, "%x ", ((uint8_t *)p)[i * min_n_col + j]); fprintf(align_score_file, "\n"); diff --git a/ksw2_extd2_sse.c b/ksw2_extd2_sse.c index 9fe9e96..184e82e 100644 --- a/ksw2_extd2_sse.c +++ b/ksw2_extd2_sse.c @@ -158,13 +158,21 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin st0 = st, en0 = en; st = st / 16 * 16, en = (en + 16) / 16 * 16 - 1; // set boundary conditions - if (st > 0) { - if (st - 1 >= last_st && st - 1 <= last_en) { - x1 = x8[st - 1], x21 = x28[st - 1], v1 = v8[st - 1]; // (r-1,s-1) calculated in the last round - } else { + if (st0 > 0) { + if (st0 - 1 >= last_st && st0 - 1 <= last_en) { + if (r == 160) + printf("this is r == 160, last_st = %d, st0 = %d\n", + last_st, st0); + x1 = x8[st0 - 1], x21 = x28[st0 - 1], + v1 = v8[st0 - 1]; // (r-1,s-1) calculated in the last round + } else { x1 = -q - e, x21 = -q2 - e2; - v1 = -q - e; - } + // v1 = -q - e; + v1 = r == 0 ? -q - e + : r < long_thres ? -e + : r == long_thres ? long_diff + : -e2; + } } else { x1 = -q - e, x21 = -q2 - e2; v1 = r == 0? -q - e : r < long_thres? -e : r == long_thres? long_diff : -e2; @@ -393,12 +401,18 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin if (r == qlen + tlen - 2 && en0 == tlen - 1) ez->score = H0; } - last_st = st, last_en = en; + last_st = st0, last_en = en0; #ifdef DEBUG fprintf(align_debug_file, "#%d (st=%d en=%d) ", r, st0, en0); - // for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging + // for (t = st0; t <= en0; ++t) + // printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], + // ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for + // debugging for (t = st0; t <= en0; ++t) { - fprintf(align_debug_file, "%d ", H[t]); + // fprintf(align_debug_file, "%d ", H[t]); + fprintf(align_debug_file, "[%d %d v1 %d s %d p %x]%d ", v8[t], + u8[t], v1, ((int8_t *)s)[t], + ((uint8_t *)p)[r * n_col_ * 16 - st + t], H[t]); } fprintf(align_debug_file, "\n"); #endif // DEBUG @@ -409,10 +423,14 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin n_col = (n_col < w + 1 ? n_col : w + 1); for (int i = 0; i < qlen + tlen - 1; i++) { int len = i + 1 < qlen + tlen - i - 1 ? i + 1 : qlen + tlen - i - 1; - len = n_col < len ? n_col : len; - fprintf(align_score_file, "#%d ", i); - for (int j = 0; j < len; j++) { - fprintf(align_score_file, "%x ", ((uint8_t*)p)[i * n_col_*16 + j]); + len = w < len ? w : len; + int st = 0; + if (st < i - qlen + 1) st = i - qlen + 1; + if (st < (i - wr + 1) >> 1) st = (i - wr + 1) >> 1; // take the ceil + fprintf(align_score_file, "#%d (%d) ", i, len, off[i] - st); + for (int j = st - off[i]; j < len + st - off[i]; j++) { + fprintf(align_score_file, "%x ", + ((uint8_t *)p)[i * n_col_ * 16 + j]); } fprintf(align_score_file, "\n"); }