Mercurial > repos > rliterman > csp2

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