Mercurial > repos > rliterman > csp2

comparison CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/aligner/MultiStateAligner9PacBioAdapter2.java @ 68:5028fdace37b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
comparison
equal deleted inserted replaced
67:0e9998148a16 68:5028fdace37b
1 package aligner;
2
3 import java.util.Arrays;
4
5 import dna.AminoAcid;
6 import dna.ChromosomeArray;
7 import dna.Data;
8 import shared.Tools;
9 import stream.Read;
10 import stream.SiteScore;
11
12 /**
13 * Based on MSA9ts, with transform scores tweaked for PacBio. */
14 public final class MultiStateAligner9PacBioAdapter2 implements Aligner {
15
16
17 public MultiStateAligner9PacBioAdapter2(){
18 //// assert(maxColumns_>=200);
19 //// assert(maxRows_>=200);
20 // maxRows=maxRows_;
21 // maxColumns=maxColumns_;
22 // packed=new int[3][maxRows+1][maxColumns+1];
23 //
24 // vertLimit=new int[maxRows+1];
25 // horizLimit=new int[maxColumns+1];
26 // Arrays.fill(vertLimit, BADoff);
27 // Arrays.fill(horizLimit, BADoff);
28 //
29 //// for(int i=0; i<maxColumns+1; i++){
30 //// scores[0][i]=0-i;
31 //// }
32 //
33 // for(int matrix=0; matrix<packed.length; matrix++){
34 // for(int i=1; i<=maxRows; i++){
35 // for(int j=0; j<packed[matrix][i].length; j++){
36 // packed[matrix][i][j]|=BADoff;
37 // }
38 //// packed[matrix][i][0]|=MODE_INS;
39 // }
40 //// for(int i=0; i<maxRows+1; i++){
41 //// scores[matrix][i][0]=(i*POINTSoff_NOREF);
42 //// }
43 // for(int i=0; i<=maxRows; i++){
44 //
45 // int prevScore=(i<2 ? 0 : packed[matrix][i-1][0]);
46 // int score=(i<2 ? (i*POINTSoff_INS) :
47 // (i<LIMIT_FOR_COST_3 ? prevScore+POINTSoff_INS2 :
48 // (i<LIMIT_FOR_COST_4 ? prevScore+POINTSoff_INS3 : prevScore+POINTSoff_INS4)));
49 //
50 // packed[matrix][i][0]=score;
51 // }
52 //// for(int i=1; i<maxColumns+1; i++){
53 //// prevState[matrix][0][i]=MODE_DEL;
54 //// }
55 //// for(int i=0; i<=maxColumns; i++){
56 //// packed[matrix][0][i]|=MODE_MS;
57 //// }
58 // }
59 }
60
61 // private void prefillTopRow(){
62 // final int[] header=packed[0];
63 // final int qlen=rows;
64 // for(int i=0; i<=columns; i++){
65 // int x=columns-i+1;
66 // int qbases=qlen-x;
67 //
68 // //Minimal points to prefer a leftmost alignment
69 // header[i]=qbases<=0 ? 0 : ((-qbases)<<STARTOFFSET);
70 //
71 // //Forces consumption of query, but does not allow for insertions...
72 //// header[i]=qbases<=0 ? 0 : calcDelScoreOffset(qbases);
73 // }
74 // }
75 //
76 // private void prefillLeftColumnStartingAt(int i){
77 // packed[0][0]=MODE_MATCH;
78 // i=Tools.max(1, i);
79 // for(int score=MODE_INS+(POINTSoff_INS*i); i<=maxRows; i++){//Fill column 0 with insertions
80 // score+=POINTSoff_INS;
81 // packed[i][0]=score;
82 // }
83 // }
84
85 private void initialize(int rows_, int columns_){
86 rows=rows_;
87 columns=columns_;
88 if(rows<=maxRows && columns<=maxColumns){
89 // prefillTopRow();
90 return;
91 }
92
93 final int maxRows0=maxRows;
94 final int maxColumns0=maxColumns;
95 final int[][][] packed0=packed;
96
97
98 //Monotonic increase
99 maxRows=Tools.max(maxRows, rows+10);
100 maxColumns=Tools.max(maxColumns, columns+10);
101
102 //TODO: Use KillSwitch
103 packed=new int[3][maxRows+1][maxColumns+1];
104
105 vertLimit=new int[maxRows+1];
106 horizLimit=new int[maxColumns+1];
107 Arrays.fill(vertLimit, BADoff);
108 Arrays.fill(horizLimit, BADoff);
109
110 for(int matrix=0; matrix<packed.length; matrix++){
111 for(int i=1; i<=maxRows; i++){
112 for(int j=0; j<packed[matrix][i].length; j++){
113 packed[matrix][i][j]|=BADoff;
114 }
115 // packed[matrix][i][0]|=MODE_INS;
116 }
117 // for(int i=0; i<maxRows+1; i++){
118 // scores[matrix][i][0]=(i*POINTSoff_NOREF);
119 // }
120 for(int i=0; i<=maxRows; i++){
121
122 int prevScore=(i<2 ? 0 : packed[matrix][i-1][0]);
123 int score=(i<2 ? (i*POINTSoff_INS) :
124 (i<LIMIT_FOR_COST_3 ? prevScore+POINTSoff_INS2 :
125 (i<LIMIT_FOR_COST_4 ? prevScore+POINTSoff_INS3 : prevScore+POINTSoff_INS4)));
126
127 packed[matrix][i][0]=score;
128 }
129 // for(int i=1; i<maxColumns+1; i++){
130 // prevState[matrix][0][i]=MODE_DEL;
131 // }
132 // for(int i=0; i<=maxColumns; i++){
133 // packed[matrix][0][i]|=MODE_MS;
134 // }
135 }
136 }
137
138 /** return new int[] {rows, maxC, maxS, max};
139 * Will not fill areas that cannot match minScore */
140 @Override
141 public final int[] fillLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
142 return fillLimitedX(read, ref, refStartLoc, refEndLoc, minScore);
143 }
144
145
146 /** return new int[] {rows, maxC, maxS, max};
147 * Will not fill areas that cannot match minScore */
148 private final int[] fillLimitedX(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
149 initialize(read.length, refEndLoc-refStartLoc+1);
150
151 if(/*read.length<40 || */false || minScore<=0 || columns>read.length+Tools.min(100, read.length)){
152 // assert(false) : minScore;
153 return fillUnlimited(read, ref, refStartLoc, refEndLoc);
154 }
155
156 minScore-=100; //Increases quality trivially
157
158 assert(rows<=maxRows) : "Check that values are in-bounds before calling this function: "+rows+", "+maxRows+"\n"+
159 refStartLoc+", "+refEndLoc+", "+rows+", "+maxRows+", "+columns+", "+maxColumns+"\n"+new String(read)+"\n";
160 assert(columns<=maxColumns) : "Check that values are in-bounds before calling this function: "+columns+", "+maxColumns+"\n"+
161 refStartLoc+", "+refEndLoc+", "+rows+", "+maxRows+", "+columns+", "+maxColumns+"\n"+new String(read)+"\n";
162
163 assert(refStartLoc>=0) : "Check that values are in-bounds before calling this function: "+refStartLoc+"\n"+
164 refStartLoc+", "+refEndLoc+", "+rows+", "+maxRows+", "+columns+", "+maxColumns+"\n"+new String(read)+"\n";
165 assert(refEndLoc<ref.length) : "Check that values are in-bounds before calling this function: "+refEndLoc+", "+ref.length+"\n"+
166 refStartLoc+", "+refEndLoc+", "+rows+", "+maxRows+", "+columns+", "+maxColumns+"\n"+new String(read)+"\n";
167
168 // for(int x=0; x<packed.length; x++){
169 // for(int y=1; y<rows+1; y++){
170 // Arrays.fill(packed[x][y], 1, columns+1, BADoff);
171 // }
172 // }
173 for(int x=0; x<packed.length; x++){
174
175 // Arrays.fill(packed[x][1], 1, columns+1, BADoff);
176 Arrays.fill(packed[x][rows], 1, columns+1, BADoff);
177 // for(int y=1; y<rows+1; y++){
178 // Arrays.fill(packed[x][y], 1, columns+1, BADoff);
179 // }
180 }
181
182 int minGoodCol=1;
183 int maxGoodCol=columns;
184
185 final int minScore_off=(minScore<<SCOREOFFSET);
186 final int maxGain=(read.length-1)*POINTSoff_MATCH2+POINTSoff_MATCH;
187 final int floor=minScore_off-maxGain;
188 // final int subfloor=Tools.max(BADoff, floor-200*POINTSoff_MATCH2);
189 final int subfloor=floor-5*POINTSoff_MATCH2;
190 assert(subfloor>BADoff); //TODO: Actually, it needs to be substantially more.
191 assert(subfloor<minScore_off) : minScore_off+", "+floor+", "+BADoff+", "+subfloor;
192
193 if(verbose2){
194 System.out.println();
195 System.out.println("minScore="+minScore+"\t"+minScore_off);
196 System.out.println("maxGain="+(maxGain>>SCOREOFFSET)+"\t"+(maxGain));
197 System.out.println("floor="+(floor>>SCOREOFFSET)+"\t"+(floor));
198 System.out.println("subfloor="+(subfloor>>SCOREOFFSET)+"\t"+(subfloor));
199 System.out.println("BADoff="+(BADoff>>SCOREOFFSET)+"\t"+(BADoff));
200 System.out.println("maxGain="+(maxGain>>SCOREOFFSET)+"\t"+(maxGain));
201 System.out.println();
202 }
203
204 vertLimit[rows]=minScore_off;
205 for(int i=rows-1; i>=0; i--){
206 vertLimit[i]=Tools.max(vertLimit[i+1]-POINTSoff_MATCH2, floor);
207 }
208
209 horizLimit[columns]=minScore_off;
210 for(int i=columns-1; i>=0; i--){
211 horizLimit[i]=Tools.max(horizLimit[i+1]-POINTSoff_MATCH2, floor);
212 }
213
214 for(int row=1; row<=rows; row++){
215
216 int colStart=minGoodCol;
217 int colStop=maxGoodCol;
218
219 minGoodCol=-1;
220 maxGoodCol=-2;
221
222 final int vlimit=vertLimit[row];
223
224 if(verbose2){
225 System.out.println();
226 System.out.println("row="+row);
227 System.out.println("colStart="+colStart);
228 System.out.println("colStop="+colStop);
229 System.out.println("vlimit="+(vlimit>>SCOREOFFSET)+"\t"+(vlimit));
230 }
231
232 if(colStart<0 || colStop<colStart){break;}
233
234
235 if(colStart>1){
236 assert(row>0);
237 packed[MODE_MS][row][colStart-1]=subfloor;
238 packed[MODE_INS][row][colStart-1]=subfloor;
239 packed[MODE_DEL][row][colStart-1]=subfloor;
240 }
241
242
243 for(int col=colStart; col<=columns; col++){
244
245
246 if(verbose2){
247 System.out.println("\ncol "+col);
248 }
249
250 final byte call0=(row<2 ? (byte)'?' : read[row-2]);
251 final byte call1=read[row-1];
252 final byte ref0=(col<2 ? (byte)'!' : ref[refStartLoc+col-2]);
253 final byte ref1=ref[refStartLoc+col-1];
254
255 // final boolean match=(read[row-1]==ref[refStartLoc+col-1]);
256 // final boolean prevMatch=(row<2 || col<2 ? false : read[row-2]==ref[refStartLoc+col-2]);
257 final boolean match=(call1==ref1 && ref1!='N');
258 final boolean prevMatch=(call0==ref0 && ref0!='N');
259
260 // System.err.println("")
261
262 iterationsLimited++;
263 final int limit=Tools.max(vlimit, horizLimit[col]);
264 final int limit3=Tools.max(floor, (match ? limit-POINTSoff_MATCH2 : limit-POINTSoff_SUB3));
265
266 final int delNeeded=Tools.max(0, row-col-1);
267 final int insNeeded=Tools.max(0, (rows-row)-(columns-col)-1);
268
269 final int delPenalty=calcDelScoreOffset(delNeeded);
270 final int insPenalty=calcInsScoreOffset(insNeeded);
271
272
273 final int scoreFromDiag_MS=packed[MODE_MS][row-1][col-1]&SCOREMASK;
274 final int scoreFromDel_MS=packed[MODE_DEL][row-1][col-1]&SCOREMASK;
275 final int scoreFromIns_MS=packed[MODE_INS][row-1][col-1]&SCOREMASK;
276
277 final int scoreFromDiag_DEL=packed[MODE_MS][row][col-1]&SCOREMASK;
278 final int scoreFromDel_DEL=packed[MODE_DEL][row][col-1]&SCOREMASK;
279
280 final int scoreFromDiag_INS=packed[MODE_MS][row-1][col]&SCOREMASK;
281 final int scoreFromIns_INS=packed[MODE_INS][row-1][col]&SCOREMASK;
282
283 // if(scoreFromDiag_MS<limit3 && scoreFromDel_MS<limit3 && scoreFromIns_MS<limit3
284 // && scoreFromDiag_DEL<limit && scoreFromDel_DEL<limit && scoreFromDiag_INS<limit && scoreFromIns_INS<limit){
285 // iterationsLimited--; //A "fast" iteration
286 // }
287
288 if((scoreFromDiag_MS<=limit3 && scoreFromDel_MS<=limit3 && scoreFromIns_MS<=limit3)){
289 packed[MODE_MS][row][col]=subfloor;
290 }else{//Calculate match and sub scores
291 final int streak=(packed[MODE_MS][row-1][col-1]&TIMEMASK);
292
293 {//Calculate match/sub score
294
295 int score;
296 int time;
297 byte prevState;
298
299 if(match){
300
301 int scoreMS=scoreFromDiag_MS+(prevMatch ? POINTSoff_MATCH2 : POINTSoff_MATCH);
302 int scoreD=scoreFromDel_MS+POINTSoff_MATCH;
303 int scoreI=scoreFromIns_MS+POINTSoff_MATCH;
304
305 // byte prevState;
306 if(scoreMS>=scoreD && scoreMS>=scoreI){
307 score=scoreMS;
308 time=(prevMatch ? streak+1 : 1);
309 prevState=MODE_MS;
310 }else if(scoreD>=scoreI){
311 score=scoreD;
312 time=1;
313 prevState=MODE_DEL;
314 }else{
315 score=scoreI;
316 time=1;
317 prevState=MODE_INS;
318 }
319
320 // assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
321 // assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
322 // assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
323 //// packed[MODE_MS][row][col]=(score|prevState|time);
324 // packed[MODE_MS][row][col]=(score|time);
325 // assert((score&SCOREMASK)==score);
326 //// assert((prevState&MODEMASK)==prevState);
327 // assert((time&TIMEMASK)==time);
328
329 }else{
330
331 int scoreMS;
332 if(ref1!='N' && call1!='N'){
333 scoreMS=scoreFromDiag_MS+(prevMatch ? (streak<=1 ? POINTSoff_SUBR : POINTSoff_SUB) :
334 (streak==0 ? POINTSoff_SUB : streak<LIMIT_FOR_COST_3 ? POINTSoff_SUB2 : POINTSoff_SUB3));
335 }else{
336 scoreMS=scoreFromDiag_MS+POINTSoff_NOCALL;
337 }
338
339 int scoreD=scoreFromDel_MS+POINTSoff_SUB; //+2 to move it as close as possible to the deletion / insertion
340 int scoreI=scoreFromIns_MS+POINTSoff_SUB;
341
342 if(scoreMS>=scoreD && scoreMS>=scoreI){
343 score=scoreMS;
344 time=(prevMatch ? 1 : streak+1);
345 // time=(prevMatch ? (streak==1 ? 3 : 1) : streak+1);
346 prevState=MODE_MS;
347 }else if(scoreD>=scoreI){
348 score=scoreD;
349 time=1;
350 prevState=MODE_DEL;
351 }else{
352 score=scoreI;
353 time=1;
354 prevState=MODE_INS;
355 }
356
357 // assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
358 // assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
359 // assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
360 //// packed[MODE_MS][row][col]=(score|prevState|time);
361 // packed[MODE_MS][row][col]=(score|time);
362 // assert((score&SCOREMASK)==score);
363 //// assert((prevState&MODEMASK)==prevState);
364 // assert((time&TIMEMASK)==time);
365 }
366
367 final int limit2;
368 if(delNeeded>0){
369 limit2=limit-delPenalty;
370 }else if(insNeeded>0){
371 limit2=limit-insPenalty;
372 }else{
373 limit2=limit;
374 }
375 assert(limit2>=limit);
376
377 if(verbose2){System.err.println("MS: \tlimit2="+(limit2>>SCOREOFFSET)+"\t, score="+(score>>SCOREOFFSET));}
378
379 if(score>=limit2){
380 maxGoodCol=col;
381 if(minGoodCol<0){minGoodCol=col;}
382 }else{
383 score=subfloor;
384 }
385
386 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
387 assert(score>=MINoff_SCORE || score==BADoff) : "Score overflow - use MSA2 instead";
388 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
389 // packed[MODE_MS][row][col]=(score|prevState|time);
390 packed[MODE_MS][row][col]=(score|time);
391 assert((score&SCOREMASK)==score);
392 // assert((prevState&MODEMASK)==prevState);
393 assert((time&TIMEMASK)==time);
394 }
395 }
396
397 if((scoreFromDiag_DEL<=limit && scoreFromDel_DEL<=limit)){
398 // assert((scoreFromDiag_DEL<=limit && scoreFromDel_DEL<=limit)) : scoreFromDiag_DEL+", "+row;
399 packed[MODE_DEL][row][col]=subfloor;
400 }else{//Calculate DEL score
401
402 final int streak=packed[MODE_DEL][row][col-1]&TIMEMASK;
403
404 int scoreMS=scoreFromDiag_DEL+POINTSoff_DEL;
405 int scoreD=scoreFromDel_DEL+(streak==0 ? POINTSoff_DEL :
406 streak<LIMIT_FOR_COST_3 ? POINTSoff_DEL2 :
407 streak<LIMIT_FOR_COST_4 ? POINTSoff_DEL3 : POINTSoff_DEL4);
408 // int scoreI=scoreFromIns+POINTSoff_DEL;
409
410
411 if(ref1=='N'){
412 scoreMS+=POINTSoff_DEL_REF_N;
413 scoreD+=POINTSoff_DEL_REF_N;
414 }
415
416 //if(match){scoreMS=subfloor;}
417
418 int score;
419 int time;
420 byte prevState;
421 if(scoreMS>=scoreD){
422 score=scoreMS;
423 time=1;
424 prevState=MODE_MS;
425 }else{
426 score=scoreD;
427 time=streak+1;
428 prevState=MODE_DEL;
429 }
430
431 final int limit2;
432 if(insNeeded>0){
433 limit2=limit-insPenalty;
434 }else if(delNeeded>0){
435 limit2=limit-calcDelScoreOffset(time+delNeeded)+calcDelScoreOffset(time);
436 }else{
437 limit2=limit;
438 }
439 assert(limit2>=limit);
440 if(verbose2){System.err.println("DEL: \tlimit2="+(limit2>>SCOREOFFSET)+"\t, score="+(score>>SCOREOFFSET));}
441
442 if(score>=limit2){
443 maxGoodCol=col;
444 if(minGoodCol<0){minGoodCol=col;}
445 }else{
446 score=subfloor;
447 }
448
449 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
450 assert(score>=MINoff_SCORE || score==BADoff) : "Score overflow - use MSA2 instead";
451 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
452 // packed[MODE_DEL][row][col]=(score|prevState|time);
453 packed[MODE_DEL][row][col]=(score|time);
454 assert((score&SCOREMASK)==score);
455 // assert((prevState&MODEMASK)==prevState);
456 assert((time&TIMEMASK)==time);
457 }
458
459 if(scoreFromDiag_INS<=limit && scoreFromIns_INS<=limit){
460 packed[MODE_INS][row][col]=subfloor;
461 }else{//Calculate INS score
462
463 final int streak=packed[MODE_INS][row-1][col]&TIMEMASK;
464
465 int scoreMS=scoreFromDiag_INS+POINTSoff_INS;
466 // int scoreD=scoreFromDel+POINTSoff_INS;
467 int scoreI=scoreFromIns_INS+(streak==0 ? POINTSoff_INS :
468 streak<LIMIT_FOR_COST_3 ? POINTSoff_INS2 :
469 streak<LIMIT_FOR_COST_4 ? POINTSoff_INS3 : POINTSoff_INS4);
470
471
472 // System.err.println("("+row+","+col+")\t"+scoreFromDiag+"+"+POINTSoff_INS+"="+scoreM+", "+
473 // scoreFromSub+"+"+POINTSoff_INS+"="+scoreS+", "
474 // +scoreD+", "+scoreFromIns+"+"+
475 // (streak==0 ? POINTSoff_INS : streak<LIMIT_FOR_COST_3 ? POINTSoff_INS2 : POINTSoff_INS3)+"="+scoreI);
476
477 //if(match){scoreMS=subfloor;}
478
479 int score;
480 int time;
481 byte prevState;
482 if(scoreMS>=scoreI){
483 score=scoreMS;
484 time=1;
485 prevState=MODE_MS;
486 }else{
487 score=scoreI;
488 time=streak+1;
489 prevState=MODE_INS;
490 }
491
492 final int limit2;
493 if(delNeeded>0){
494 limit2=limit-delPenalty;
495 }else if(insNeeded>0){
496 limit2=limit-calcInsScoreOffset(time+insNeeded)+calcInsScoreOffset(time);
497 }else{
498 limit2=limit;
499 }
500 assert(limit2>=limit);
501
502 if(verbose2){System.err.println("INS: \tlimit2="+(limit2>>SCOREOFFSET)+"\t, score="+(score>>SCOREOFFSET));}
503 if(score>=limit2){
504 maxGoodCol=col;
505 if(minGoodCol<0){minGoodCol=col;}
506 }else{
507 score=subfloor;
508 }
509
510 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
511 assert(score>=MINoff_SCORE || score==BADoff) : "Score overflow - use MSA2 instead";
512 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
513 // packed[MODE_INS][row][col]=(score|prevState|time);
514 packed[MODE_INS][row][col]=(score|time);
515 assert((score&SCOREMASK)==score);
516 // assert((prevState&MODEMASK)==prevState);
517 assert((time&TIMEMASK)==time);
518 }
519
520
521 if(col>=colStop){
522 if(col>colStop && maxGoodCol<col){break;}
523 if(row>1){
524 packed[MODE_MS][row-1][col+1]=subfloor;
525 packed[MODE_INS][row-1][col+1]=subfloor;
526 packed[MODE_DEL][row-1][col+1]=subfloor;
527 }
528 }
529 }
530 }
531
532
533 int maxCol=-1;
534 int maxState=-1;
535 int maxScore=Integer.MIN_VALUE;
536
537 for(int state=0; state<packed.length; state++){
538 for(int col=1; col<=columns; col++){
539 int x=packed[state][rows][col]&SCOREMASK;
540 if(x>maxScore){
541 maxScore=x;
542 maxCol=col;
543 maxState=state;
544 }
545 }
546 }
547
548 assert(maxScore>=BADoff);
549 // if(maxScore==BADoff){
550 // return null;
551 // }
552 // if(maxScore<floor){
553 // return null;
554 // }
555 if(maxScore<minScore_off){
556 return null;
557 }
558
559 maxScore>>=SCOREOFFSET;
560
561 // System.err.println("Returning "+rows+", "+maxCol+", "+maxState+", "+maxScore);
562 return new int[] {rows, maxCol, maxState, maxScore};
563 }
564
565 @Override
566 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
567 return fillUnlimited(read, ref, refStartLoc, refEndLoc);
568 }
569
570 /** return new int[] {rows, maxCol, maxState, maxScore, maxStart}; */
571 @Override
572 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc){
573 initialize(read.length, refEndLoc-refStartLoc+1);
574
575 final int maxGain=(read.length-1)*POINTSoff_MATCH2+POINTSoff_MATCH;
576 final int subfloor=0-2*maxGain;
577 assert(subfloor>BADoff && subfloor*2>BADoff) : subfloor+", "+BADoff+", "+read.length; //TODO: Actually, it needs to be substantially more.
578
579 //temporary, for finding a bug
580 if(rows>maxRows || columns>maxColumns){
581 throw new RuntimeException("rows="+rows+", maxRows="+maxRows+", cols="+columns+", maxCols="+maxColumns+"\n"+new String(read)+"\n");
582 }
583
584 assert(rows<=maxRows) : "Check that values are in-bounds before calling this function: "+rows+", "+maxRows;
585 assert(columns<=maxColumns) : "Check that values are in-bounds before calling this function: "+columns+", "+maxColumns;
586
587 assert(refStartLoc>=0) : "Check that values are in-bounds before calling this function: "+refStartLoc;
588 assert(refEndLoc<ref.length) : "Check that values are in-bounds before calling this function: "+refEndLoc+", "+ref.length;
589
590 for(int row=1; row<=rows; row++){
591
592 // int minc=max(1, row-20);
593 // int maxc=min(columns, row+20);
594
595 for(int col=1; col<=columns; col++){
596 iterationsUnlimited++;
597
598 // final boolean match=(read[row-1]==ref[refStartLoc+col-1]);
599 // final boolean prevMatch=(row<2 || col<2 ? false : read[row-2]==ref[refStartLoc+col-2]);
600
601 final byte call0=(row<2 ? (byte)'?' : read[row-2]);
602 final byte call1=read[row-1];
603 final byte ref0=(col<2 ? (byte)'!' : ref[refStartLoc+col-2]);
604 final byte ref1=ref[refStartLoc+col-1];
605
606 final boolean match=(call1==ref1 && ref1!='N');
607 final boolean prevMatch=(call0==ref0 && ref0!='N');
608
609 {//Calculate match and sub scores
610
611 final int scoreFromDiag=packed[MODE_MS][row-1][col-1]&SCOREMASK;
612 final int scoreFromDel=packed[MODE_DEL][row-1][col-1]&SCOREMASK;
613 final int scoreFromIns=packed[MODE_INS][row-1][col-1]&SCOREMASK;
614 final int streak=(packed[MODE_MS][row-1][col-1]&TIMEMASK);
615
616 {//Calculate match/sub score
617
618 if(match){
619
620 int scoreMS=scoreFromDiag+(prevMatch ? POINTSoff_MATCH2 : POINTSoff_MATCH);
621 int scoreD=scoreFromDel+POINTSoff_MATCH;
622 int scoreI=scoreFromIns+POINTSoff_MATCH;
623
624 int score;
625 int time;
626 // byte prevState;
627 if(scoreMS>=scoreD && scoreMS>=scoreI){
628 score=scoreMS;
629 time=(prevMatch ? streak+1 : 1);
630 // prevState=MODE_MS;
631 }else if(scoreD>=scoreI){
632 score=scoreD;
633 time=1;
634 // prevState=MODE_DEL;
635 }else{
636 score=scoreI;
637 time=1;
638 // prevState=MODE_INS;
639 }
640
641 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
642 assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
643 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
644 // packed[MODE_MS][row][col]=(score|prevState|time);
645 packed[MODE_MS][row][col]=(score|time);
646 assert((score&SCOREMASK)==score);
647 // assert((prevState&MODEMASK)==prevState);
648 assert((time&TIMEMASK)==time);
649
650 }else{
651
652 int scoreMS;
653 if(ref1!='N' && call1!='N'){
654 scoreMS=scoreFromDiag+(prevMatch ? (streak<=1 ? POINTSoff_SUBR : POINTSoff_SUB) :
655 (streak==0 ? POINTSoff_SUB : streak<LIMIT_FOR_COST_3 ? POINTSoff_SUB2 : POINTSoff_SUB3));
656 }else{
657 scoreMS=scoreFromDiag+POINTSoff_NOCALL;
658 }
659
660 int scoreD=scoreFromDel+POINTSoff_SUB; //+2 to move it as close as possible to the deletion / insertion
661 int scoreI=scoreFromIns+POINTSoff_SUB;
662
663 int score;
664 int time;
665 byte prevState;
666 if(scoreMS>=scoreD && scoreMS>=scoreI){
667 score=scoreMS;
668 time=(prevMatch ? 1 : streak+1);
669 // time=(prevMatch ? (streak==1 ? 3 : 1) : streak+1);
670 prevState=MODE_MS;
671 }else if(scoreD>=scoreI){
672 score=scoreD;
673 time=1;
674 prevState=MODE_DEL;
675 }else{
676 score=scoreI;
677 time=1;
678 prevState=MODE_INS;
679 }
680
681 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
682 assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
683 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
684 // packed[MODE_MS][row][col]=(score|prevState|time);
685 packed[MODE_MS][row][col]=(score|time);
686 assert((score&SCOREMASK)==score);
687 // assert((prevState&MODEMASK)==prevState);
688 assert((time&TIMEMASK)==time);
689 }
690 }
691 }
692
693 {//Calculate DEL score
694
695 final int streak=packed[MODE_DEL][row][col-1]&TIMEMASK;
696
697 final int scoreFromDiag=packed[MODE_MS][row][col-1]&SCOREMASK;
698 final int scoreFromDel=packed[MODE_DEL][row][col-1]&SCOREMASK;
699
700 int scoreMS=scoreFromDiag+POINTSoff_DEL;
701 int scoreD=scoreFromDel+(streak==0 ? POINTSoff_DEL :
702 streak<LIMIT_FOR_COST_3 ? POINTSoff_DEL2 :
703 streak<LIMIT_FOR_COST_4 ? POINTSoff_DEL3 : POINTSoff_DEL4);
704 // int scoreI=scoreFromIns+POINTSoff_DEL;
705
706 if(ref1=='N'){
707 scoreMS+=POINTSoff_DEL_REF_N;
708 scoreD+=POINTSoff_DEL_REF_N;
709 }
710
711 //if(match){scoreMS=subfloor;}
712
713 int score;
714 int time;
715 byte prevState;
716 if(scoreMS>=scoreD){
717 score=scoreMS;
718 time=1;
719 prevState=MODE_MS;
720 }else{
721 score=scoreD;
722 time=streak+1;
723 prevState=MODE_DEL;
724 }
725
726 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
727 assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
728 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
729 // packed[MODE_DEL][row][col]=(score|prevState|time);
730 packed[MODE_DEL][row][col]=(score|time);
731 assert((score&SCOREMASK)==score);
732 // assert((prevState&MODEMASK)==prevState);
733 assert((time&TIMEMASK)==time);
734 }
735
736 {//Calculate INS score
737
738 final int streak=packed[MODE_INS][row-1][col]&TIMEMASK;
739
740 final int scoreFromDiag=packed[MODE_MS][row-1][col]&SCOREMASK;
741 final int scoreFromIns=packed[MODE_INS][row-1][col]&SCOREMASK;
742
743 int scoreMS=scoreFromDiag+POINTSoff_INS;
744 // int scoreD=scoreFromDel+POINTSoff_INS;
745 int scoreI=scoreFromIns+(streak==0 ? POINTSoff_INS :
746 streak<LIMIT_FOR_COST_3 ? POINTSoff_INS2 :
747 streak<LIMIT_FOR_COST_4 ? POINTSoff_INS3 : POINTSoff_INS4);
748
749 // System.err.println("("+row+","+col+")\t"+scoreFromDiag+"+"+POINTSoff_INS+"="+scoreM+", "+
750 // scoreFromSub+"+"+POINTSoff_INS+"="+scoreS+", "
751 // +scoreD+", "+scoreFromIns+"+"+
752 // (streak==0 ? POINTSoff_INS : streak<LIMIT_FOR_COST_3 ? POINTSoff_INS2 : POINTSoff_INS3)+"="+scoreI);
753
754 //if(match){scoreMS=subfloor;}
755
756 int score;
757 int time;
758 byte prevState;
759 if(scoreMS>=scoreI){
760 score=scoreMS;
761 time=1;
762 prevState=MODE_MS;
763 }else{
764 score=scoreI;
765 time=streak+1;
766 prevState=MODE_INS;
767 }
768
769 assert(time<=MAX_TIME);//if(time>MAX_TIME){time=MAX_TIME-3;}
770 assert(score>=MINoff_SCORE) : "Score overflow - use MSA2 instead";
771 assert(score<=MAXoff_SCORE) : "Score overflow - use MSA2 instead";
772 // packed[MODE_INS][row][col]=(score|prevState|time);
773 packed[MODE_INS][row][col]=(score|time);
774 assert((score&SCOREMASK)==score);
775 // assert((prevState&MODEMASK)==prevState);
776 assert((time&TIMEMASK)==time);
777 }
778 }
779 }
780
781
782 int maxCol=-1;
783 int maxState=-1;
784 int maxScore=Integer.MIN_VALUE;
785
786 for(int state=0; state<packed.length; state++){
787 for(int col=1; col<=columns; col++){
788 int x=packed[state][rows][col]&SCOREMASK;
789 if(x>maxScore){
790 maxScore=x;
791 maxCol=col;
792 maxState=state;
793 }
794 }
795 }
796 maxScore>>=SCOREOFFSET;
797
798 // System.err.println("Returning "+rows+", "+maxCol+", "+maxState+", "+maxScore);
799 return new int[] {rows, maxCol, maxState, maxScore};
800 }
801
802 @Override
803 public final byte[] traceback(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state){
804 // assert(false);
805 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
806 assert(row==rows);
807
808 byte[] out=new byte[row+col-1]; //TODO if an out of bound crash occurs, try removing the "-1".
809 int outPos=0;
810
811 if(state==MODE_INS){
812 //TODO ? Maybe not needed.
813 }
814
815 while(row>0 && col>0){
816
817 // byte prev0=(byte)(packed[state][row][col]&MODEMASK);
818
819 final int time=packed[state][row][col]&TIMEMASK;
820 final byte prev;
821
822 // System.err.println("state="+state+", prev="+prev+", row="+row+", col="+col+", score="+scores[state][row][col]);
823
824 if(state==MODE_MS){
825 if(time>1){prev=(byte)state;}
826 else{
827 final int scoreFromDiag=packed[MODE_MS][row-1][col-1]&SCOREMASK;
828 final int scoreFromDel=packed[MODE_DEL][row-1][col-1]&SCOREMASK;
829 final int scoreFromIns=packed[MODE_INS][row-1][col-1]&SCOREMASK;
830 if(scoreFromDiag>=scoreFromDel && scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
831 else if(scoreFromDel>=scoreFromIns){prev=MODE_DEL;}
832 else{prev=MODE_INS;}
833 }
834
835 byte c=read[row-1];
836 byte r=ref[refStartLoc+col-1];
837 if(c==r){
838 out[outPos]='m';
839 }else{
840 if(!AminoAcid.isFullyDefined(c)){
841 out[outPos]='N';
842 }else if(!AminoAcid.isFullyDefined(r)){
843 // out[outPos]='X';
844 out[outPos]='N';
845 }else{
846 out[outPos]='S';
847 }
848 }
849
850 row--;
851 col--;
852 }else if(state==MODE_DEL){
853 if(time>1){prev=(byte)state;}
854 else{
855 final int scoreFromDiag=packed[MODE_MS][row][col-1]&SCOREMASK;
856 final int scoreFromDel=packed[MODE_DEL][row][col-1]&SCOREMASK;
857 if(scoreFromDiag>=scoreFromDel){prev=MODE_MS;}
858 else{prev=MODE_DEL;}
859 }
860
861 byte r=ref[refStartLoc+col-1];
862 out[outPos]='D';
863
864 col--;
865 }else{
866 if(time>1){prev=(byte)state;}
867 else{
868 final int scoreFromDiag=packed[MODE_MS][row-1][col]&SCOREMASK;
869 final int scoreFromIns=packed[MODE_INS][row-1][col]&SCOREMASK;
870 if(scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
871 else{prev=MODE_INS;}
872 }
873
874 assert(state==MODE_INS) : state;
875 if(col==0){
876 out[outPos]='X';
877 }else if(col>=columns){
878 out[outPos]='Y';
879 }else{
880 out[outPos]='I';
881 }
882 row--;
883 }
884
885 // assert(prev==prev0);
886 state=prev;
887 outPos++;
888 }
889
890 assert(row==0 || col==0);
891 if(col!=row){
892 while(row>0){
893 out[outPos]='X';
894 outPos++;
895 row--;
896 col--;
897 }
898 if(col>0){
899 //do nothing
900 }
901 }
902
903
904 //Shrink and reverse the string
905 byte[] out2=new byte[outPos];
906 for(int i=0; i<outPos; i++){
907 out2[i]=out[outPos-i-1];
908 }
909 out=null;
910
911 return out2;
912 }
913
914 @Override
915 /** Generates identity;
916 * fills 'extra' with {match, sub, del, ins, N, clip} if present */
917 public float tracebackIdentity(byte[] query, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state, int[] extra){
918
919 // assert(false);
920 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
921 assert(row==rows);
922
923 // assert(state==(packed[row][col]&MODEMASK));
924 int match=0, sub=0, del=0, ins=0, noref=0, clip=0;
925
926 while(row>0 && col>0){
927
928 // byte prev0=(byte)(packed[state][row][col]&MODEMASK);
929
930 final int time=packed[state][row][col]&TIMEMASK;
931 final byte prev;
932
933 // System.err.println("state="+state+", prev="+prev+", row="+row+", col="+col+", score="+scores[state][row][col]);
934
935 if(state==MODE_MS){
936 if(time>1){prev=(byte)state;}
937 else{
938 final int scoreFromDiag=packed[MODE_MS][row-1][col-1]&SCOREMASK;
939 final int scoreFromDel=packed[MODE_DEL][row-1][col-1]&SCOREMASK;
940 final int scoreFromIns=packed[MODE_INS][row-1][col-1]&SCOREMASK;
941 if(scoreFromDiag>=scoreFromDel && scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
942 else if(scoreFromDel>=scoreFromIns){prev=MODE_DEL;}
943 else{prev=MODE_INS;}
944 }
945
946 byte q=query[row-1];
947 byte r=ref[refStartLoc+col-1];
948 boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
949
950 col--;
951 row--;
952 match+=(defined ? 1 : 0);
953 noref+=(defined ? 0 : 1);
954 }else if(state==MODE_DEL){
955 if(time>1){prev=(byte)state;}
956 else{
957 final int scoreFromDiag=packed[MODE_MS][row][col-1]&SCOREMASK;
958 final int scoreFromDel=packed[MODE_DEL][row][col-1]&SCOREMASK;
959 if(scoreFromDiag>=scoreFromDel){prev=MODE_MS;}
960 else{prev=MODE_DEL;}
961 }
962
963 col--;
964 del++;
965 }else{
966 if(time>1){prev=(byte)state;}
967 else{
968 final int scoreFromDiag=packed[MODE_MS][row-1][col]&SCOREMASK;
969 final int scoreFromIns=packed[MODE_INS][row-1][col]&SCOREMASK;
970 if(scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
971 else{prev=MODE_INS;}
972 }
973
974 assert(state==MODE_INS) : state;
975 row--;
976 boolean edge=(col<=1 || col>=columns);
977 ins+=(edge ? 0 : 1);
978 clip+=(edge ? 1 : 0);
979 }
980 state=prev;
981 }
982
983 assert(row==0 || col==0);
984 if(col!=row){//Not sure what this is doing
985 while(row>0){
986 clip++;
987 row--;
988 col--;
989 }
990 if(col>0){
991 //do nothing
992 }
993 }
994
995 if(extra!=null){
996 assert(extra.length==5);
997 extra[0]=match;
998 extra[1]=sub;
999 extra[2]=del;
1000 extra[3]=ins;
1001 extra[4]=noref;
1002 extra[5]=clip;
1003 }
1004
1005 int len=match+sub+ins+del;
1006 float id=match/Tools.max(1.0f, len);
1007 return id;
1008 }
1009
1010 @Override
1011 /** @return {score, bestRefStart, bestRefStop} */
1012 public final int[] score(final byte[] read, final byte[] ref, final int refStartLoc, final int refEndLoc,
1013 final int maxRow, final int maxCol, final int maxState){
1014
1015 int row=maxRow;
1016 int col=maxCol;
1017 int state=maxState;
1018
1019 assert(maxState>=0 && maxState<packed.length) :
1020 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
1021 assert(maxRow>=0 && maxRow<packed[0].length) :
1022 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
1023 assert(maxCol>=0 && maxCol<packed[0][0].length) :
1024 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
1025
1026 int score=packed[maxState][maxRow][maxCol]&SCOREMASK; //Or zero, if it is to be recalculated
1027
1028 if(row<rows){
1029 int difR=rows-row;
1030 int difC=columns-col;
1031
1032 while(difR>difC){
1033 score+=POINTSoff_NOREF;
1034 difR--;
1035 }
1036
1037 row+=difR;
1038 col+=difR;
1039
1040 }
1041
1042 assert(refStartLoc<=refEndLoc);
1043 assert(row==rows);
1044
1045
1046 final int bestRefStop=refStartLoc+col-1;
1047
1048 while(row>0 && col>0){
1049 // System.err.println("state="+state+", row="+row+", col="+col);
1050
1051
1052
1053 // byte prev0=(byte)(packed[state][row][col]&MODEMASK);
1054
1055 final int time=packed[state][row][col]&TIMEMASK;
1056 final byte prev;
1057
1058 if(state==MODE_MS){
1059 if(time>1){prev=(byte)state;}
1060 else{
1061 final int scoreFromDiag=packed[MODE_MS][row-1][col-1]&SCOREMASK;
1062 final int scoreFromDel=packed[MODE_DEL][row-1][col-1]&SCOREMASK;
1063 final int scoreFromIns=packed[MODE_INS][row-1][col-1]&SCOREMASK;
1064 if(scoreFromDiag>=scoreFromDel && scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
1065 else if(scoreFromDel>=scoreFromIns){prev=MODE_DEL;}
1066 else{prev=MODE_INS;}
1067 }
1068 row--;
1069 col--;
1070 }else if(state==MODE_DEL){
1071 if(time>1){prev=(byte)state;}
1072 else{
1073 final int scoreFromDiag=packed[MODE_MS][row][col-1]&SCOREMASK;
1074 final int scoreFromDel=packed[MODE_DEL][row][col-1]&SCOREMASK;
1075 if(scoreFromDiag>=scoreFromDel){prev=MODE_MS;}
1076 else{prev=MODE_DEL;}
1077 }
1078 col--;
1079 }else{
1080 assert(state==MODE_INS);
1081 if(time>1){prev=(byte)state;}
1082 else{
1083 final int scoreFromDiag=packed[MODE_MS][row-1][col]&SCOREMASK;
1084 final int scoreFromIns=packed[MODE_INS][row-1][col]&SCOREMASK;
1085 if(scoreFromDiag>=scoreFromIns){prev=MODE_MS;}
1086 else{prev=MODE_INS;}
1087 }
1088 row--;
1089 }
1090
1091 if(col<0){
1092 System.err.println(row);
1093 break; //prevents an out of bounds access
1094
1095 }
1096
1097 // assert(prev==prev0);
1098 state=prev;
1099
1100 // System.err.println("state2="+state+", row2="+row+", col2="+col+"\n");
1101 }
1102 // assert(false) : row+", "+col;
1103 if(row>col){
1104 col-=row;
1105 }
1106
1107 final int bestRefStart=refStartLoc+col;
1108
1109 score>>=SCOREOFFSET;
1110 int[] rvec;
1111 if(bestRefStart<refStartLoc || bestRefStop>refEndLoc){ //Suggest extra padding in cases of overflow
1112 int padLeft=Tools.max(0, refStartLoc-bestRefStart);
1113 int padRight=Tools.max(0, bestRefStop-refEndLoc);
1114 rvec=new int[] {score, bestRefStart, bestRefStop, padLeft, padRight};
1115 }else{
1116 rvec=new int[] {score, bestRefStart, bestRefStop};
1117 }
1118 return rvec;
1119 }
1120
1121
1122 /** Will not fill areas that cannot match minScore.
1123 * @return {score, bestRefStart, bestRefStop} */
1124 @Override
1125 public final int[] fillAndScoreLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
1126 initialize(read.length, refEndLoc-refStartLoc+1);
1127 int a=Tools.max(0, refStartLoc);
1128 int b=Tools.min(ref.length-1, refEndLoc);
1129 assert(b>=a);
1130
1131 int[] score;
1132
1133 if(b-a>=maxColumns){
1134 System.err.println("Warning: Max alignment columns exceeded; restricting range. "+(b-a+1)+" > "+maxColumns);//TODO: For testing; remove
1135 assert(false) : refStartLoc+", "+refEndLoc;
1136 b=Tools.min(ref.length-1, a+maxColumns-1);
1137 }
1138 int[] max=fillLimited(read, ref, a, b, minScore);
1139 score=(max==null ? null : score(read, ref, a, b, max[0], max[1], max[2]));
1140
1141 return score;
1142 }
1143
1144
1145
1146 public final static int scoreNoIndels(byte[] read, SiteScore ss){
1147 ChromosomeArray cha=Data.getChromosome(ss.chrom);
1148 return scoreNoIndels(read, cha.array, ss.start, ss);
1149 }
1150
1151 public final static int scoreNoIndels(byte[] read, final int chrom, final int refStart){
1152 ChromosomeArray cha=Data.getChromosome(chrom);
1153 return scoreNoIndels(read, cha.array, refStart, null);
1154 }
1155
1156 public final static int scoreNoIndels(byte[] read, SiteScore ss, byte[] baseScores){
1157 ChromosomeArray cha=Data.getChromosome(ss.chrom);
1158 return scoreNoIndels(read, cha.array, baseScores, ss.start, ss);
1159 }
1160
1161 public final static int scoreNoIndels(byte[] read, final int chrom, final int refStart, byte[] baseScores){
1162 ChromosomeArray cha=Data.getChromosome(chrom);
1163 return scoreNoIndels(read, cha.array, baseScores, refStart, null);
1164 }
1165
1166
1167
1168 public final static int scoreNoIndels(byte[] read, byte[] ref, final int refStart, final SiteScore ss){
1169
1170 int score=0;
1171 int mode=-1;
1172 int timeInMode=0;
1173
1174 //This block handles cases where the read runs outside the reference
1175 //Of course, padding the reference with 'N' would be better, but...
1176 int readStart=0;
1177 int readStop=read.length;
1178 final int refStop=refStart+read.length;
1179 boolean semiperfect=true;
1180 int norefs=0;
1181
1182 if(refStart<0){
1183 readStart=0-refStart;
1184 score+=POINTS_NOREF*readStart;
1185 norefs+=readStart;
1186 }
1187 if(refStop>ref.length){
1188 int dif=(refStop-ref.length);
1189 readStop-=dif;
1190 score+=POINTS_NOREF*dif;
1191 norefs+=dif;
1192 }
1193
1194 // if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
1195
1196 for(int i=readStart; i<readStop; i++){
1197 byte c=read[i];
1198 byte r=ref[refStart+i];
1199
1200 if(c==r && c!='N'){
1201 if(mode==MODE_MS){
1202 timeInMode++;
1203 score+=POINTS_MATCH2;
1204 }else{
1205 timeInMode=0;
1206 score+=POINTS_MATCH;
1207 }
1208 mode=MODE_MS;
1209 }else if(c<0 || c=='N'){
1210 score+=POINTS_NOCALL;
1211 semiperfect=false;
1212 }else if(r<0 || r=='N'){
1213 score+=POINTS_NOREF;
1214 norefs++;
1215 }else{
1216 if(mode==MODE_SUB){timeInMode++;}
1217 else{timeInMode=0;}
1218
1219 if(timeInMode==0){score+=POINTS_SUB;}
1220 else if(timeInMode<LIMIT_FOR_COST_3){score+=POINTS_SUB2;}
1221 else{score+=POINTS_SUB3;}
1222 mode=MODE_SUB;
1223 semiperfect=false;
1224 }
1225 }
1226
1227 if(semiperfect && ss!=null){ss.semiperfect=((ss.stop==ss.start+read.length-1) && (norefs<=read.length/2));}
1228
1229 return score;
1230 }
1231
1232
1233 public final static int scoreNoIndels(byte[] read, byte[] ref, byte[] baseScores, final int refStart, SiteScore ss){
1234
1235 int score=0;
1236 int mode=-1;
1237 int timeInMode=0;
1238 int norefs=0;
1239
1240 //This block handles cases where the read runs outside the reference
1241 //Of course, padding the reference with 'N' would be better, but...
1242 int readStart=0;
1243 int readStop=read.length;
1244 final int refStop=refStart+read.length;
1245 boolean semiperfect=true;
1246
1247 if(refStart<0){
1248 readStart=0-refStart;
1249 score+=POINTS_NOREF*readStart;
1250 norefs+=readStart;
1251 }
1252 if(refStop>ref.length){
1253 int dif=(refStop-ref.length);
1254 readStop-=dif;
1255 score+=POINTS_NOREF*dif;
1256 norefs+=dif;
1257 }
1258
1259 // if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
1260
1261 for(int i=readStart; i<readStop; i++){
1262 byte c=read[i];
1263 byte r=ref[refStart+i];
1264
1265 if(c==r && c!='N'){
1266 if(mode==MODE_MS){
1267 timeInMode++;
1268 score+=POINTS_MATCH2;
1269 }else{
1270 timeInMode=0;
1271 score+=POINTS_MATCH;
1272 }
1273 score+=baseScores[i];
1274 mode=MODE_MS;
1275 }else if(c<0 || c=='N'){
1276 score+=POINTS_NOCALL;
1277 semiperfect=false;
1278 }else if(r<0 || r=='N'){
1279 score+=POINTS_NOREF;
1280 norefs++;
1281 }else{
1282 if(mode==MODE_SUB){timeInMode++;}
1283 else{timeInMode=0;}
1284
1285 if(timeInMode==0){score+=POINTS_SUB;}
1286 else if(timeInMode<LIMIT_FOR_COST_3){score+=POINTS_SUB2;}
1287 else{score+=POINTS_SUB3;}
1288 mode=MODE_SUB;
1289 semiperfect=false;
1290 }
1291 }
1292
1293 if(semiperfect && ss!=null){ss.semiperfect=((ss.stop==ss.start+read.length-1) && (norefs<=read.length/2));}
1294 assert(Read.CHECKSITE(ss, read, -1));
1295
1296 return score;
1297 }
1298
1299
1300 public final static int scoreNoIndelsAndMakeMatchString(byte[] read, byte[] ref, byte[] baseScores, final int refStart, byte[][] matchReturn){
1301 int score=0;
1302 int mode=-1;
1303 int timeInMode=0;
1304
1305 assert(refStart<=ref.length) : refStart+", "+ref.length;
1306
1307 //This block handles cases where the read runs outside the reference
1308 //Of course, padding the reference with 'N' would be better, but...
1309 int readStart=0;
1310 int readStop=read.length;
1311 final int refStop=refStart+read.length;
1312 if(refStart<0){
1313 readStart=0-refStart;
1314 score+=POINTS_NOREF*readStart;
1315 }
1316 if(refStop>ref.length){
1317 int dif=(refStop-ref.length);
1318 System.err.println("dif="+dif+", ref.length="+ref.length+", refStop="+refStop);
1319 readStop-=dif;
1320 score+=POINTS_NOREF*dif;
1321 }
1322 assert(refStart+readStop<=ref.length) : "readStart="+readStart+", readStop="+readStop+
1323 ", refStart="+refStart+", refStop="+refStop+", ref.length="+ref.length+", read.length="+read.length;
1324
1325 assert(matchReturn!=null);
1326 assert(matchReturn.length==1);
1327 if(matchReturn[0]==null || matchReturn[0].length!=read.length){
1328 assert(matchReturn[0]==null || matchReturn[0].length<read.length) : matchReturn[0].length+"!="+read.length;
1329 matchReturn[0]=new byte[read.length];
1330 }
1331 final byte[] match=matchReturn[0];
1332
1333 // if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
1334
1335 for(int i=readStart; i<readStop; i++){
1336 byte c=read[i];
1337 byte r=ref[refStart+i];
1338
1339 assert(r!='.' && c!='.');
1340
1341 if(c==r && c!='N'){
1342 if(mode==MODE_MS){
1343 timeInMode++;
1344 score+=POINTS_MATCH2;
1345 }else{
1346 timeInMode=0;
1347 score+=POINTS_MATCH;
1348 }
1349 score+=baseScores[i];
1350 match[i]='m';
1351 mode=MODE_MS;
1352 }else if(c<0 || c=='N'){
1353 score+=POINTS_NOCALL;
1354 match[i]='N';
1355 }else if(r<0 || r=='N'){
1356 score+=POINTS_NOREF;
1357 // match[i]='m';
1358 match[i]='N';
1359 }else{
1360 match[i]='S';
1361 if(mode==MODE_SUB){timeInMode++;}
1362 else{timeInMode=0;}
1363
1364 if(timeInMode==0){score+=POINTS_SUB;}
1365 else if(timeInMode<LIMIT_FOR_COST_3){score+=POINTS_SUB2;}
1366 else{score+=POINTS_SUB3;}
1367 mode=MODE_SUB;
1368 }
1369 }
1370
1371 return score;
1372 }
1373
1374
1375 public final static int scoreNoIndelsAndMakeMatchString(byte[] read, byte[] ref, final int refStart, byte[][] matchReturn){
1376 int score=0;
1377 int mode=-1;
1378 int timeInMode=0;
1379
1380 assert(refStart<=ref.length) : refStart+", "+ref.length;
1381
1382 //This block handles cases where the read runs outside the reference
1383 //Of course, padding the reference with 'N' would be better, but...
1384 int readStart=0;
1385 int readStop=read.length;
1386 final int refStop=refStart+read.length;
1387 if(refStart<0){
1388 readStart=0-refStart;
1389 score+=POINTS_NOREF*readStart;
1390 }
1391 if(refStop>ref.length){
1392 int dif=(refStop-ref.length);
1393 System.err.println("dif="+dif+", ref.length="+ref.length+", refStop="+refStop);
1394 readStop-=dif;
1395 score+=POINTS_NOREF*dif;
1396 }
1397 assert(refStart+readStop<=ref.length) : "readStart="+readStart+", readStop="+readStop+
1398 ", refStart="+refStart+", refStop="+refStop+", ref.length="+ref.length+", read.length="+read.length;
1399
1400 assert(matchReturn!=null);
1401 assert(matchReturn.length==1);
1402 if(matchReturn[0]==null || matchReturn[0].length!=read.length){
1403 assert(matchReturn[0]==null || matchReturn[0].length<read.length) : matchReturn[0].length+"!="+read.length;
1404 matchReturn[0]=new byte[read.length];
1405 }
1406 final byte[] match=matchReturn[0];
1407
1408 // if(refStart<0 || refStart+read.length>ref.length){return -99999;} //No longer needed.
1409
1410 for(int i=readStart; i<readStop; i++){
1411 byte c=read[i];
1412 byte r=ref[refStart+i];
1413
1414 assert(r!='.' && c!='.');
1415
1416 if(c==r && c!='N'){
1417 if(mode==MODE_MS){
1418 timeInMode++;
1419 score+=POINTS_MATCH2;
1420 }else{
1421 timeInMode=0;
1422 score+=POINTS_MATCH;
1423 }
1424 match[i]='m';
1425 mode=MODE_MS;
1426 }else if(c<0 || c=='N'){
1427 score+=POINTS_NOCALL;
1428 match[i]='N';
1429 }else if(r<0 || r=='N'){
1430 score+=POINTS_NOREF;
1431 // match[i]='m';
1432 match[i]='N';
1433 }else{
1434 match[i]='S';
1435 if(mode==MODE_SUB){timeInMode++;}
1436 else{timeInMode=0;}
1437
1438 if(timeInMode==0){score+=POINTS_SUB;}
1439 else if(timeInMode<LIMIT_FOR_COST_3){score+=POINTS_SUB2;}
1440 else{score+=POINTS_SUB3;}
1441 mode=MODE_SUB;
1442 }
1443 }
1444
1445 return score;
1446 }
1447
1448 public static final int maxQuality(int numBases){
1449 return POINTS_MATCH+(numBases-1)*(POINTS_MATCH2);
1450 }
1451
1452 public static final int maxQuality(byte[] baseScores){
1453 return POINTS_MATCH+(baseScores.length-1)*(POINTS_MATCH2)+Tools.sumInt(baseScores);
1454 }
1455
1456 public static final int maxImperfectScore(int numBases){
1457 // int maxQ=maxQuality(numBases);
1458 //// maxImperfectSwScore=maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB);
1459 // int maxI=maxQ+POINTS_DEL;
1460 // maxI=Tools.max(maxI, maxQ+POINTS_INS-POINTS_MATCH2);
1461 // maxI=Tools.min(maxI, maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB));
1462
1463 int maxQ=maxQuality(numBases);
1464 int maxI=maxQ+Tools.min(POINTS_DEL, POINTS_INS-POINTS_MATCH2);
1465 assert(maxI<(maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB)));
1466 return maxI;
1467 }
1468
1469 public static final int maxImperfectScore(byte[] baseScores){
1470 // int maxQ=maxQuality(numBases);
1471 //// maxImperfectSwScore=maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB);
1472 // int maxI=maxQ+POINTS_DEL;
1473 // maxI=Tools.max(maxI, maxQ+POINTS_INS-POINTS_MATCH2);
1474 // maxI=Tools.min(maxI, maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB));
1475
1476 int maxQ=maxQuality(baseScores);
1477 int maxI=maxQ+Tools.min(POINTS_DEL, POINTS_INS-POINTS_MATCH2);
1478 assert(maxI<(maxQ-(POINTS_MATCH2+POINTS_MATCH2)+(POINTS_MATCH+POINTS_SUB)));
1479 return maxI;
1480 }
1481
1482 public static final String toString(int[] a){
1483
1484 int width=7;
1485
1486 StringBuilder sb=new StringBuilder((a.length+1)*width+2);
1487 for(int num : a){
1488 String s=" "+num;
1489 int spaces=width-s.length();
1490 assert(spaces>=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
1491 for(int i=0; i<spaces; i++){sb.append(' ');}
1492 sb.append(s);
1493 }
1494
1495 return sb.toString();
1496 }
1497
1498 public static final String toTimePacked(int[] a){
1499 int width=7;
1500
1501 StringBuilder sb=new StringBuilder((a.length+1)*width+2);
1502 for(int num_ : a){
1503 int num=num_&TIMEMASK;
1504 String s=" "+num;
1505 int spaces=width-s.length();
1506 assert(spaces>=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
1507 for(int i=0; i<spaces; i++){sb.append(' ');}
1508 sb.append(s);
1509 }
1510
1511 return sb.toString();
1512 }
1513
1514 public static final String toScorePacked(int[] a){
1515 int width=7;
1516
1517 String minString=" -";
1518 String maxString=" ";
1519 while(minString.length()<width){minString+='9';}
1520 while(maxString.length()<width){maxString+='9';}
1521
1522 StringBuilder sb=new StringBuilder((a.length+1)*width+2);
1523 for(int num_ : a){
1524 int num=num_>>SCOREOFFSET;
1525 String s=" "+num;
1526 if(s.length()>width){s=num>0 ? maxString : minString;}
1527 int spaces=width-s.length();
1528 assert(spaces>=0) : width+", "+s.length()+", "+s+", "+num+", "+spaces;
1529 for(int i=0; i<spaces; i++){sb.append(' ');}
1530 sb.append(s);
1531 }
1532
1533 return sb.toString();
1534 }
1535
1536 public static final String toString(byte[] a){
1537
1538 int width=6;
1539
1540 StringBuilder sb=new StringBuilder((a.length+1)*width+2);
1541 for(int num : a){
1542 String s=" "+num;
1543 int spaces=width-s.length();
1544 assert(spaces>=0);
1545 for(int i=0; i<spaces; i++){sb.append(' ');}
1546 sb.append(s);
1547 }
1548
1549 return sb.toString();
1550 }
1551
1552 public static final String toString(byte[] ref, int startLoc, int stopLoc){
1553 StringBuilder sb=new StringBuilder(stopLoc-startLoc+1);
1554 for(int i=startLoc; i<=stopLoc; i++){sb.append((char)ref[i]);}
1555 return sb.toString();
1556 }
1557
1558 // public int maxScoreByIdentity(int len, float identity){
1559 // assert(identity>=0 && identity<=1);
1560 // return (int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
1561 // }
1562
1563 @Override
1564 public int minScoreByIdentity(int len, float identity){
1565 assert(identity>=0 && identity<=1);
1566
1567 int a=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
1568 int b=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_INS));
1569 int c=(int)(len*(1*POINTS_MATCH+((1/(Tools.max(identity, 0.000001f)))-1)*POINTS_DEL));
1570 return Tools.min(a, b, c);
1571 }
1572
1573 public static int calcDelScore(int len){
1574 if(len<=0){return 0;}
1575 int score=POINTS_DEL;
1576
1577 if(len>LIMIT_FOR_COST_4){
1578 score+=(len-LIMIT_FOR_COST_4)*POINTS_DEL4;
1579 len=LIMIT_FOR_COST_4;
1580 }
1581 if(len>LIMIT_FOR_COST_3){
1582 score+=(len-LIMIT_FOR_COST_3)*POINTS_DEL3;
1583 len=LIMIT_FOR_COST_3;
1584 }
1585 if(len>1){
1586 score+=(len-1)*POINTS_DEL2;
1587 }
1588 return score;
1589 }
1590
1591 private static int calcDelScoreOffset(int len){
1592 if(len<=0){return 0;}
1593 int score=POINTSoff_DEL;
1594
1595 if(len>LIMIT_FOR_COST_4){
1596 score+=(len-LIMIT_FOR_COST_4)*POINTSoff_DEL4;
1597 len=LIMIT_FOR_COST_4;
1598 }
1599 if(len>LIMIT_FOR_COST_3){
1600 score+=(len-LIMIT_FOR_COST_3)*POINTSoff_DEL3;
1601 len=LIMIT_FOR_COST_3;
1602 }
1603 if(len>1){
1604 score+=(len-1)*POINTSoff_DEL2;
1605 }
1606 return score;
1607 }
1608
1609 public static int calcInsScore(int len){
1610 if(len<=0){return 0;}
1611 int score=POINTS_INS;
1612
1613 if(len>LIMIT_FOR_COST_4){
1614 score+=(len-LIMIT_FOR_COST_4)*POINTS_INS4;
1615 len=LIMIT_FOR_COST_4;
1616 }
1617 if(len>LIMIT_FOR_COST_3){
1618 score+=(len-LIMIT_FOR_COST_3)*POINTS_INS3;
1619 len=LIMIT_FOR_COST_3;
1620 }
1621 if(len>1){
1622 score+=(len-1)*POINTS_INS2;
1623 }
1624 return score;
1625 }
1626
1627 private static int calcInsScoreOffset(int len){
1628 if(len<=0){return 0;}
1629 int score=POINTSoff_INS;
1630 if(len>LIMIT_FOR_COST_4){
1631 score+=(len-LIMIT_FOR_COST_4)*POINTSoff_INS4;
1632 len=LIMIT_FOR_COST_4;
1633 }
1634 if(len>LIMIT_FOR_COST_3){
1635 score+=(len-LIMIT_FOR_COST_3)*POINTSoff_INS3;
1636 len=LIMIT_FOR_COST_3;
1637 }
1638 if(len>1){
1639 score+=(len-1)*POINTSoff_INS2;
1640 }
1641 return score;
1642 }
1643
1644 @Override
1645 public int rows(){return rows;}
1646 @Override
1647 public int columns(){return columns;}
1648
1649
1650 public int maxRows;
1651 public int maxColumns;
1652
1653 private int[][][] packed;
1654
1655 public int[] vertLimit;
1656 public int[] horizLimit;
1657
1658 CharSequence showVertLimit(){
1659 StringBuilder sb=new StringBuilder();
1660 for(int i=0; i<=rows; i++){sb.append(vertLimit[i]>>SCOREOFFSET).append(",");}
1661 return sb;
1662 }
1663 CharSequence showHorizLimit(){
1664 StringBuilder sb=new StringBuilder();
1665 for(int i=0; i<=columns; i++){sb.append(horizLimit[i]>>SCOREOFFSET).append(",");}
1666 return sb;
1667 }
1668
1669 // public static final int MODEBITS=2;
1670 public static final int TIMEBITS=12;
1671 public static final int SCOREBITS=32-TIMEBITS;
1672 public static final int MAX_TIME=((1<<TIMEBITS)-1);
1673 public static final int MAX_SCORE=((1<<(SCOREBITS-1))-1)-2000;
1674 public static final int MIN_SCORE=0-MAX_SCORE; //Keeps it 1 point above "BAD".
1675
1676 // public static final int MODEOFFSET=0; //Always zero.
1677 // public static final int TIMEOFFSET=0;
1678 public static final int SCOREOFFSET=TIMEBITS;
1679
1680 // public static final int MODEMASK=~((-1)<<MODEBITS);
1681 // public static final int TIMEMASK=(~((-1)<<TIMEBITS))<<TIMEOFFSET;
1682 public static final int TIMEMASK=~((-1)<<TIMEBITS);
1683 public static final int SCOREMASK=(~((-1)<<SCOREBITS))<<SCOREOFFSET;
1684
1685 private static final byte MODE_MS=0;
1686 private static final byte MODE_DEL=1;
1687 private static final byte MODE_INS=2;
1688 private static final byte MODE_SUB=3;
1689
1690 // public static final int POINTS_NOREF=-12;
1691 // public static final int POINTS_NOCALL=-12;
1692 // public static final int POINTS_MATCH=90;
1693 // public static final int POINTS_MATCH2=100;
1694 // public static final int POINTS_COMPATIBLE=50;
1695 // public static final int POINTS_SUB=-145;
1696 // public static final int POINTS_SUBR=-163; //increased penalty if prior match streak was at most 1
1697 // public static final int POINTS_SUB2=-59;
1698 // public static final int POINTS_SUB3=-35;
1699 // public static final int POINTS_MATCHSUB=-10;
1700 // public static final int POINTS_INS=-214;
1701 // public static final int POINTS_INS2=-55;
1702 // public static final int POINTS_INS3=-39;
1703 // public static final int POINTS_INS4=-17;
1704 // public static final int POINTS_DEL=-258;
1705 // public static final int POINTS_DEL2=-35;
1706 // public static final int POINTS_DEL3=-26;
1707 // public static final int POINTS_DEL4=-14;
1708 // public static final int POINTS_DEL_REF_N=-10;
1709
1710 public static final int POINTS_NOREF=-6;
1711 public static final int POINTS_NOCALL=-6;
1712 public static final int POINTS_MATCH=45;
1713 public static final int POINTS_MATCH2=50;
1714 public static final int POINTS_COMPATIBLE=25;
1715 public static final int POINTS_SUB=-72;
1716 public static final int POINTS_SUBR=-81; //increased penalty if prior match streak was at most 1
1717 public static final int POINTS_SUB2=-30;
1718 public static final int POINTS_SUB3=-18;
1719 public static final int POINTS_MATCHSUB=-5;
1720 public static final int POINTS_INS=-107;
1721 public static final int POINTS_INS2=-27;
1722 public static final int POINTS_INS3=-20;
1723 public static final int POINTS_INS4=-9;
1724 public static final int POINTS_DEL=-129;
1725 public static final int POINTS_DEL2=-18;
1726 public static final int POINTS_DEL3=-13;
1727 public static final int POINTS_DEL4=-7;
1728 public static final int POINTS_DEL_REF_N=-5;
1729
1730
1731 public static final int LIMIT_FOR_COST_3=5;
1732 public static final int LIMIT_FOR_COST_4=30;
1733
1734 public static final int BAD=MIN_SCORE-1;
1735
1736
1737 public static final int POINTSoff_NOREF=(POINTS_NOREF<<SCOREOFFSET);
1738 public static final int POINTSoff_NOCALL=(POINTS_NOCALL<<SCOREOFFSET);
1739 public static final int POINTSoff_MATCH=(POINTS_MATCH<<SCOREOFFSET);
1740 public static final int POINTSoff_MATCH2=(POINTS_MATCH2<<SCOREOFFSET);
1741 public static final int POINTSoff_COMPATIBLE=(POINTS_COMPATIBLE<<SCOREOFFSET);
1742 public static final int POINTSoff_SUB=(POINTS_SUB<<SCOREOFFSET);
1743 public static final int POINTSoff_SUBR=(POINTS_SUBR<<SCOREOFFSET);
1744 public static final int POINTSoff_SUB2=(POINTS_SUB2<<SCOREOFFSET);
1745 public static final int POINTSoff_SUB3=(POINTS_SUB3<<SCOREOFFSET);
1746 public static final int POINTSoff_MATCHSUB=(POINTS_MATCHSUB<<SCOREOFFSET);
1747 public static final int POINTSoff_INS=(POINTS_INS<<SCOREOFFSET);
1748 public static final int POINTSoff_INS2=(POINTS_INS2<<SCOREOFFSET);
1749 public static final int POINTSoff_INS3=(POINTS_INS3<<SCOREOFFSET);
1750 public static final int POINTSoff_INS4=(POINTS_INS4<<SCOREOFFSET);
1751 public static final int POINTSoff_DEL=(POINTS_DEL<<SCOREOFFSET);
1752 public static final int POINTSoff_DEL2=(POINTS_DEL2<<SCOREOFFSET);
1753 public static final int POINTSoff_DEL3=(POINTS_DEL3<<SCOREOFFSET);
1754 public static final int POINTSoff_DEL4=(POINTS_DEL4<<SCOREOFFSET);
1755 public static final int POINTSoff_DEL_REF_N=(POINTS_DEL_REF_N<<SCOREOFFSET);
1756 public static final int BADoff=(BAD<<SCOREOFFSET);
1757 public static final int MAXoff_SCORE=MAX_SCORE<<SCOREOFFSET;
1758 public static final int MINoff_SCORE=MIN_SCORE<<SCOREOFFSET;
1759
1760 private int rows;
1761 private int columns;
1762
1763 public long iterationsLimited=0;
1764 public long iterationsUnlimited=0;
1765
1766 public boolean verbose=false;
1767 public boolean verbose2=false;
1768
1769 }