Skip to content

Commit 3e3499a

Browse files
y3tsengy3tseng
y3tseng
authored and
y3tseng
committed
fix bug on tb exceed band
1 parent 5a93860 commit 3e3499a

File tree

1 file changed

+35
-17
lines changed

1 file changed

+35
-17
lines changed

‎src/TALCO-XDrop.cpp‎

Lines changed: 35 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@
2626
#include "TALCO-XDrop.hpp"
2727
#endif
2828

29+
using scoreType = float;
30+
2931
void Talco_xdrop::Align_freq (
3032
Params* params,
3133
const std::vector<std::vector<float>>& freqRef,
@@ -65,8 +67,10 @@ void Talco_xdrop::Align_freq (
6567
}
6668
for (int i= tile_aln.size()-1; i>=0; i--){
6769
if (i == tile_aln.size()-1 && tile>0) continue;
70+
// if (freqRef.size() == 1776 && freqQry.size() == 6124) std::cout << (tile_aln[i] & 0xFFFF);
6871
aln.push_back(tile_aln[i]);
6972
}
73+
// if (freqRef.size() == 1776 && freqQry.size() == 6124) std::cout << '\n';
7074
tile_aln.clear();
7175
tile++;
7276
}
@@ -202,14 +206,14 @@ void Talco_xdrop::Tile (
202206
) {
203207

204208
// Initialising variables
205-
int32_t inf = 2*param->xdrop + 1;
209+
scoreType inf = 2.0*param->xdrop + 1.0;
206210
int32_t marker = param->marker;
207211
bool converged = false; bool conv_logic = false;
208212
float refNum = num.first, qryNum = num.second;
209213
int32_t reference_length = reference.size() - reference_idx;
210214
int32_t query_length = query.size() - query_idx;
211215
int32_t fLen = std::min(param->fLen, std::min(reference_length, query_length));
212-
int32_t score = 0; int32_t max_score = 0; int32_t max_score_prime = -inf; // int32_t max_score_ref_idx = 0; int32_t max_score_query_idx = 0;
216+
scoreType score = 0; scoreType max_score = 0; scoreType max_score_prime = -inf; // int32_t max_score_ref_idx = 0; int32_t max_score_query_idx = 0;
213217
int32_t conv_score = 0; int32_t conv_value = 0; int32_t conv_ref_idx = 0; int32_t conv_query_idx = 0;
214218
int32_t tb_start_addr = 0; int32_t tb_start_ftr = 0; // int32_t max_score_start_addr = 0; int32_t max_score_start_ftr = 0;
215219
int8_t tb_state = 0;
@@ -255,7 +259,7 @@ void Talco_xdrop::Tile (
255259
// float gapExtendTail_qry = 0;
256260

257261
int32_t L[3], U[3];
258-
int32_t *S[3], *I[2], *D[2];
262+
scoreType *S[3], *I[2], *D[2];
259263
int32_t *CS[3], *CI[2], *CD[2];
260264
std::vector<int8_t> tb;
261265
std::vector<int32_t> ftr_length;
@@ -265,11 +269,11 @@ void Talco_xdrop::Tile (
265269
int32_t prev_conv_s = -1;
266270

267271
for (size_t sIndx=0; sIndx<3; sIndx++) { // Allocate memory for S, I, D, and CS array
268-
S[sIndx] = new int32_t [fLen];
272+
S[sIndx] = new scoreType [fLen];
269273
CS[sIndx] = new int32_t [fLen];
270274
if (sIndx < 2) {
271-
I[sIndx] = new int32_t [fLen];
272-
D[sIndx] = new int32_t [fLen];
275+
I[sIndx] = new scoreType [fLen];
276+
D[sIndx] = new scoreType [fLen];
273277
CI[sIndx] = new int32_t [fLen];
274278
CD[sIndx] = new int32_t [fLen];
275279
}
@@ -285,13 +289,13 @@ void Talco_xdrop::Tile (
285289
I[sIndx][sLenIndex] = -1;
286290
D[sIndx][sLenIndex] = -1;
287291
CI[sIndx][sLenIndex] = -1;
288-
CD[sIndx][sLenIndex] = -1;
292+
CD[sIndx][sLenIndex] = -2;
289293
}
290294
}
291295
}
292296

293297
if ((reference_length < 0) || (query_length < 0)) {
294-
std::cout << reference_length << " " << query_length << std::endl;
298+
// std::cout << reference_length << " " << query_length << std::endl;
295299
fprintf(stderr, "ERROR: Reference/Query index exceeded limit!\n");
296300
errorType = 3;
297301
aln.clear();
@@ -309,7 +313,7 @@ void Talco_xdrop::Tile (
309313
}
310314
for (int32_t k = 0; k < reference_length + query_length - 1; k++){
311315

312-
// if (reference.size() == 3304 && query.size() == 1948) printf("Tile: %d, k: %d, L: %d, U: %d, (%d, %d)\n", tile, k, L[k%3], U[k%3]+1, reference_length, query_length);
316+
// if (reference.size() == 1776 && query.size() == 6124) printf("Tile: %d, k: %d, L: %d, U: %d, (%d, %d)\n", tile, k, L[k%3], U[k%3]+1, reference_length, query_length);
313317
if (L[k%3] >= U[k%3]+1) { // No more cells to compute based on x-drop critieria
314318
last_tile = true;
315319
errorType = 1;
@@ -360,7 +364,7 @@ void Talco_xdrop::Tile (
360364
exit(1);
361365
}
362366

363-
int32_t match = -inf, insOp = -inf, delOp = -inf, insExt = -inf, delExt = -inf;
367+
scoreType match = -inf, insOp = -inf, delOp = -inf, insExt = -inf, delExt = -inf;
364368
int32_t offset = i-L[k%3];
365369
int32_t offsetDiag = L[k%3]-L[(k+1)%3]+offset-1; // L[0] - L[1] + 0 - 1
366370
int32_t offsetUp = L[k%3]-L[(k+2)%3]+offset;
@@ -370,13 +374,13 @@ void Talco_xdrop::Tile (
370374
if ((k==0) ||
371375
((offsetDiag >= 0) && (offsetDiag <= U[(k+1)%3]-L[(k+1)%3])) ||
372376
(tile == 0 && (i == 0 || j == 0 ))) {
373-
int32_t similarScore = 0;
377+
scoreType similarScore = 0;
374378
float numerator = 0;
375379
if (type == 0) {
376380
for (int l = 0; l < 6; ++l) {
377381
for (int m = 0; m < 6; ++m) {
378382
if (m == 5 && l == 5) numerator += 0;
379-
else if (m == 5 || l == 5) numerator += reference[reference_idx+j][l]*query[query_idx+i][m]*gapExtend;
383+
else if (m == 5 || l == 5) numerator += reference[reference_idx+j][l]*query[query_idx+i][m]*0;
380384
else numerator += reference[reference_idx+j][l]*query[query_idx+i][m]*param->scoreMatrix[m][l];
381385
}
382386
}
@@ -392,7 +396,8 @@ void Talco_xdrop::Tile (
392396
}
393397

394398

395-
similarScore = static_cast<int32_t>(std::nearbyint(numerator/denominator));
399+
// similarScore = static_cast<int32_t>(std::nearbyint(numerator/denominator));
400+
similarScore = numerator/denominator;
396401
if (tile == 0 && (i == 0 || j == 0 )) {
397402
if (i == 0 && j > 0) match = similarScore + gapOpenHead_ref + gapExtendHead_ref * (reference_idx+j - 1);
398403
else if (i > 0 && j == 0) match = similarScore + gapOpenHead_qry + gapExtendHead_qry * (query_idx+i - 1);
@@ -405,14 +410,18 @@ void Talco_xdrop::Tile (
405410
else match = S[(k+1)%3][offsetDiag] + similarScore;
406411
}
407412

408-
int32_t pos_gapOpen_ref = static_cast<int32_t>(std::nearbyint(gapOp[0][reference_idx+j]));
409-
int32_t pos_gapOpen_qry = static_cast<int32_t>(std::nearbyint(gapOp[1][query_idx+i]));
413+
// int32_t pos_gapOpen_ref = static_cast<int32_t>(std::nearbyint(gapOp[0][reference_idx+j]));
414+
// int32_t pos_gapOpen_qry = static_cast<int32_t>(std::nearbyint(gapOp[1][query_idx+i]));
415+
scoreType pos_gapOpen_ref = gapOp[0][reference_idx+j];
416+
scoreType pos_gapOpen_qry = gapOp[1][query_idx+i];
410417
// float gapRef = reference[reference_idx+j][reference[0].size()-1];
411418
// float gapQry = query[query_idx+i][query[0].size()-1];
412419

413420

414-
int32_t pos_gapExtend_ref = static_cast<int32_t>(std::nearbyint(gapEx[0][reference_idx+j]));
415-
int32_t pos_gapExtend_qry = static_cast<int32_t>(std::nearbyint(gapEx[1][query_idx+i]));
421+
// int32_t pos_gapExtend_ref = static_cast<int32_t>(std::nearbyint(gapEx[0][reference_idx+j]));
422+
// int32_t pos_gapExtend_qry = static_cast<int32_t>(std::nearbyint(gapEx[1][query_idx+i]));
423+
scoreType pos_gapExtend_ref = gapEx[0][reference_idx+j];
424+
scoreType pos_gapExtend_qry = gapEx[1][query_idx+i];
416425
// int32_t pos_gapExtend_ref = (reference[reference_idx+j][21] > 0) ? gapExtend / 2 : gapExtend;
417426
// int32_t pos_gapExtend_qry = (query[query_idx+i][21] > 0) ? gapExtend / 2 : gapExtend;
418427

@@ -648,6 +657,15 @@ void Talco_xdrop::Tile (
648657
tb_start_ftr = (tb_state == 3) ? ftr_length.size() - 2: ftr_length.size() - 1;
649658
}
650659
}
660+
661+
if (conv_query_idx == 65535) {
662+
conv_query_idx = 0;
663+
conv_ref_idx = param->marker;
664+
}
665+
else if (conv_query_idx == 65534) {
666+
conv_query_idx = param->marker;
667+
conv_ref_idx = 0;
668+
}
651669

652670
reference_idx += conv_ref_idx;
653671
query_idx += conv_query_idx;

0 commit comments

Comments
 (0)