Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/icecream/IceCreamAlignerJNI.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 icecream;
jpayne@68 2
jpayne@68 3 import aligner.AlignmentResult;
jpayne@68 4 import shared.KillSwitch;
jpayne@68 5 import shared.Shared;
jpayne@68 6
jpayne@68 7 public final class IceCreamAlignerJNI extends IceCreamAligner {
jpayne@68 8
jpayne@68 9 static {
jpayne@68 10 Shared.loadJNI();
jpayne@68 11 }
jpayne@68 12
jpayne@68 13 IceCreamAlignerJNI(){}
jpayne@68 14
jpayne@68 15 /**
jpayne@68 16 * @param query
jpayne@68 17 * @param ref
jpayne@68 18 * @param qstart
jpayne@68 19 * @param rstart
jpayne@68 20 * @param rstop
jpayne@68 21 * @param minScore Quit early if score drops below this
jpayne@68 22 * @param minRatio Don't return results if max score is less than this fraction of max possible score
jpayne@68 23 * @return
jpayne@68 24 */
jpayne@68 25 @Override
jpayne@68 26 public AlignmentResult alignForward(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore,
jpayne@68 27 final float minRatio) {
jpayne@68 28
jpayne@68 29 final int qlen=query.length;
jpayne@68 30 final int rlen=rstop-rstart+1;
jpayne@68 31 final int[] retVec=KillSwitch.allocInt1D(4);
jpayne@68 32
jpayne@68 33 // alignForwardJNI(qInt, rInt, retVec, qlen, rlen, minScore, minRatio);
jpayne@68 34 // alignForwardPseudo(qInt, rInt, retVec, qlen, rlen, minScore, minRatio);
jpayne@68 35
jpayne@68 36 if((qlen+rlen+32)*2>Short.MAX_VALUE) {
jpayne@68 37 final int[] qInt=new int[query.length];
jpayne@68 38 final int[] rInt=new int[rlen];
jpayne@68 39 for(int i=0; i<query.length; i++){qInt[i]=query[i];}
jpayne@68 40 for(int i=0; i<rlen; i++){rInt[i]=ref[i+rstart];}
jpayne@68 41
jpayne@68 42 // alignForwardPseudo(qInt, rInt, retVec, qlen, rlen, minScore, minRatio);
jpayne@68 43 alignForwardJNI(qInt, rInt, retVec, qlen, rlen, minScore, minRatio);
jpayne@68 44 // alignForwardShortJNI(qInt, rInt, retVec, qlen, rlen);
jpayne@68 45 }else{
jpayne@68 46 final short[] qInt=new short[query.length];
jpayne@68 47 final short[] rInt=new short[rlen];
jpayne@68 48 for(int i=0; i<query.length; i++){qInt[i]=query[i];}
jpayne@68 49 for(int i=0; i<rlen; i++){rInt[i]=ref[i+rstart];}
jpayne@68 50
jpayne@68 51 alignForward16JNI(qInt, rInt, retVec, (short)qlen, (short)rlen, (short)minScore, minRatio);
jpayne@68 52 // alignForwardShort16JNI(qInt, rInt, retVec, (short)qlen, (short)rlen);
jpayne@68 53 }
jpayne@68 54
jpayne@68 55 int maxScore=retVec[0];
jpayne@68 56 final int maxQpos=retVec[1];
jpayne@68 57 final int maxRpos=retVec[2]+rstart;
jpayne@68 58 final int jniIters=retVec[3];
jpayne@68 59
jpayne@68 60 iters+=jniIters;
jpayne@68 61
jpayne@68 62 maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
jpayne@68 63 final float ratio=maxScore/(float)query.length;
jpayne@68 64
jpayne@68 65 if(ratio<minRatio){return null;}
jpayne@68 66
jpayne@68 67 return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
jpayne@68 68 }
jpayne@68 69
jpayne@68 70 /**
jpayne@68 71 * @param query
jpayne@68 72 * @param ref
jpayne@68 73 * @param qstart
jpayne@68 74 * @param rstart
jpayne@68 75 * @param rstop
jpayne@68 76 * @param minScore Quit early if score drops below this
jpayne@68 77 * @param minRatio Don't return results if max score is less than this fraction of max possible score
jpayne@68 78 * @return
jpayne@68 79 */
jpayne@68 80 @Override
jpayne@68 81 public AlignmentResult alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore,
jpayne@68 82 final float minRatio) {
jpayne@68 83
jpayne@68 84 final int qlen=query.length;
jpayne@68 85 final int rlen=rstop-rstart+1;
jpayne@68 86 final int[] retVec=KillSwitch.allocInt1D(4);
jpayne@68 87
jpayne@68 88 if((qlen+rlen+32)*2>Short.MAX_VALUE) {
jpayne@68 89 final int[] qInt=new int[query.length];
jpayne@68 90 final int[] rInt=new int[rlen];
jpayne@68 91 for(int i=0; i<query.length; i++){qInt[i]=query[i];}
jpayne@68 92 for(int i=0; i<rlen; i++){rInt[i]=ref[i+rstart];}
jpayne@68 93
jpayne@68 94 // alignForwardShortPseudo(qInt, rInt, retVec, qlen, rlen);
jpayne@68 95 alignForwardShortJNI(qInt, rInt, retVec, qlen, rlen);
jpayne@68 96 }else{
jpayne@68 97 final short[] qInt=new short[query.length];
jpayne@68 98 final short[] rInt=new short[rlen];
jpayne@68 99 for(int i=0; i<query.length; i++){qInt[i]=query[i];}
jpayne@68 100 for(int i=0; i<rlen; i++){rInt[i]=ref[i+rstart];}
jpayne@68 101
jpayne@68 102 alignForwardShort16JNI(qInt, rInt, retVec, (short)qlen, (short)rlen);
jpayne@68 103 }
jpayne@68 104 int maxScore=retVec[0];
jpayne@68 105 final int maxQpos=retVec[1];
jpayne@68 106 final int maxRpos=retVec[2]+rstart;
jpayne@68 107 final int jniIters=retVec[3];
jpayne@68 108
jpayne@68 109 itersShort+=jniIters;
jpayne@68 110
jpayne@68 111 maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
jpayne@68 112 final float ratio=maxScore/(float)query.length;
jpayne@68 113
jpayne@68 114 if(ratio<minRatio){return null;}
jpayne@68 115
jpayne@68 116 return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
jpayne@68 117 }
jpayne@68 118
jpayne@68 119 /*--------------------------------------------------------------*/
jpayne@68 120 /*---------------- JNI ----------------*/
jpayne@68 121 /*--------------------------------------------------------------*/
jpayne@68 122
jpayne@68 123 private static void alignForwardPseudo(final int[] query, final int[] ref, final int[] retArray,
jpayne@68 124 final int qlen, final int rlen, final int minScore, final float minRatio) {
jpayne@68 125
jpayne@68 126 final int arrayLength=rlen;
jpayne@68 127 final int arrayLength2=rlen+1;
jpayne@68 128
jpayne@68 129 //Stack allocated; faster.
jpayne@68 130 int[] array1=new int[arrayLength2];
jpayne@68 131 int[] array2=new int[arrayLength2];
jpayne@68 132
jpayne@68 133 int[] prev=array1;
jpayne@68 134 int[] next=array2;
jpayne@68 135 int maxScore=-32000;
jpayne@68 136 int maxQpos=-1;
jpayne@68 137 int maxRpos=-1;
jpayne@68 138 int iters=0;
jpayne@68 139
jpayne@68 140 final int minPassingScore=(int)(qlen*minRatio*pointsMatch);
jpayne@68 141 final int minPassingScore3=minPassingScore-qlen*pointsMatch;
jpayne@68 142
jpayne@68 143 for(int i=0; i<=arrayLength-qlen; i++){prev[i]=0;}
jpayne@68 144 for(int i=arrayLength-qlen, score=0; i<=arrayLength; i++, score+=pointsDel) {
jpayne@68 145 prev[i]=score;
jpayne@68 146 next[i]=0;
jpayne@68 147 }
jpayne@68 148
jpayne@68 149 for(int qpos=0; qpos<qlen; qpos++){
jpayne@68 150 prev[0]=pointsIns*qpos;
jpayne@68 151
jpayne@68 152 final int q=query[qpos];
jpayne@68 153 int maxScoreThisPass=-32000;
jpayne@68 154 final int remainingBases=(qlen-qpos);
jpayne@68 155 final int remainingPoints=remainingBases*pointsMatch;
jpayne@68 156 final int minViableScore=minPassingScore3-remainingPoints;
jpayne@68 157
jpayne@68 158 // for(int rpos=0, apos=1; rpos<rlen; rpos++, apos++){
jpayne@68 159 // final int r=ref[rpos];
jpayne@68 160 // final boolean match=(q==r);
jpayne@68 161 // final int vScore=prev[apos]+pointsIns;
jpayne@68 162 // final int hScore=next[apos-1]+pointsDel;
jpayne@68 163 // final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68 164 //
jpayne@68 165 // //Should be branchless conditional moves
jpayne@68 166 // int score=(dScore>=vScore ? dScore : vScore);
jpayne@68 167 // score=(hScore>score ? hScore : score);
jpayne@68 168 // next[apos]=score;
jpayne@68 169 // maxScoreThisPass=(score>maxScoreThisPass ? score : maxScoreThisPass);
jpayne@68 170 // }
jpayne@68 171
jpayne@68 172 for(int rpos=0, apos=1; rpos<rlen; rpos++, apos++){
jpayne@68 173 final int r=ref[rpos];
jpayne@68 174 final boolean match=(q==r);
jpayne@68 175 final int vScore=prev[apos]+pointsIns;
jpayne@68 176 final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68 177
jpayne@68 178 //Should be branchless conditional moves
jpayne@68 179 int score=(dScore>=vScore ? dScore : vScore);
jpayne@68 180 next[apos]=score;
jpayne@68 181 }
jpayne@68 182
jpayne@68 183 for(int apos=1; apos<arrayLength2; apos++){
jpayne@68 184 final int hScore=next[apos-1]+pointsDel;
jpayne@68 185
jpayne@68 186 //Should be branchless conditional moves
jpayne@68 187 int score=next[apos];
jpayne@68 188 score=(hScore>score ? hScore : score);
jpayne@68 189 next[apos]=score;
jpayne@68 190 maxScoreThisPass=(score>maxScoreThisPass ? score : maxScoreThisPass);
jpayne@68 191 }
jpayne@68 192 iters+=arrayLength;
jpayne@68 193
jpayne@68 194 //Aggressive early exit
jpayne@68 195 if(maxScoreThisPass<minScore){return;}
jpayne@68 196
jpayne@68 197 //Safe early exit
jpayne@68 198 if(maxScoreThisPass<minViableScore){return;}
jpayne@68 199
jpayne@68 200 int[] temp=prev;
jpayne@68 201 prev=next;
jpayne@68 202 next=temp;
jpayne@68 203 }
jpayne@68 204
jpayne@68 205 maxScore=-32000;
jpayne@68 206 for(int apos=1; apos<arrayLength2; apos++){//Grab high score from last iteration
jpayne@68 207 int score=prev[apos];
jpayne@68 208 if(score>=maxScore){
jpayne@68 209 maxScore=score;
jpayne@68 210 maxQpos=qlen;
jpayne@68 211 maxRpos=apos-1;
jpayne@68 212 }
jpayne@68 213 }
jpayne@68 214 retArray[0]=maxScore;
jpayne@68 215 retArray[1]=maxQpos;
jpayne@68 216 retArray[2]=maxRpos;
jpayne@68 217 retArray[3]=iters;
jpayne@68 218 }
jpayne@68 219
jpayne@68 220 private static void alignForwardShortPseudo(int[] query, int[] ref, int[] retArray, int qlen, int rlen) {
jpayne@68 221
jpayne@68 222 final int arrayLength=qlen;
jpayne@68 223 final int arrayLength2=qlen+1;
jpayne@68 224
jpayne@68 225 //Stack allocated; faster.
jpayne@68 226 int[] array1=new int[arrayLength2];
jpayne@68 227 int[] array2=new int[arrayLength2];
jpayne@68 228
jpayne@68 229 int[] prev=array1;
jpayne@68 230 int[] next=array2;
jpayne@68 231 int maxScore=-999999;
jpayne@68 232 int maxQpos=-1;
jpayne@68 233 int maxRpos=-1;
jpayne@68 234 int itersShort=0;
jpayne@68 235
jpayne@68 236 for(int i=0; i<arrayLength2; i++) {
jpayne@68 237 prev[i]=pointsIns*i;
jpayne@68 238 next[i]=0;//For C version initialization
jpayne@68 239 }
jpayne@68 240
jpayne@68 241 for(int rpos=0; rpos<rlen; rpos++){
jpayne@68 242 if(rlen-rpos<arrayLength){prev[0]=next[0]+pointsDel;}
jpayne@68 243
jpayne@68 244 final int r=ref[rpos];
jpayne@68 245 int score=0;
jpayne@68 246
jpayne@68 247 // //Inner loop
jpayne@68 248 // for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68 249 // final int q=query[qpos];
jpayne@68 250 // final boolean match=(q==r);
jpayne@68 251 // final int vScore=prev[apos]+pointsIns;
jpayne@68 252 // final int hScore=next[apos-1]+pointsDel;
jpayne@68 253 // final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68 254 //
jpayne@68 255 // score=(dScore>=vScore ? dScore : vScore);
jpayne@68 256 // score=(hScore>score ? hScore : score);
jpayne@68 257 //
jpayne@68 258 // next[apos]=score;
jpayne@68 259 // }
jpayne@68 260
jpayne@68 261 //Inner DV loop
jpayne@68 262 for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68 263 final int q=query[qpos];
jpayne@68 264 final boolean match=(q==r);
jpayne@68 265 final int vScore=prev[apos]+pointsIns;
jpayne@68 266 final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68 267
jpayne@68 268 score=(dScore>=vScore ? dScore : vScore);
jpayne@68 269
jpayne@68 270 next[apos]=score;
jpayne@68 271 }
jpayne@68 272
jpayne@68 273 //Inner I loop
jpayne@68 274 for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68 275 final int hScore=next[apos-1]+pointsDel;
jpayne@68 276
jpayne@68 277 score=next[apos];
jpayne@68 278 score=(hScore>score ? hScore : score);
jpayne@68 279
jpayne@68 280 next[apos]=score;
jpayne@68 281 }
jpayne@68 282
jpayne@68 283 itersShort+=arrayLength;
jpayne@68 284
jpayne@68 285 if(score>=maxScore){
jpayne@68 286 maxScore=score;
jpayne@68 287 maxQpos=arrayLength-1;
jpayne@68 288 maxRpos=rpos;
jpayne@68 289 }
jpayne@68 290
jpayne@68 291 int[] temp=prev;
jpayne@68 292 prev=next;
jpayne@68 293 next=temp;
jpayne@68 294 }
jpayne@68 295 retArray[0]=maxScore;
jpayne@68 296 retArray[1]=maxQpos;
jpayne@68 297 retArray[2]=maxRpos;
jpayne@68 298 retArray[3]=itersShort;
jpayne@68 299 }
jpayne@68 300
jpayne@68 301 private static native void alignForwardJNI(int[] query, int[] ref, int[] retArray, int qlen, int rlen, int minScore, float minRatio);
jpayne@68 302 private static native void alignForward16JNI(short[] query, short[] ref, int[] retArray, short qlen, short rlen, short minScore, float minRatio);
jpayne@68 303 private static native void alignForwardShortJNI(int[] query, int[] ref, int[] retArray, int qlen, int rlen);
jpayne@68 304 private static native void alignForwardShort16JNI(short[] query, short[] ref, int[] retArray, short qlen, short rlen);
jpayne@68 305
jpayne@68 306 /*--------------------------------------------------------------*/
jpayne@68 307 /*---------------- Getters ----------------*/
jpayne@68 308 /*--------------------------------------------------------------*/
jpayne@68 309
jpayne@68 310 @Override
jpayne@68 311 long iters(){return iters;}
jpayne@68 312
jpayne@68 313 @Override
jpayne@68 314 long itersShort(){return itersShort;}
jpayne@68 315
jpayne@68 316 /*--------------------------------------------------------------*/
jpayne@68 317 /*---------------- Fields ----------------*/
jpayne@68 318 /*--------------------------------------------------------------*/
jpayne@68 319
jpayne@68 320 long iters = 0;
jpayne@68 321 long itersShort = 0;
jpayne@68 322
jpayne@68 323 /*--------------------------------------------------------------*/
jpayne@68 324 /*---------------- Constants ----------------*/
jpayne@68 325 /*--------------------------------------------------------------*/
jpayne@68 326
jpayne@68 327 public static final int pointsMatch = 1;
jpayne@68 328 public static final int pointsSub = -1;
jpayne@68 329 public static final int pointsDel = -2;
jpayne@68 330 public static final int pointsIns = -2;
jpayne@68 331
jpayne@68 332 }