1 #include <jni.h>
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <limits.h>
5 #include "align2_MultiStateAligner11tsJNI.h"
6
7 // C doesn't have min() or max() so we define our own
8 #define max(a,b) \
9 ({ __typeof__ (a) _a = (a); \
10 __typeof__ (b) _b = (b); \
11 _a > _b ? _a : _b; })
12
13 #define min(a,b) \
14 ({ __typeof__ (a) _a = (a); \
15 __typeof__ (b) _b = (b); \
16 _a < _b ? _a : _b; })
17
18 #define GAPLEN 128
19 #define GAPCOST max(1,GAPLEN/64)
20
21 #define GAPC '-'
22
23 #define MODE_MS 0
24 #define MODE_DEL 1
25 #define MODE_INS 2
26 #define MODE_SUB 3
27
28 #define AFFINE_ARRAYS 1
29
30 #define TIMEBITS 11
31 #define SCOREBITS (32-TIMEBITS)
32 #define MAX_TIME ((1<<TIMEBITS)-1)
33 #define MAX_SCORE (((1<<(SCOREBITS-1))-1)-2000)
34 #define MIN_SCORE (0-MAX_SCORE) //Keeps it 1 point above "BAD".
35
36 #define SCOREOFFSET TIMEBITS
37
38 #define TIMEMASK (~((-1)<<TIMEBITS))
39 #define SCOREMASK ((~((-1)<<SCOREBITS))<<SCOREOFFSET)
40
41 #define POINTS_NOREF 0 //default -110
42 #define POINTS_NOCALL 0
43 #define POINTS_MATCH 70 //default 50
44 #define POINTS_MATCH2 100 //Note: Changing to 90 substantially reduces false positives
45 #define POINTS_COMPATIBLE 50
46 #define POINTS_SUB -127 //default -133
47 #define POINTS_SUBR -147 //increased penalty if prior match streak was at most 1 (I have no idea why this improves things)
48 #define POINTS_SUB2 -51 //default -47
49 #define POINTS_SUB3 -25
50 #define POINTS_MATCHSUB -10
51 #define POINTS_INS -395 //default -251
52 #define POINTS_INS2 -39 //default -61
53 #define POINTS_INS3 -23 //default -20
54 #define POINTS_INS4 -8 //default -20
55 #define POINTS_DEL -472 //default -239
56 #define POINTS_DEL2 -33 //default -30
57 #define POINTS_DEL3 -9 //default -7
58 #define POINTS_DEL4 -1 //default -1
59 #define POINTS_DEL5 -1 //default -1
60 #define POINTS_DEL_REF_N -10 //default -10
61 #define POINTS_GAP (0-GAPCOST) //default -10
62
63 #define TIMESLIP 4
64 #define MASK5 (TIMESLIP-1)
65
66 #define BARRIER_I1 2
67 #define BARRIER_D1 3
68
69 #define LIMIT_FOR_COST_3 5
70 #define LIMIT_FOR_COST_4 20
71 #define LIMIT_FOR_COST_5 80
72
73 #define BAD (MIN_SCORE-1)
74
75 #define POINTSoff_NOREF (POINTS_NOREF<<SCOREOFFSET)
76 #define POINTSoff_NOCALL (POINTS_NOCALL<<SCOREOFFSET)
77 #define POINTSoff_MATCH (POINTS_MATCH<<SCOREOFFSET)
78 #define POINTSoff_MATCH2 (POINTS_MATCH2<<SCOREOFFSET)
79 #define POINTSoff_COMPATIBLE (POINTS_COMPATIBLE<<SCOREOFFSET)
80 #define POINTSoff_SUB (POINTS_SUB<<SCOREOFFSET)
81 #define POINTSoff_SUBR (POINTS_SUBR<<SCOREOFFSET)
82 #define POINTSoff_SUB2 (POINTS_SUB2<<SCOREOFFSET)
83 #define POINTSoff_SUB3 (POINTS_SUB3<<SCOREOFFSET)
84 #define POINTSoff_MATCHSUB (POINTS_MATCHSUB<<SCOREOFFSET)
85 #define POINTSoff_INS (POINTS_INS<<SCOREOFFSET)
86 #define POINTSoff_INS2 (POINTS_INS2<<SCOREOFFSET)
87 #define POINTSoff_INS3 (POINTS_INS3<<SCOREOFFSET)
88 #define POINTSoff_INS4 (POINTS_INS4<<SCOREOFFSET)
89 #define POINTSoff_DEL (POINTS_DEL<<SCOREOFFSET)
90 #define POINTSoff_DEL2 (POINTS_DEL2<<SCOREOFFSET)
91 #define POINTSoff_DEL3 (POINTS_DEL3<<SCOREOFFSET)
92 #define POINTSoff_DEL4 (POINTS_DEL4<<SCOREOFFSET)
93 #define POINTSoff_DEL5 (POINTS_DEL5<<SCOREOFFSET)
94 #define POINTSoff_GAP (POINTS_GAP<<SCOREOFFSET)
95 #define POINTSoff_DEL_REF_N (POINTS_DEL_REF_N<<SCOREOFFSET)
96 #define BADoff (BAD<<SCOREOFFSET)
97 #define MAXoff_SCORE (MAX_SCORE<<SCOREOFFSET)
98 #define MINoff_SCORE (MIN_SCORE<<SCOREOFFSET)
99
fillUnlimited(jbyte * read,jbyte * ref,jsize read_length,jsize ref_length,jint refStartLoc,jint refEndLoc,jint * result,jlong * iterationsUnlimited,jint * packed,jint * POINTSoff_SUB_ARRAY,jint * POINTSoff_INS_ARRAY,jint maxRows,jint maxColumns)100 void fillUnlimited(
101 jbyte * read,
102 jbyte * ref,
103 jsize read_length,
104 jsize ref_length,
105 jint refStartLoc,
106 jint refEndLoc,
107 jint * result,
108 jlong * iterationsUnlimited,
109 jint * packed,
110 jint * POINTSoff_SUB_ARRAY,
111 jint * POINTSoff_INS_ARRAY,
112 jint maxRows,
113 jint maxColumns
114 ) {
115
116 const jint rows=read_length;
117 const jint columns=refEndLoc-refStartLoc+1;
118
119 const jint maxGain=(read_length-1)*POINTSoff_MATCH2+POINTSoff_MATCH;
120 const jint subfloor=0-2*maxGain;
121 const jint BARRIER_I2=rows-BARRIER_I1, BARRIER_I2b=columns-1;
122 const jint BARRIER_D2=rows-BARRIER_D1;
123
124 const int sizeXY=(maxRows+1)*(maxColumns+1);
125 const int idxMS=MODE_MS*sizeXY;
126 const int idxDEL=MODE_DEL*sizeXY;
127 const int idxINS=MODE_INS*sizeXY;
128
129 //temporary, for finding a bug
130 if(rows>maxRows || columns>maxColumns){
131 printf("error\n"); exit(0);
132 }
133
134 for(int row=1; row<=rows; row++){
135 const int tmp1=(row-1)*(maxColumns+1);
136 const int tmp2=(row)*(maxColumns+1);
137 for(int col=1; col<=columns; col++){
138 (*iterationsUnlimited)++;
139
140 const jbyte call0=(row<2 ? (jbyte)'?' : read[row-2]);
141 const jbyte call1=read[row-1];
142 const jbyte ref0=(col<2 ? (jbyte)'!' : ref[refStartLoc+col-2]);
143 const jbyte ref1=ref[refStartLoc+col-1];
144
145 const jboolean match=(call1==ref1 && ref1!='N');
146 const jboolean prevMatch=(call0==ref0 && ref0!='N');
147
148 const jboolean gap=(ref1==GAPC);
149
150 if(gap){
151 packed[idxMS+tmp2+col]=subfloor;
152 }else{//Calculate match and sub scores
153
154 const jint scoreFromDiag=packed[idxMS+tmp1+col-1]&SCOREMASK;
155 const jint scoreFromDel=packed[idxDEL+tmp1+col-1]&SCOREMASK;
156 const jint scoreFromIns=packed[idxINS+tmp1+col-1]&SCOREMASK;
157 const jint streak=(packed[idxMS+tmp1+col-1]&TIMEMASK);
158
159 {//Calculate match/sub score
160
161 if(match){
162
163 const jint scoreMS=scoreFromDiag+(prevMatch ? POINTSoff_MATCH2 : POINTSoff_MATCH);
164 const jint scoreD=scoreFromDel+POINTSoff_MATCH;
165 const jint scoreI=scoreFromIns+POINTSoff_MATCH;
166
167 jint score;
168 jint time;
169 if(scoreMS>=scoreD && scoreMS>=scoreI){
170 score=scoreMS;
171 time=(prevMatch ? streak+1 : 1);
172 }else if(scoreD>=scoreI){
173 score=scoreD;
174 time=1;
175 }else{
176 score=scoreI;
177 time=1;
178 }
179
180 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
181 packed[idxMS+tmp2+col]=(score|time);
182
183 }else{
184
185 jint scoreMS;
186 if(ref1!='N' && call1!='N'){
187 scoreMS=scoreFromDiag+(prevMatch ? (streak<=1 ? POINTSoff_SUBR : POINTSoff_SUB) :
188 POINTSoff_SUB_ARRAY[streak+1]);
189 }else{
190 scoreMS=scoreFromDiag+POINTSoff_NOCALL;
191 }
192
193 const jint scoreD=scoreFromDel+POINTSoff_SUB; //+2 to move it as close as possible to the deletion / insertion
194 const jint scoreI=scoreFromIns+POINTSoff_SUB;
195
196 jint score;
197 jint time;
198 jbyte prevState;
199 if(scoreMS>=scoreD && scoreMS>=scoreI){
200 score=scoreMS;
201 time=(prevMatch ? 1 : streak+1);
202 prevState=MODE_MS;
203 }else if(scoreD>=scoreI){
204 score=scoreD;
205 time=1;
206 prevState=MODE_DEL;
207 }else{
208 score=scoreI;
209 time=1;
210 prevState=MODE_INS;
211 }
212
213 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
214 packed[idxMS+tmp2+col]=(score|time);
215 }
216 }
217 }
218
219 if(row<BARRIER_D1 || row>BARRIER_D2){
220 packed[idxDEL+tmp2+col]=subfloor;
221 }else{//Calculate DEL score
222
223 const jint streak=packed[idxDEL+tmp2+col-1]&TIMEMASK;
224
225 const jint scoreFromDiag=packed[idxMS+tmp2+col-1]&SCOREMASK;
226 const jint scoreFromDel=packed[idxDEL+tmp2+col-1]&SCOREMASK;
227
228 jint scoreMS=scoreFromDiag+POINTSoff_DEL;
229 jint scoreD=scoreFromDel+(streak==0 ? POINTSoff_DEL :
230 streak<LIMIT_FOR_COST_3 ? POINTSoff_DEL2 :
231 streak<LIMIT_FOR_COST_4 ? POINTSoff_DEL3 :
232 streak<LIMIT_FOR_COST_5 ? POINTSoff_DEL4 :
233 ((streak&MASK5)==0 ? POINTSoff_DEL5 : 0));
234
235 if(ref1=='N'){
236 scoreMS+=POINTSoff_DEL_REF_N;
237 scoreD+=POINTSoff_DEL_REF_N;
238 }else if(gap){
239 scoreMS+=POINTSoff_GAP;
240 scoreD+=POINTSoff_GAP;
241 }
242
243 jint score;
244 jint time;
245 jbyte prevState;
246 if(scoreMS>=scoreD){
247 score=scoreMS;
248 time=1;
249 prevState=MODE_MS;
250 }else{
251 score=scoreD;
252 time=streak+1;
253 prevState=MODE_DEL;
254 }
255
256 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
257 packed[idxDEL+tmp2+col]=(score|time);
258 }
259
260 //Calculate INS score
261 if(gap || (row<BARRIER_I1 && col>1) || (row>BARRIER_I2 && col<BARRIER_I2b)){
262 packed[idxINS+tmp2+col]=subfloor;
263 }else{//Calculate INS score
264
265 const jint streak=packed[idxINS+tmp1+col]&TIMEMASK;
266
267 const jint scoreFromDiag=packed[idxMS+tmp1+col]&SCOREMASK;
268 const jint scoreFromIns=packed[idxINS+tmp1+col]&SCOREMASK;
269
270 const jint scoreMS=scoreFromDiag+POINTSoff_INS;
271 const jint scoreI=scoreFromIns+POINTSoff_INS_ARRAY[streak+1];
272
273 jint score;
274 jint time;
275 jbyte prevState;
276 if(scoreMS>=scoreI){
277 score=scoreMS;
278 time=1;
279 prevState=MODE_MS;
280 }else{
281 score=scoreI;
282 time=streak+1;
283 prevState=MODE_INS;
284 }
285
286 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
287 packed[idxINS+tmp2+col]=(score|time);
288 }
289 }
290 }
291
292 jint maxCol=-1;
293 jint maxState=-1;
294 jint maxScore=INT_MIN;
295
296 const int tmp=rows*(maxColumns+1);
297 for(int state=0; state<3; state++){
298 for(int col=1; col<=columns; col++){
299 const int x=packed[(state)*sizeXY+tmp+col]&SCOREMASK;
300 if(x>maxScore){
301 maxScore=x;
302 maxCol=col;
303 maxState=state;
304 }
305 }
306 }
307 maxScore>>=SCOREOFFSET;
308
309 result[0]=rows;
310 result[1]=maxCol;
311 result[2]=maxState;
312 result[3]=maxScore;
313 return;
314 }
315
calcDelScoreOffset(jint len)316 jint calcDelScoreOffset(jint len){
317 if(len<=0){return 0;}
318 jint score=POINTSoff_DEL;
319
320 if(len>LIMIT_FOR_COST_5){
321 score+=((len-LIMIT_FOR_COST_5+MASK5)/TIMESLIP)*POINTSoff_DEL5;
322 len=LIMIT_FOR_COST_5;
323 }
324 if(len>LIMIT_FOR_COST_4){
325 score+=(len-LIMIT_FOR_COST_4)*POINTSoff_DEL4;
326 len=LIMIT_FOR_COST_4;
327 }
328 if(len>LIMIT_FOR_COST_3){
329 score+=(len-LIMIT_FOR_COST_3)*POINTSoff_DEL3;
330 len=LIMIT_FOR_COST_3;
331 }
332 if(len>1){
333 score+=(len-1)*POINTSoff_DEL2;
334 }
335 return score;
336 }
337
calcInsScoreOffset(jint len,jint * POINTSoff_INS_ARRAY_C)338 jint calcInsScoreOffset(jint len,
339 jint * POINTSoff_INS_ARRAY_C
340 ){
341 if(len<=0){return 0;}
342 if(AFFINE_ARRAYS){
343 return POINTSoff_INS_ARRAY_C[len];
344 }else{
345 jint score=POINTSoff_INS;
346 if(len>LIMIT_FOR_COST_4){
347 score+=(len-LIMIT_FOR_COST_4)*POINTSoff_INS4;
348 len=LIMIT_FOR_COST_4;
349 }
350 if(len>LIMIT_FOR_COST_3){
351 score+=(len-LIMIT_FOR_COST_3)*POINTSoff_INS3;
352 len=LIMIT_FOR_COST_3;
353 }
354 if(len>1){
355 score+=(len-1)*POINTSoff_INS2;
356 }
357 return score;
358 }
359 }
360
fillLimitedX(jbyte * read,jbyte * ref,jsize read_length,jsize ref_length,jint refStartLoc,jint refEndLoc,jint minScore,jint * result,jlong * iterationsLimited,jint * packed,jint * POINTSoff_SUB_ARRAY,jint * POINTSoff_INS_ARRAY,jint maxRows,jint maxColumns,jint bandwidth,jfloat bandwidthRatio,jint * vertLimit,jint * horizLimit,jbyte * baseToNumber,jint * POINTSoff_INS_ARRAY_C)361 void fillLimitedX(
362 jbyte * read,
363 jbyte * ref,
364 jsize read_length,
365 jsize ref_length,
366 jint refStartLoc,
367 jint refEndLoc,
368 jint minScore,
369 jint * result,
370 jlong * iterationsLimited,
371 jint * packed,
372 jint * POINTSoff_SUB_ARRAY,
373 jint * POINTSoff_INS_ARRAY,
374 jint maxRows,
375 jint maxColumns,
376 jint bandwidth,
377 jfloat bandwidthRatio,
378 jint * vertLimit,
379 jint * horizLimit,
380 jbyte * baseToNumber,
381 jint * POINTSoff_INS_ARRAY_C
382 ) {
383
384 const jint rows=read_length;
385 const jint columns=refEndLoc-refStartLoc+1;
386
387 const int sizeXY=(maxRows+1)*(maxColumns+1);
388 const int idxMS=MODE_MS*sizeXY;
389 const int idxDEL=MODE_DEL*sizeXY;
390 const int idxINS=MODE_INS*sizeXY;
391
392 const jint halfband=(bandwidth<1 && bandwidthRatio<=0) ? 0 :
393 max(min(bandwidth<1 ? 9999999 : bandwidth, bandwidthRatio<=0 ? 9999999 : 8+(jint)(rows*bandwidthRatio)), (columns-rows+8))/2;
394
395 const jint BARRIER_I2=rows-BARRIER_I1, BARRIER_I2b=columns-1;
396 const jint BARRIER_D2=rows-BARRIER_D1;
397
398 const int tmp=rows*(maxColumns+1);
399 for(int x=0; x<3; x++){
400 for(int i=1; i<columns+1; i++) {
401 packed[x*sizeXY+tmp+i]=BADoff;
402 }
403 }
404
405 jint minGoodCol=1;
406 jint maxGoodCol=columns;
407
408 const jint minScore_off=(minScore<<SCOREOFFSET);
409 const jint maxGain=(read_length-1)*POINTSoff_MATCH2+POINTSoff_MATCH;
410 const jint floor=minScore_off-maxGain;
411 const jint subfloor=floor-5*POINTSoff_MATCH2;
412
413 vertLimit[rows]=minScore_off;
414 jboolean prevDefined=0;
415 for(int i=rows-1; i>=0; i--){
416 jbyte c=read[i];
417 //if(AminoAcid.isFullyDefined(c)){
418 if(baseToNumber[c]>=0){
419 vertLimit[i]=max(vertLimit[i+1]-(prevDefined ? POINTSoff_MATCH2 : POINTSoff_MATCH), floor);
420 prevDefined=1;
421 }else{
422 vertLimit[i]=max(vertLimit[i+1]-POINTSoff_NOCALL, floor);
423 prevDefined=0;
424 }
425 }
426
427 horizLimit[columns]=minScore_off;
428 prevDefined=0;
429 for(int i=columns-1; i>=0; i--){
430 jbyte c=ref[refStartLoc+i];
431 if(baseToNumber[c]>=0){
432 horizLimit[i]=max(horizLimit[i+1]-(prevDefined ? POINTSoff_MATCH2 : POINTSoff_MATCH), floor);
433 prevDefined=1;
434 }else{
435 horizLimit[i]=max(horizLimit[i+1]-(prevDefined && c==GAPC ? POINTSoff_DEL : POINTSoff_NOREF), floor);
436 prevDefined=0;
437 }
438 }
439
440 for(int row=1; row<=rows; row++){
441 const jint colStart=(halfband<1 ? minGoodCol : max(minGoodCol, row-halfband));
442 const jint colStop=(halfband<1 ? maxGoodCol : min(maxGoodCol, row+halfband*2-1));
443
444 minGoodCol=-1;
445 maxGoodCol=-2;
446
447 const jint vlimit=vertLimit[row];
448
449 if(colStart<0 || colStop<colStart){break;}
450
451 if(colStart>1){
452 const int tmp3=(row)*(maxColumns+1)+(colStart-1);
453 packed[idxMS+tmp3]=subfloor;
454 packed[idxINS+tmp3]=subfloor;
455 packed[idxDEL+tmp3]=subfloor;
456 }
457
458 const int tmp1=(row-1)*(maxColumns+1);
459 const int tmp2=(row)*(maxColumns+1);
460 for(int col=colStart; col<=columns; col++){
461 const jbyte call0=(row<2 ? (jbyte)'?' : read[row-2]);
462 const jbyte call1=read[row-1];
463 const jbyte ref0=(col<2 ? (jbyte)'!' : ref[refStartLoc+col-2]);
464 const jbyte ref1=ref[refStartLoc+col-1];
465
466 const jboolean gap=(ref1==GAPC);
467
468 const jboolean match=(call1==ref1 && ref1!='N');
469 const jboolean prevMatch=(call0==ref0 && ref0!='N');
470
471 (*iterationsLimited)++;
472 const jint limit=max(vlimit, horizLimit[col]);
473 const jint limit3=max(floor, (match ? limit-POINTSoff_MATCH2 : limit-POINTSoff_SUB3));
474
475 const jint delNeeded=max(0, row-col-1);
476 const jint insNeeded=max(0, (rows-row)-(columns-col)-1);
477
478 const jint delPenalty=calcDelScoreOffset(delNeeded);
479 const jint insPenalty=calcInsScoreOffset(insNeeded,POINTSoff_INS_ARRAY_C);
480
481 const jint scoreFromDiag_MS=packed[idxMS+tmp1+col-1]&SCOREMASK;
482 const jint scoreFromDel_MS=packed[idxDEL+tmp1+col-1]&SCOREMASK;
483 const jint scoreFromIns_MS=packed[idxINS+tmp1+col-1]&SCOREMASK;
484
485 const jint scoreFromDiag_DEL=packed[idxMS+tmp2+col-1]&SCOREMASK;
486 const jint scoreFromDel_DEL=packed[idxDEL+tmp2+col-1]&SCOREMASK;
487
488 const jint scoreFromDiag_INS=packed[idxMS+tmp1+col]&SCOREMASK;
489 const jint scoreFromIns_INS=packed[idxINS+tmp1+col]&SCOREMASK;
490
491 if(gap || (scoreFromDiag_MS<=limit3 && scoreFromDel_MS<=limit3 && scoreFromIns_MS<=limit3)){
492 packed[idxMS+tmp2+col]=subfloor;
493 }else{//Calculate match and sub scores
494 const jint streak=(packed[idxMS+tmp1+col-1]&TIMEMASK);
495 {//Calculate match/sub score
496 jint score;
497 jint time;
498 jbyte prevState;
499
500 if(match){
501 const jint scoreMS=scoreFromDiag_MS+(prevMatch ? POINTSoff_MATCH2 : POINTSoff_MATCH);
502 const jint scoreD=scoreFromDel_MS+POINTSoff_MATCH;
503 const jint scoreI=scoreFromIns_MS+POINTSoff_MATCH;
504
505 if(scoreMS>=scoreD && scoreMS>=scoreI){
506 score=scoreMS;
507 time=(prevMatch ? streak+1 : 1);
508 prevState=MODE_MS;
509 }else if(scoreD>=scoreI){
510 score=scoreD;
511 time=1;
512 prevState=MODE_DEL;
513 }else{
514 score=scoreI;
515 time=1;
516 prevState=MODE_INS;
517 }
518 }else{
519 jint scoreMS;
520 if(ref1!='N' && call1!='N'){
521 scoreMS=scoreFromDiag_MS+(prevMatch ? (streak<=1 ? POINTSoff_SUBR : POINTSoff_SUB) :
522 POINTSoff_SUB_ARRAY[streak+1]);
523 }else{
524 scoreMS=scoreFromDiag_MS+POINTSoff_NOCALL;
525 }
526
527 const jint scoreD=scoreFromDel_MS+POINTSoff_SUB; //+2 to move it as close as possible to the deletion / insertion
528 const jint scoreI=scoreFromIns_MS+POINTSoff_SUB;
529
530 if(scoreMS>=scoreD && scoreMS>=scoreI){
531 score=scoreMS;
532 time=(prevMatch ? 1 : streak+1);
533 prevState=MODE_MS;
534 }else if(scoreD>=scoreI){
535 score=scoreD;
536 time=1;
537 prevState=MODE_DEL;
538 }else{
539 score=scoreI;
540 time=1;
541 prevState=MODE_INS;
542 }
543 }
544
545 jint limit2;
546 if(delNeeded>0){
547 limit2=limit-delPenalty;
548 }else if(insNeeded>0){
549 limit2=limit-insPenalty;
550 }else{
551 limit2=limit;
552 }
553
554 if(score>=limit2){
555 maxGoodCol=col;
556 if(minGoodCol<0){minGoodCol=col;}
557 }else{
558 score=subfloor;
559 }
560
561 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
562 packed[idxMS+tmp2+col]=(score|time);
563 }
564 }
565
566 if((scoreFromDiag_DEL<=limit && scoreFromDel_DEL<=limit) || row<BARRIER_D1 || row>BARRIER_D2){
567 packed[idxDEL+tmp2+col]=subfloor;
568 }else{//Calculate DEL score
569 const jint streak=packed[idxDEL+tmp2+col-1]&TIMEMASK;
570
571 jint scoreMS=scoreFromDiag_DEL+POINTSoff_DEL;
572 jint scoreD=scoreFromDel_DEL+(streak==0 ? POINTSoff_DEL :
573 streak<LIMIT_FOR_COST_3 ? POINTSoff_DEL2 :
574 streak<LIMIT_FOR_COST_4 ? POINTSoff_DEL3 :
575 streak<LIMIT_FOR_COST_5 ? POINTSoff_DEL4 :
576 ((streak&MASK5)==0 ? POINTSoff_DEL5 : 0));
577
578 if(ref1=='N'){
579 scoreMS+=POINTSoff_DEL_REF_N;
580 scoreD+=POINTSoff_DEL_REF_N;
581 }else if(gap){
582 scoreMS+=POINTSoff_GAP;
583 scoreD+=POINTSoff_GAP;
584 }
585
586 jint score;
587 jint time;
588 jbyte prevState;
589 if(scoreMS>=scoreD){
590 score=scoreMS;
591 time=1;
592 prevState=MODE_MS;
593 }else{
594 score=scoreD;
595 time=streak+1;
596 prevState=MODE_DEL;
597 }
598
599 jint limit2;
600 if(insNeeded>0){
601 limit2=limit-insPenalty;
602 }else if(delNeeded>0){
603 limit2=limit-calcDelScoreOffset(time+delNeeded)+calcDelScoreOffset(time);
604 }else{
605 limit2=limit;
606 }
607
608 if(score>=limit2){
609 maxGoodCol=col;
610 if(minGoodCol<0){minGoodCol=col;}
611 }else{
612 score=subfloor;
613 }
614
615 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
616 packed[idxDEL+tmp2+col]=(score|time);
617 }
618
619 if(gap || (scoreFromDiag_INS<=limit && scoreFromIns_INS<=limit) || (row<BARRIER_I1 && col>1) || (row>BARRIER_I2 && col<BARRIER_I2b)){
620 packed[idxINS+tmp2+col]=subfloor;
621 }else{//Calculate INS score
622 const jint streak=packed[idxINS+tmp1+col]&TIMEMASK;
623
624 const jint scoreMS=scoreFromDiag_INS+POINTSoff_INS;
625 const jint scoreI=scoreFromIns_INS+POINTSoff_INS_ARRAY[streak+1];
626
627 jint score;
628 jint time;
629 jbyte prevState;
630 if(scoreMS>=scoreI){
631 score=scoreMS;
632 time=1;
633 prevState=MODE_MS;
634 }else{
635 score=scoreI;
636 time=streak+1;
637 prevState=MODE_INS;
638 }
639
640 jint limit2;
641 if(delNeeded>0){
642 limit2=limit-delPenalty;
643 }else if(insNeeded>0){
644 limit2=limit-calcInsScoreOffset(time+insNeeded,POINTSoff_INS_ARRAY_C)+calcInsScoreOffset(time,POINTSoff_INS_ARRAY_C);
645 }else{
646 limit2=limit;
647 }
648
649 if(score>=limit2){
650 maxGoodCol=col;
651 if(minGoodCol<0){minGoodCol=col;}
652 }else{
653 score=subfloor;
654 }
655
656 if(time>MAX_TIME){time=MAX_TIME-MASK5;}
657 packed[idxINS+tmp2+col]=(score|time);
658 }
659
660 if(col>=colStop){
661 if(col>colStop && (maxGoodCol<col || halfband>0)){break;}
662 if(row>1){
663 const int tmp3=tmp1+(col+1);
664 packed[idxMS+tmp3]=subfloor;
665 packed[idxINS+tmp3]=subfloor;
666 packed[idxDEL+tmp3]=subfloor;
667 }
668 }
669 }
670 }
671
672 jint maxCol=-1;
673 jint maxState=-1;
674 jint maxScore=INT_MIN;
675
676 //const int tmp=rows*(maxColumns+1);
677 for(int state=0; state<3; state++){
678 for(int col=1; col<=columns; col++){
679 const int x=packed[(state)*sizeXY+tmp+col]&SCOREMASK;
680 if(x>maxScore){
681 maxScore=x;
682 maxCol=col;
683 maxState=state;
684 }
685 }
686 }
687
688 if(maxScore<minScore_off){
689 result[0]=rows;
690 result[1]=maxCol;
691 result[2]=maxState;
692 result[3]=maxScore;
693 result[4]=1;
694 return;
695 }
696
697 maxScore>>=SCOREOFFSET;
698 result[0]=rows;
699 result[1]=maxCol;
700 result[2]=maxState;
701 result[3]=maxScore;
702 result[4]=0;
703 return;
704 }
705
706
Java_align2_MultiStateAligner11tsJNI_fillUnlimitedJNI(JNIEnv * env,jobject obj,jbyteArray read,jbyteArray ref,jint refStartLoc,jint refEndLoc,jintArray result,jlongArray iterationsUnlimited,jintArray packed,jintArray POINTSoff_SUB_ARRAY,jintArray POINTSoff_INS_ARRAY,jint maxRows,jint maxColumns)707 JNIEXPORT void JNICALL Java_align2_MultiStateAligner11tsJNI_fillUnlimitedJNI(
708 JNIEnv *env,
709 jobject obj,
710 jbyteArray read,
711 jbyteArray ref,
712 jint refStartLoc,
713 jint refEndLoc,
714 jintArray result,
715 jlongArray iterationsUnlimited,
716 jintArray packed,
717 jintArray POINTSoff_SUB_ARRAY,
718 jintArray POINTSoff_INS_ARRAY,
719 jint maxRows,
720 jint maxColumns
721 ) {
722 // Get the size of the read and the reference arrays
723 jsize read_length=(*env)->GetArrayLength(env, read);
724 jsize ref_length=(*env)->GetArrayLength(env, ref);
725
726 // Copy arrays from Java
727 jint * jPOINTSoff_INS_ARRAY=(jint*)(*env)->GetPrimitiveArrayCritical(env, POINTSoff_INS_ARRAY, NULL);
728 jint * jPOINTSoff_SUB_ARRAY=(jint*)(*env)->GetPrimitiveArrayCritical(env, POINTSoff_SUB_ARRAY, NULL);
729 jint * jpacked=(jint*)(*env)->GetPrimitiveArrayCritical(env, packed, NULL);
730 jbyte * jread=(jbyte*)(*env)->GetPrimitiveArrayCritical(env, read, NULL);
731 jbyte * jref=(jbyte*)(*env)->GetPrimitiveArrayCritical(env, ref, NULL);
732 jint * jresult=(jint*)(*env)->GetPrimitiveArrayCritical(env, result, NULL);
733 jlong * jiterationsUnlimited=(jlong*)(*env)->GetPrimitiveArrayCritical(env, iterationsUnlimited, NULL);
734
735 // Using pointers for variables that need to be passed back to Java so they can be updated by the called functions
736 jlong * iterationsUnlimitedPointer=&jiterationsUnlimited[0];
737
738 // Call the fillUnlimited function in C; the 4 return values will be in jresult[]
739 fillUnlimited(jread,jref,read_length,ref_length,refStartLoc,refEndLoc,jresult,iterationsUnlimitedPointer,jpacked,jPOINTSoff_SUB_ARRAY,jPOINTSoff_INS_ARRAY,maxRows,maxColumns);
740
741 // Release Java arrays; 0 copies the array back to Java, JNI_ABORT does not copy the current array values to Java
742 (*env)->ReleasePrimitiveArrayCritical(env, result, jresult, 0);
743 (*env)->ReleasePrimitiveArrayCritical(env, iterationsUnlimited, jiterationsUnlimited, 0);
744 (*env)->ReleasePrimitiveArrayCritical(env, read, jread, JNI_ABORT);
745 (*env)->ReleasePrimitiveArrayCritical(env, ref, jref, JNI_ABORT);
746 (*env)->ReleasePrimitiveArrayCritical(env, packed, jpacked, 0);
747 (*env)->ReleasePrimitiveArrayCritical(env, POINTSoff_SUB_ARRAY, jPOINTSoff_SUB_ARRAY, JNI_ABORT);
748 (*env)->ReleasePrimitiveArrayCritical(env, POINTSoff_INS_ARRAY, jPOINTSoff_INS_ARRAY, JNI_ABORT);
749
750 return;
751 }
752
Java_align2_MultiStateAligner11tsJNI_fillLimitedXJNI(JNIEnv * env,jobject obj,jbyteArray read,jbyteArray ref,jint refStartLoc,jint refEndLoc,jint minScore,jintArray result,jlongArray iterationsLimited,jintArray packed,jintArray POINTSoff_SUB_ARRAY,jintArray POINTSoff_INS_ARRAY,jint maxRows,jint maxColumns,jint bandwidth,jfloat bandwidthRatio,jintArray vertLimit,jintArray horizLimit,jbyteArray baseToNumber,jintArray POINTSoff_INS_ARRAY_C)753 JNIEXPORT void JNICALL Java_align2_MultiStateAligner11tsJNI_fillLimitedXJNI(
754 JNIEnv *env,
755 jobject obj,
756 jbyteArray read,
757 jbyteArray ref,
758 jint refStartLoc,
759 jint refEndLoc,
760 jint minScore,
761 jintArray result,
762 jlongArray iterationsLimited,
763 jintArray packed,
764 jintArray POINTSoff_SUB_ARRAY,
765 jintArray POINTSoff_INS_ARRAY,
766 jint maxRows,
767 jint maxColumns,
768 jint bandwidth,
769 jfloat bandwidthRatio,
770 jintArray vertLimit,
771 jintArray horizLimit,
772 jbyteArray baseToNumber,
773 jintArray POINTSoff_INS_ARRAY_C
774 ) {
775 // Get the size of the read and the reference arrays
776 jsize read_length=(*env)->GetArrayLength(env, read);
777 jsize ref_length=(*env)->GetArrayLength(env, ref);
778
779 // Copy arrays from Java
780 jint * jPOINTSoff_INS_ARRAY=(jint*)(*env)->GetPrimitiveArrayCritical(env, POINTSoff_INS_ARRAY, NULL);
781 jint * jPOINTSoff_SUB_ARRAY=(jint*)(*env)->GetPrimitiveArrayCritical(env, POINTSoff_SUB_ARRAY, NULL);
782 jint * jpacked=(jint*)(*env)->GetPrimitiveArrayCritical(env, packed, NULL);
783 jbyte * jread=(jbyte*)(*env)->GetPrimitiveArrayCritical(env, read, NULL);
784 jbyte * jref=(jbyte*)(*env)->GetPrimitiveArrayCritical(env, ref, NULL);
785 jint * jresult=(jint*)(*env)->GetPrimitiveArrayCritical(env, result, NULL);
786 jlong * jiterationsLimited=(jlong*)(*env)->GetPrimitiveArrayCritical(env, iterationsLimited, NULL);
787 jint * jvertLimit=(jint*)(*env)->GetPrimitiveArrayCritical(env, vertLimit, NULL);
788 jint * jhorizLimit=(jint*)(*env)->GetPrimitiveArrayCritical(env, horizLimit, NULL);
789 jbyte * jbaseToNumber=(jbyte*)(*env)->GetPrimitiveArrayCritical(env, baseToNumber, NULL);
790 jint * jPOINTSoff_INS_ARRAY_C=(jint*)(*env)->GetPrimitiveArrayCritical(env, POINTSoff_INS_ARRAY_C, NULL);
791
792 // Using pointers for variables that need to be passed back to Java so they can be updated by the called functions
793 jlong * iterationsLimitedPointer=&jiterationsLimited[0];
794
795 // Call the fillLimitedX function in C; the 5 return values will be in jresult[]
796 fillLimitedX(jread,jref,read_length,ref_length,refStartLoc,refEndLoc,minScore,jresult,iterationsLimitedPointer,jpacked,jPOINTSoff_SUB_ARRAY,jPOINTSoff_INS_ARRAY,maxRows,maxColumns,bandwidth,bandwidthRatio,jvertLimit,jhorizLimit,jbaseToNumber,jPOINTSoff_INS_ARRAY_C);
797
798 // Release Java arrays; 0 copies the array back to Java, JNI_ABORT does not copy the current array values to Java
799 (*env)->ReleasePrimitiveArrayCritical(env, result, jresult, 0);
800 (*env)->ReleasePrimitiveArrayCritical(env, iterationsLimited, jiterationsLimited, 0);
801 (*env)->ReleasePrimitiveArrayCritical(env, read, jread, JNI_ABORT);
802 (*env)->ReleasePrimitiveArrayCritical(env, ref, jref, JNI_ABORT);
803 (*env)->ReleasePrimitiveArrayCritical(env, packed, jpacked, 0);
804 (*env)->ReleasePrimitiveArrayCritical(env, POINTSoff_SUB_ARRAY, jPOINTSoff_SUB_ARRAY, JNI_ABORT);
805 (*env)->ReleasePrimitiveArrayCritical(env, POINTSoff_INS_ARRAY, jPOINTSoff_INS_ARRAY, JNI_ABORT);
806 (*env)->ReleasePrimitiveArrayCritical(env, vertLimit, jvertLimit, 0);
807 (*env)->ReleasePrimitiveArrayCritical(env, horizLimit, jhorizLimit, 0);
808 (*env)->ReleasePrimitiveArrayCritical(env, baseToNumber, jbaseToNumber, JNI_ABORT);
809 (*env)->ReleasePrimitiveArrayCritical(env, POINTSoff_INS_ARRAY_C, jPOINTSoff_INS_ARRAY_C, JNI_ABORT);
810
811 return;
812 }
813
814