Separate kaln and kprobaln as I am preparing further changes. At

[samtools.git] / kaln.c

diff --git a/kaln.c b/kaln.c
index 65f6c497298e3ee5983307e4a011b18add99a54c..4cba98d3ae4bcc9d9e3ab9a30d537412dbd0f2ed 100644 (file)
--- a/kaln.c
+++ b/kaln.c
@@ -27,6 +27,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <stdint.h>
+#include <math.h>
 #include "kaln.h"
 
 #define FROM_M 0
@@ -72,8 +73,8 @@ int aln_sm_blast[] = {
        -2, -2, -2, -2, -2
 };
 
-ka_param_t ka_param_blast = {  5,  2,  2, aln_sm_blast, 5, 50 };
-ka_param_t ka_param_aa2aa = { 10,  2,  2, aln_sm_blosum62, 22, 50 };
+ka_param_t ka_param_blast = {  5,  2,   5, 2, aln_sm_blast, 5, 50 };
+ka_param_t ka_param_aa2aa = { 10,  2,  10, 2, aln_sm_blosum62, 22, 50 };
 
 static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
 {
@@ -141,13 +142,13 @@ static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
 }
 #define set_end_I(II, cur, p)                                                  \
 {                                                                                                              \
-       if (gap_end >= 0) {                                                                     \
-               if ((p)->M - gap_open > (p)->I) {                               \
+       if (gap_end_ext >= 0) {                                                         \
+               if ((p)->M - gap_end_open > (p)->I) {                   \
                        (cur)->It = FROM_M;                                                     \
-                       (II) = (p)->M - gap_open - gap_end;                     \
+                       (II) = (p)->M - gap_end_open - gap_end_ext;     \
                } else {                                                                                \
                        (cur)->It = FROM_I;                                                     \
-                       (II) = (p)->I - gap_end;                                        \
+                       (II) = (p)->I - gap_end_ext;                            \
                }                                                                                               \
        } else set_I(II, cur, p);                                                       \
 }
@@ -163,13 +164,13 @@ static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
 }
 #define set_end_D(DD, cur, p)                                                  \
 {                                                                                                              \
-       if (gap_end >= 0) {                                                                     \
-               if ((p)->M - gap_open > (p)->D) {                               \
+       if (gap_end_ext >= 0) {                                                         \
+               if ((p)->M - gap_end_open > (p)->D) {                   \
                        (cur)->Dt = FROM_M;                                                     \
-                       (DD) = (p)->M - gap_open - gap_end;                     \
+                       (DD) = (p)->M - gap_end_open - gap_end_ext;     \
                } else {                                                                                \
                        (cur)->Dt = FROM_D;                                                     \
-                       (DD) = (p)->D - gap_end;                                        \
+                       (DD) = (p)->D - gap_end_ext;                            \
                }                                                                                               \
        } else set_D(DD, cur, p);                                                       \
 }
@@ -182,17 +183,6 @@ typedef struct {
        int M, I, D;
 } dpscore_t;
 
-/* build score profile for accelerating alignment, in theory */
-static void aln_init_score_array(uint8_t *seq, int len, int row, int *score_matrix, int **s_array)
-{
-       int *tmp, *tmp2, i, k;
-       for (i = 0; i != row; ++i) {
-               tmp = score_matrix + i * row;
-               tmp2 = s_array[i];
-               for (k = 0; k != len; ++k)
-                       tmp2[k] = tmp[seq[k]];
-       }
-}
 /***************************
  * banded global alignment *
  ***************************/
@@ -206,17 +196,19 @@ uint32_t *ka_global_core(uint8_t *seq1, int len1, uint8_t *seq2, int len2, const
        uint8_t type, ctype;
        uint32_t *cigar = 0;
 
-       int gap_open, gap_ext, gap_end, b;
+       int gap_open, gap_ext, gap_end_open, gap_end_ext, b;
        int *score_matrix, N_MATRIX_ROW;
 
        /* initialize some align-related parameters. just for compatibility */
        gap_open = ap->gap_open;
        gap_ext = ap->gap_ext;
-       gap_end = ap->gap_end;
+       gap_end_open = ap->gap_end_open;
+       gap_end_ext = ap->gap_end_ext;
        b = ap->band_width;
        score_matrix = ap->matrix;
        N_MATRIX_ROW = ap->row;
-       
+
+       *n_cigar = 0;
        if (len1 == 0 || len2 == 0) return 0;
 
        /* calculate b1 and b2 */