Mercurial > repos > rliterman > csp2

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