Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/aligner/SingleStateAlignerFlatFloat.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 dna.AminoAcid;
jpayne@68 4 import shared.KillSwitch;
jpayne@68 5 import shared.Tools;
jpayne@68 6
jpayne@68 7 /**
jpayne@68 8 * Based on SSAFlat, but with previous state pointers removed. */
jpayne@68 9 public final class SingleStateAlignerFlatFloat implements Aligner {
jpayne@68 10
jpayne@68 11
jpayne@68 12 public SingleStateAlignerFlatFloat(){}
jpayne@68 13
jpayne@68 14 private void prefillTopRow(){
jpayne@68 15 final float[] header=packed[0];
jpayne@68 16 final int qlen=rows;
jpayne@68 17 for(int i=0; i<=columns; i++){
jpayne@68 18 int x=columns-i+1;
jpayne@68 19 int qbases=qlen-x;
jpayne@68 20
jpayne@68 21 //Minimal points to prefer a leftmost alignment
jpayne@68 22 header[i]=qbases<=0 ? 0 : -qbases;
jpayne@68 23
jpayne@68 24 //Forces consumption of query, but does not allow for insertions...
jpayne@68 25 // header[i]=qbases<=0 ? 0 : calcDelScoreOffset(qbases);
jpayne@68 26 }
jpayne@68 27 }
jpayne@68 28
jpayne@68 29 private void prefillLeftColumnStartingAt(int i){
jpayne@68 30 packed[0][0]=MODE_MATCH;
jpayne@68 31 i=Tools.max(1, i);
jpayne@68 32 for(float score=MODE_INS+(POINTS_INS*i); i<=maxRows; i++){//Fill column 0 with insertions
jpayne@68 33 score+=POINTS_INS;
jpayne@68 34 packed[i][0]=score;
jpayne@68 35 }
jpayne@68 36 }
jpayne@68 37
jpayne@68 38 private void initialize(int rows_, int columns_){
jpayne@68 39 rows=rows_;
jpayne@68 40 columns=columns_;
jpayne@68 41 if(rows<=maxRows && columns<=maxColumns){
jpayne@68 42 prefillTopRow();
jpayne@68 43 // prefillLeftColumn();
jpayne@68 44 return;
jpayne@68 45 }
jpayne@68 46
jpayne@68 47 final int maxRows0=maxRows;
jpayne@68 48 final int maxColumns0=maxColumns;
jpayne@68 49 final float[][] packed0=packed;
jpayne@68 50
jpayne@68 51 //Monotonic increase
jpayne@68 52 maxRows=Tools.max(maxRows, rows+10);
jpayne@68 53 maxColumns=Tools.max(maxColumns, columns+10);
jpayne@68 54
jpayne@68 55 if(packed==null || maxColumns>maxColumns0){//Make a new matrix
jpayne@68 56 packed=KillSwitch.allocFloat2D(maxRows+1, maxColumns+1);
jpayne@68 57 prefillLeftColumnStartingAt(1);
jpayne@68 58 }else{//Copy old rows
jpayne@68 59 assert(maxRows0>0 && maxColumns0>0);
jpayne@68 60 assert(maxRows>maxRows0 && maxColumns<=maxColumns0) : "rows="+rows+",maxRows="+maxRows+
jpayne@68 61 ",maxRows0="+maxRows0+",columns="+columns+",maxColumns="+maxColumns+",maxColumns0="+maxColumns0;
jpayne@68 62 packed=KillSwitch.allocFloat2D(maxRows+1);
jpayne@68 63 for(int i=0; i<packed.length; i++){
jpayne@68 64 if(i<packed0.length){
jpayne@68 65 packed[i]=packed0[i];
jpayne@68 66 }else{
jpayne@68 67 packed[i]=KillSwitch.allocFloat1D(maxColumns+1);
jpayne@68 68 }
jpayne@68 69 }
jpayne@68 70 //Fill column 0 with insertions
jpayne@68 71 prefillLeftColumnStartingAt(maxRows0);
jpayne@68 72 }
jpayne@68 73 prefillTopRow();
jpayne@68 74 }
jpayne@68 75
jpayne@68 76 /** return new int[] {rows, maxCol, maxState, maxScore, maxStart};
jpayne@68 77 * Will not fill areas that cannot match minScore */
jpayne@68 78 @Override
jpayne@68 79 public final int[] fillLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
jpayne@68 80 return fillUnlimited(read, ref, refStartLoc, refEndLoc, minScore);
jpayne@68 81 }
jpayne@68 82
jpayne@68 83 @Override
jpayne@68 84 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc){
jpayne@68 85 return fillUnlimited(read, ref, refStartLoc, refEndLoc, -999999);
jpayne@68 86 }
jpayne@68 87
jpayne@68 88 /** return new int[] {rows, maxCol, maxState, maxScore, maxStart};
jpayne@68 89 * Min score is optional */
jpayne@68 90 @Override
jpayne@68 91 public final int[] fillUnlimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
jpayne@68 92 initialize(read.length, refEndLoc-refStartLoc+1);
jpayne@68 93
jpayne@68 94 assert(refWeights==null || refWeights.length==ref.length);
jpayne@68 95
jpayne@68 96 //temporary, for finding a bug
jpayne@68 97 if(rows>maxRows || columns>maxColumns){
jpayne@68 98 throw new RuntimeException("rows="+rows+", maxRows="+maxRows+", cols="+columns+", maxCols="+maxColumns+"\n"+new String(read)+"\n");
jpayne@68 99 }
jpayne@68 100
jpayne@68 101 assert(rows<=maxRows) : "Check that values are in-bounds before calling this function: "+rows+", "+maxRows;
jpayne@68 102 assert(columns<=maxColumns) : "Check that values are in-bounds before calling this function: "+columns+", "+maxColumns;
jpayne@68 103
jpayne@68 104 assert(refStartLoc>=0) : "Check that values are in-bounds before calling this function: "+refStartLoc;
jpayne@68 105 assert(refEndLoc<ref.length) : "Check that values are in-bounds before calling this function: "+refEndLoc+", "+ref.length;
jpayne@68 106
jpayne@68 107 final int refOffset=refStartLoc-1;
jpayne@68 108 for(int row=1; row<=rows; row++){
jpayne@68 109
jpayne@68 110 final byte qBase=read[row-1];
jpayne@68 111 for(int col=1; col<=columns; col++){
jpayne@68 112
jpayne@68 113 final int rpos=refOffset+col;
jpayne@68 114 final byte rBase=ref[rpos];
jpayne@68 115
jpayne@68 116 final boolean match=(qBase==rBase);
jpayne@68 117 final boolean defined=(qBase!='N' && rBase!='N');
jpayne@68 118
jpayne@68 119 final float scoreFromDiag=packed[row-1][col-1];
jpayne@68 120 final float scoreFromDel=packed[row][col-1];
jpayne@68 121 final float scoreFromIns=packed[row-1][col];
jpayne@68 122
jpayne@68 123 final float diagScoreM=POINTS_MATCH;
jpayne@68 124 final float diagScoreS=POINTS_SUB;
jpayne@68 125 final float delScore=scoreFromDel+POINTS_DEL;//*delWeights[rpos];
jpayne@68 126 final float insScore=scoreFromIns+POINTS_INS;//*insWeights[rpos];
jpayne@68 127
jpayne@68 128 // assert(delWeights[rpos]==1f) : Arrays.toString(delWeights);
jpayne@68 129 // assert(insWeights[rpos]==1f);
jpayne@68 130 // assert(refWeights[rpos]==1f);
jpayne@68 131
jpayne@68 132 // final int diagScore=scoreFromDiag+(defined ? (match ? diagScoreM : diagScoreS) : POINTS_NOREF);
jpayne@68 133 float diagScore=(match ? diagScoreM : diagScoreS);
jpayne@68 134 diagScore=scoreFromDiag+(defined ? diagScore : POINTS_NOREF)*refWeights[rpos];
jpayne@68 135
jpayne@68 136 float score=diagScore>=delScore ? diagScore : delScore;
jpayne@68 137 score=score>=insScore ? score : insScore;
jpayne@68 138
jpayne@68 139 packed[row][col]=score;
jpayne@68 140 }
jpayne@68 141 //iterationsUnlimited+=columns;
jpayne@68 142 }
jpayne@68 143
jpayne@68 144
jpayne@68 145 int maxCol=-1;
jpayne@68 146 int maxState=-1;
jpayne@68 147 float maxStart=-1;
jpayne@68 148 float maxScore=Integer.MIN_VALUE;
jpayne@68 149
jpayne@68 150 for(int col=1; col<=columns; col++){
jpayne@68 151 float x=packed[rows][col];
jpayne@68 152 if(x>maxScore){
jpayne@68 153 maxScore=x;
jpayne@68 154 maxCol=col;
jpayne@68 155
jpayne@68 156 // assert(rows-1<read.length) : (rows-1)+", "+read.length;
jpayne@68 157 // assert(refOffset+col<ref.length) : refOffset+", "+col+", "+ref.length;
jpayne@68 158 // maxState=getState(rows, col, read[rows-1], ref[refOffset+col], refWeights[refOffset+col], insWeights[refOffset+col], delWeights[refOffset+col]);
jpayne@68 159 maxState=getState(rows, col, read[rows-1], ref[refOffset+col], refWeights[refOffset+col]);
jpayne@68 160 maxStart=x;
jpayne@68 161 }
jpayne@68 162 }
jpayne@68 163
jpayne@68 164 // System.err.println("Returning "+rows+", "+maxCol+", "+maxState+", "+maxScore+"; minScore="+minScore);
jpayne@68 165 return maxScore<minScore ? null : new int[] {rows, maxCol, maxState, (int)maxScore, (int)maxStart};
jpayne@68 166 }
jpayne@68 167
jpayne@68 168 int getState(int row, int col, byte q, byte r, float refWeight, float insWeight, float delWeight){
jpayne@68 169 final boolean match=(q==r);
jpayne@68 170 final boolean defined=(q!='N' && r!='N');
jpayne@68 171
jpayne@68 172 final float scoreFromDiag=packed[row-1][col-1];
jpayne@68 173 final float scoreFromDel=packed[row][col-1];
jpayne@68 174 final float scoreFromIns=packed[row-1][col];
jpayne@68 175
jpayne@68 176 final float diagScoreM=POINTS_MATCH;
jpayne@68 177 final float diagScoreS=POINTS_SUB;
jpayne@68 178 final float delScore=scoreFromDel+POINTS_DEL*delWeight;
jpayne@68 179 final float insScore=scoreFromIns+POINTS_INS*insWeight;
jpayne@68 180
jpayne@68 181 final float diagScore=scoreFromDiag+(defined ? (match ? diagScoreM : diagScoreS) : POINTS_NOREF)*refWeight;
jpayne@68 182
jpayne@68 183 // int score2=diagScore>=delScore ? diagScore : delScore;
jpayne@68 184 // score2=score>=insScore ? score : insScore;
jpayne@68 185
jpayne@68 186 // assert(score==score2) : score+", "+score2;
jpayne@68 187
jpayne@68 188 if(diagScore>=delScore && diagScore>=insScore){
jpayne@68 189 return defined ? match ? MODE_MATCH : MODE_SUB : MODE_N;
jpayne@68 190 }else if(delScore>=insScore){
jpayne@68 191 return MODE_DEL;
jpayne@68 192 }
jpayne@68 193 return MODE_INS;
jpayne@68 194 }
jpayne@68 195
jpayne@68 196 int getState(int row, int col, byte q, byte r, float refWeight){
jpayne@68 197 final boolean match=(q==r);
jpayne@68 198 final boolean defined=(q!='N' && r!='N');
jpayne@68 199
jpayne@68 200 final float scoreFromDiag=packed[row-1][col-1];
jpayne@68 201 final float scoreFromDel=packed[row][col-1];
jpayne@68 202 final float scoreFromIns=packed[row-1][col];
jpayne@68 203
jpayne@68 204 final float diagScoreM=POINTS_MATCH;
jpayne@68 205 final float diagScoreS=POINTS_SUB;
jpayne@68 206 final float delScore=scoreFromDel+POINTS_DEL;
jpayne@68 207 final float insScore=scoreFromIns+POINTS_INS;
jpayne@68 208
jpayne@68 209 final float diagScore=scoreFromDiag+(defined ? (match ? diagScoreM : diagScoreS) : POINTS_NOREF)*refWeight;
jpayne@68 210
jpayne@68 211 // int score2=diagScore>=delScore ? diagScore : delScore;
jpayne@68 212 // score2=score>=insScore ? score : insScore;
jpayne@68 213
jpayne@68 214 // assert(score==score2) : score+", "+score2;
jpayne@68 215
jpayne@68 216 if(diagScore>=delScore && diagScore>=insScore){
jpayne@68 217 return defined ? match ? MODE_MATCH : MODE_SUB : MODE_N;
jpayne@68 218 }else if(delScore>=insScore){
jpayne@68 219 return MODE_DEL;
jpayne@68 220 }
jpayne@68 221 return MODE_INS;
jpayne@68 222 }
jpayne@68 223
jpayne@68 224
jpayne@68 225 /** Generates the match string.
jpayne@68 226 * State is NOT used. */
jpayne@68 227 @Override
jpayne@68 228 public final byte[] traceback(byte[] query, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state){
jpayne@68 229 // assert(false);
jpayne@68 230 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
jpayne@68 231 assert(row==rows);
jpayne@68 232
jpayne@68 233 byte[] out=new byte[row+col-1]; //TODO if an out of bound crash occurs, try removing the "-1".
jpayne@68 234 int outPos=0;
jpayne@68 235
jpayne@68 236 // assert(state==(packed[row][col]&MODEMASK));
jpayne@68 237
jpayne@68 238 while(row>0 && col>0){
jpayne@68 239 byte q=query[row-1];
jpayne@68 240 int rpos=refStartLoc+col-1;
jpayne@68 241 byte r=ref[rpos];
jpayne@68 242 boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
jpayne@68 243 // state=getState(row, col, q, r, refWeights[rpos], insWeights[rpos], delWeights[rpos]);
jpayne@68 244 state=getState(row, col, q, r, refWeights[rpos]);
jpayne@68 245 if(state==MODE_MATCH){
jpayne@68 246 col--;
jpayne@68 247 row--;
jpayne@68 248 out[outPos]=defined ? (byte)'m' : (byte)'N';
jpayne@68 249 }else if(state==MODE_SUB){
jpayne@68 250 col--;
jpayne@68 251 row--;
jpayne@68 252 out[outPos]=defined ? (byte)'S' : (byte)'N';
jpayne@68 253 }else if(state==MODE_N){
jpayne@68 254 col--;
jpayne@68 255 row--;
jpayne@68 256 out[outPos]='N';
jpayne@68 257 }else if(state==MODE_DEL){
jpayne@68 258 col--;
jpayne@68 259 out[outPos]='D';
jpayne@68 260 }else if(state==MODE_INS){
jpayne@68 261 row--;
jpayne@68 262 // out[outPos]='I';
jpayne@68 263 // out[outPos]=(col<0 || col>=columns ? (byte)'C' : (byte)'I');
jpayne@68 264 if(col>=0 && col<columns){
jpayne@68 265 out[outPos]='I';
jpayne@68 266 }else{
jpayne@68 267 out[outPos]='C';
jpayne@68 268 col--;
jpayne@68 269 }
jpayne@68 270 }else{
jpayne@68 271 assert(false) : state;
jpayne@68 272 }
jpayne@68 273 outPos++;
jpayne@68 274 }
jpayne@68 275
jpayne@68 276 assert(row==0 || col==0);
jpayne@68 277 if(col!=row){
jpayne@68 278 while(row>0){
jpayne@68 279 out[outPos]='C';
jpayne@68 280 outPos++;
jpayne@68 281 row--;
jpayne@68 282 col--;
jpayne@68 283 }
jpayne@68 284 if(col>0){
jpayne@68 285 //do nothing
jpayne@68 286 }
jpayne@68 287 }
jpayne@68 288
jpayne@68 289 //Shrink and reverse the string
jpayne@68 290 byte[] out2=new byte[outPos];
jpayne@68 291 for(int i=0; i<outPos; i++){
jpayne@68 292 out2[i]=out[outPos-i-1];
jpayne@68 293 }
jpayne@68 294 out=null;
jpayne@68 295
jpayne@68 296 return out2;
jpayne@68 297 }
jpayne@68 298
jpayne@68 299 @Override
jpayne@68 300 /** Generates identity;
jpayne@68 301 * fills 'extra' with {match, sub, del, ins, N, clip} if present */
jpayne@68 302 public float tracebackIdentity(byte[] query, byte[] ref, int refStartLoc, int refEndLoc, int row, int col, int state, int[] extra){
jpayne@68 303
jpayne@68 304 // assert(false);
jpayne@68 305 assert(refStartLoc<=refEndLoc) : refStartLoc+", "+refEndLoc;
jpayne@68 306 assert(row==rows);
jpayne@68 307
jpayne@68 308 // assert(state==(packed[row][col]&MODEMASK));
jpayne@68 309 int match=0, sub=0, del=0, ins=0, noref=0, clip=0;
jpayne@68 310
jpayne@68 311 while(row>0 && col>0){
jpayne@68 312 byte q=query[row-1];
jpayne@68 313 int rpos=refStartLoc+col-1;
jpayne@68 314 byte r=ref[rpos];
jpayne@68 315 boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
jpayne@68 316 // state=getState(row, col, q, r, refWeights[rpos], insWeights[rpos], delWeights[rpos]);
jpayne@68 317 state=getState(row, col, q, r, refWeights[rpos]);
jpayne@68 318 if(state==MODE_MATCH){
jpayne@68 319 col--;
jpayne@68 320 row--;
jpayne@68 321 match+=(defined ? 1 : 0);
jpayne@68 322 noref+=(defined ? 0 : 1);
jpayne@68 323 }else if(state==MODE_SUB){
jpayne@68 324 col--;
jpayne@68 325 row--;
jpayne@68 326 sub+=(defined ? 1 : 0);
jpayne@68 327 noref+=(defined ? 0 : 1);
jpayne@68 328 }else if(state==MODE_N){
jpayne@68 329 col--;
jpayne@68 330 row--;
jpayne@68 331 noref++;
jpayne@68 332 }else if(state==MODE_DEL){
jpayne@68 333 col--;
jpayne@68 334 del++;
jpayne@68 335 }else if(state==MODE_INS){
jpayne@68 336 row--;
jpayne@68 337 boolean edge=(col<=1 || col>=columns);
jpayne@68 338 ins+=(edge ? 0 : 1);
jpayne@68 339 clip+=(edge ? 1 : 0);
jpayne@68 340 }else{
jpayne@68 341 assert(false) : state;
jpayne@68 342 }
jpayne@68 343 }
jpayne@68 344
jpayne@68 345 assert(row==0 || col==0);
jpayne@68 346 if(col!=row){//Not sure what this is doing
jpayne@68 347 while(row>0){
jpayne@68 348 clip++;
jpayne@68 349 row--;
jpayne@68 350 col--;
jpayne@68 351 }
jpayne@68 352 if(col>0){
jpayne@68 353 //do nothing
jpayne@68 354 }
jpayne@68 355 }
jpayne@68 356
jpayne@68 357 if(extra!=null){
jpayne@68 358 assert(extra.length==5);
jpayne@68 359 extra[0]=match;
jpayne@68 360 extra[1]=sub;
jpayne@68 361 extra[2]=del;
jpayne@68 362 extra[3]=ins;
jpayne@68 363 extra[4]=noref;
jpayne@68 364 extra[5]=clip;
jpayne@68 365 }
jpayne@68 366
jpayne@68 367 float len=match+sub+ins+del+noref*0.1f;
jpayne@68 368 float id=match/Tools.max(1.0f, len);
jpayne@68 369 return id;
jpayne@68 370 }
jpayne@68 371
jpayne@68 372 /** @return {score, bestRefStart, bestRefStop} */
jpayne@68 373 @Override
jpayne@68 374 public final int[] score(final byte[] read, final byte[] ref, final int refStartLoc, final int refEndLoc,
jpayne@68 375 final int maxRow, final int maxCol, final int maxState/*, final int maxScore, final int maxStart*/){
jpayne@68 376
jpayne@68 377 int row=maxRow;
jpayne@68 378 int col=maxCol;
jpayne@68 379 int state=maxState;
jpayne@68 380
jpayne@68 381 assert(maxState>=0 && maxState<packed.length) :
jpayne@68 382 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
jpayne@68 383 assert(maxRow>=0 && maxRow<packed.length) :
jpayne@68 384 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
jpayne@68 385 assert(maxCol>=0 && maxCol<packed[0].length) :
jpayne@68 386 maxState+", "+maxRow+", "+maxCol+"\n"+new String(read)+"\n"+toString(ref, refStartLoc, refEndLoc);
jpayne@68 387
jpayne@68 388 float score=packed[maxRow][maxCol]; //Or zero, if it is to be recalculated
jpayne@68 389
jpayne@68 390 if(row<rows){
jpayne@68 391 int difR=rows-row;
jpayne@68 392 int difC=columns-col;
jpayne@68 393
jpayne@68 394 while(difR>difC){
jpayne@68 395 score+=POINTS_NOREF;
jpayne@68 396 difR--;
jpayne@68 397 }
jpayne@68 398
jpayne@68 399 row+=difR;
jpayne@68 400 col+=difR;
jpayne@68 401
jpayne@68 402 }
jpayne@68 403
jpayne@68 404 assert(refStartLoc<=refEndLoc);
jpayne@68 405 assert(row==rows);
jpayne@68 406
jpayne@68 407
jpayne@68 408 final int bestRefStop=refStartLoc+col-1;
jpayne@68 409
jpayne@68 410 while(row>0 && col>0){
jpayne@68 411 final byte q=read[row-1];
jpayne@68 412 int rpos=refStartLoc+col-1;
jpayne@68 413 final byte r=ref[rpos];
jpayne@68 414 // final boolean defined=(AminoAcid.isFullyDefined(q) && AminoAcid.isFullyDefined(r));
jpayne@68 415 // state=getState(row, col, q, r, refWeights[rpos], insWeights[rpos], delWeights[rpos]);
jpayne@68 416 state=getState(row, col, q, r, refWeights[rpos]);
jpayne@68 417 if(state==MODE_MATCH){
jpayne@68 418 col--;
jpayne@68 419 row--;
jpayne@68 420 }else if(state==MODE_SUB){
jpayne@68 421 col--;
jpayne@68 422 row--;
jpayne@68 423 }else if(state==MODE_N){
jpayne@68 424 col--;
jpayne@68 425 row--;
jpayne@68 426 }else if(state==MODE_DEL){
jpayne@68 427 col--;
jpayne@68 428 }else if(state==MODE_INS){
jpayne@68 429 row--;
jpayne@68 430 }else{
jpayne@68 431 assert(false) : state;
jpayne@68 432 }
jpayne@68 433 }
jpayne@68 434 // assert(false) : row+", "+col;
jpayne@68 435 if(row>col){
jpayne@68 436 col-=row;
jpayne@68 437 }
jpayne@68 438
jpayne@68 439 final int bestRefStart=refStartLoc+col;
jpayne@68 440
jpayne@68 441 // System.err.println("t2\t"+score+", "+maxScore+", "+maxStart+", "+bestRefStart);
jpayne@68 442 int[] rvec;
jpayne@68 443 if(bestRefStart<refStartLoc || bestRefStop>refEndLoc){ //Suggest extra padding in cases of overflow
jpayne@68 444 int padLeft=Tools.max(0, refStartLoc-bestRefStart);
jpayne@68 445 int padRight=Tools.max(0, bestRefStop-refEndLoc);
jpayne@68 446 rvec=new int[] {(int)score, bestRefStart, bestRefStop, padLeft, padRight};
jpayne@68 447 }else{
jpayne@68 448 rvec=new int[] {(int)score, bestRefStart, bestRefStop};
jpayne@68 449 }
jpayne@68 450 return rvec;
jpayne@68 451 }
jpayne@68 452
jpayne@68 453
jpayne@68 454 /** Will not fill areas that cannot match minScore.
jpayne@68 455 * @return {score, bestRefStart, bestRefStop} */
jpayne@68 456 @Override
jpayne@68 457 public final int[] fillAndScoreLimited(byte[] read, byte[] ref, int refStartLoc, int refEndLoc, int minScore){
jpayne@68 458 int a=Tools.max(0, refStartLoc);
jpayne@68 459 int b=Tools.min(ref.length-1, refEndLoc);
jpayne@68 460 assert(b>=a);
jpayne@68 461
jpayne@68 462 if(b-a>=maxColumns){
jpayne@68 463 System.err.println("Warning: Max alignment columns exceeded; restricting range. "+(b-a+1)+" > "+maxColumns);
jpayne@68 464 assert(false) : refStartLoc+", "+refEndLoc;
jpayne@68 465 b=Tools.min(ref.length-1, a+maxColumns-1);
jpayne@68 466 }
jpayne@68 467 int[] max=fillLimited(read, ref, a, b, minScore);
jpayne@68 468 // return max==null ? null : new int[] {max[3], 0, max[1]};
jpayne@68 469
jpayne@68 470 int[] score=(max==null ? null : score(read, ref, a, b, max[0], max[1], max[2]/*, max[3], max[4]*/));
jpayne@68 471
jpayne@68 472 return score;
jpayne@68 473 }
jpayne@68 474
jpayne@68 475 public static final String toString(byte[] ref, int startLoc, int stopLoc){
jpayne@68 476 StringBuilder sb=new StringBuilder(stopLoc-startLoc+1);
jpayne@68 477 for(int i=startLoc; i<=stopLoc; i++){sb.append((char)ref[i]);}
jpayne@68 478 return sb.toString();
jpayne@68 479 }
jpayne@68 480
jpayne@68 481 // public static int calcDelScore(int len){
jpayne@68 482 // if(len<=0){return 0;}
jpayne@68 483 // int score=POINTS_DEL;
jpayne@68 484 // if(len>1){
jpayne@68 485 // score+=(len-1)*POINTS_DEL2;
jpayne@68 486 // }
jpayne@68 487 // return score;
jpayne@68 488 // }
jpayne@68 489
jpayne@68 490 // public int maxScoreByIdentity(int len, float identity){
jpayne@68 491 // assert(identity>=0 && identity<=1);
jpayne@68 492 // return (int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
jpayne@68 493 // }
jpayne@68 494
jpayne@68 495 @Override
jpayne@68 496 public int minScoreByIdentity(int len, float identity){
jpayne@68 497 assert(identity>=0 && identity<=1);
jpayne@68 498
jpayne@68 499 int a=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_SUB));
jpayne@68 500 int b=(int)(len*(identity*POINTS_MATCH+(1-identity)*POINTS_INS));
jpayne@68 501 int c=(int)(len*(1*POINTS_MATCH+((1/(Tools.max(identity, 0.000001f)))-1)*POINTS_DEL));
jpayne@68 502 return Tools.min(a, b, c);
jpayne@68 503 }
jpayne@68 504
jpayne@68 505 private static float calcDelScore(int len){
jpayne@68 506 if(len<=0){return 0;}
jpayne@68 507 float score=POINTS_DEL*len;
jpayne@68 508 return score;
jpayne@68 509 }
jpayne@68 510 //
jpayne@68 511 // public static int calcInsScore(int len){
jpayne@68 512 // if(len<=0){return 0;}
jpayne@68 513 // int score=POINTS_INS;
jpayne@68 514 //
jpayne@68 515 // if(len>1){
jpayne@68 516 // score+=(len-1)*POINTS_INS2;
jpayne@68 517 // }
jpayne@68 518 // return score;
jpayne@68 519 // }
jpayne@68 520 //
jpayne@68 521 // private static int calcInsScoreOffset(int len){
jpayne@68 522 // if(len<=0){return 0;}
jpayne@68 523 // int score=POINTS_INS;
jpayne@68 524 //
jpayne@68 525 // if(len>1){
jpayne@68 526 // score+=(len-1)*POINTS_INS2;
jpayne@68 527 // }
jpayne@68 528 // return score;
jpayne@68 529 // }
jpayne@68 530
jpayne@68 531 @Override
jpayne@68 532 public int rows(){return rows;}
jpayne@68 533 @Override
jpayne@68 534 public int columns(){return columns;}
jpayne@68 535
jpayne@68 536 public void setWeights(float[] refWeights_, float[] insWeights_, float[] delWeights_){
jpayne@68 537 refWeights=refWeights_;
jpayne@68 538 // insWeights=insWeights_;
jpayne@68 539 // delWeights=delWeights_;
jpayne@68 540 }
jpayne@68 541
jpayne@68 542 public void setWeights(float[] refWeights_){
jpayne@68 543 refWeights=refWeights_;
jpayne@68 544 }
jpayne@68 545
jpayne@68 546 private int maxRows;
jpayne@68 547 private int maxColumns;
jpayne@68 548
jpayne@68 549 private float[][] packed;
jpayne@68 550
jpayne@68 551 private float[] refWeights;
jpayne@68 552 //These don't seem to help.
jpayne@68 553 // private float[] insWeights;
jpayne@68 554 // private float[] delWeights;
jpayne@68 555
jpayne@68 556 public static final int MAX_SCORE=Integer.MAX_VALUE-2000;
jpayne@68 557 public static final int MIN_SCORE=0-MAX_SCORE; //Keeps it 1 point above "BAD".
jpayne@68 558
jpayne@68 559 //For some reason changing MODE_DEL from 1 to 0 breaks everything
jpayne@68 560 private static final byte MODE_DEL=1;
jpayne@68 561 private static final byte MODE_INS=2;
jpayne@68 562 private static final byte MODE_SUB=3;
jpayne@68 563 private static final byte MODE_MATCH=4;
jpayne@68 564 private static final byte MODE_N=5;
jpayne@68 565
jpayne@68 566 public static final float POINTS_NOREF=-20;
jpayne@68 567 public static final float POINTS_MATCH=100;
jpayne@68 568 public static final float POINTS_SUB=-50;
jpayne@68 569 public static final float POINTS_INS=-121;
jpayne@68 570 public static final float POINTS_DEL=-111;
jpayne@68 571
jpayne@68 572 public static final int BAD=MIN_SCORE-1;
jpayne@68 573
jpayne@68 574 private int rows;
jpayne@68 575 private int columns;
jpayne@68 576
jpayne@68 577 // public long iterationsLimited=0;
jpayne@68 578 // public long iterationsUnlimited=0;
jpayne@68 579
jpayne@68 580 public boolean verbose=false;
jpayne@68 581 public boolean verbose2=false;
jpayne@68 582
jpayne@68 583 }