Mercurial > repos > rliterman > csp2

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
comparison
equal deleted inserted replaced
67:0e9998148a16 68:5028fdace37b
1 package sketch;
2
3 import java.util.Comparator;
4 import java.util.Locale;
5
6 import aligner.Aligner;
7 import prok.GeneCaller;
8 import shared.Tools;
9 import tax.TaxNode;
10
11 public final class Comparison extends SketchObject implements Comparable<Comparison> {
12
13 /*--------------------------------------------------------------*/
14 /*---------------- Constructors ----------------*/
15 /*--------------------------------------------------------------*/
16
17 // public Comparison(CompareBuffer buffer){
18 // this(buffer, null, null);
19 // }
20
21 // public Comparison(Sketch a_, Sketch b_){
22 // this(null, a_, b_);
23 // }
24
25 public Comparison(CompareBuffer buffer, Sketch a_, Sketch b_){
26
27 a=a_;
28 b=b_;
29
30 if(buffer!=null){setFrom(buffer);}
31
32 if(b!=null){
33 taxName=b.taxName();
34 taxID=b.taxID;
35 }
36
37 // System.err.println(this);
38 // System.err.println(b.present);
39 }
40
41 /*--------------------------------------------------------------*/
42 /*---------------- Mutators ----------------*/
43 /*--------------------------------------------------------------*/
44
45 public void setFrom(CompareBuffer buffer){
46 hits=buffer.hits();
47 multiHits=buffer.multiHits();
48 unique2=buffer.unique2();
49 unique3=buffer.unique3();
50 noHits=buffer.noHits();
51
52 contamHits=buffer.contamHits();
53 contam2Hits=buffer.contam2Hits();
54 multiContamHits=buffer.multiContamHits();
55
56 refDivisor=buffer.refDivisor();
57 queryDivisor=buffer.queryDivisor();
58
59 refSize=buffer.refSize();
60 querySize=buffer.querySize();
61
62 depth=buffer.depth();
63 depth2=buffer.depth2();
64 float x=buffer.avgRefHits();
65 if(x>0){avgRefHits=x;}
66 // volume=volume0();
67 score=score0();
68
69 hits1=buffer.hits1();
70 qSeen1=buffer.qSeen1();
71 rSeen1=buffer.rSeen1();
72 }
73
74 public void recompare(CompareBuffer buffer, int[][] taxHits, int contamLevel){
75
76 // for(int[] row : taxHits){
77 // if(row!=null){
78 // System.err.println(Arrays.toString(row));
79 // }
80 // }//Tested; correctly indicates most rows have octopus but some have other things.
81
82 assert(a.merged());
83 // int oldContam2=contam2Hits;
84 int x=a.countMatches(b, buffer, a.compareBitSet(), false, taxHits, contamLevel);
85 assert(x==hits);
86 setFrom(buffer);
87 // contam2Hits=oldContam2;
88 }
89
90 /*--------------------------------------------------------------*/
91 /*---------------- Methods ----------------*/
92 /*--------------------------------------------------------------*/
93
94 @Override
95 public boolean equals(Object b){
96 if(b==null || b.getClass()!=this.getClass()){return false;}
97 return scoreComparator.compare(this, (Comparison)b)==0;
98 }
99
100 // //WKID
101 // public float wkid(){return idMinDivisor();}
102 // public float idMinDivisor(){
103 // return hits/(float)minDivisor();
104 // }
105
106 // public float k1Fraction(){
107 // return a.k1Fraction();
108 // }
109
110 // //KID
111 // public float kid(){return idMaxDivisor();}
112 // public float idMaxDivisor(){
113 // return hits/(float)maxDivisor();
114 // }
115
116 public float idQueryDivisor(){
117 return hits/(float)(Tools.max(1, refDivisor));
118 }
119
120 public float idRefDivisor(){
121 return hits/(float)(Tools.max(1, refDivisor));
122 }
123
124 public float completeness(){
125 float complt=Tools.min(1, (queryDivisor-contamHits)/(float)Tools.max(1, refDivisor));
126 return complt;
127 // float c2=hits/(float)refDivisor;
128 // assert(queryDivisor-contamHits>=hits);
129 // assert(c1>=c2);
130 // System.err.println(hits+", "+contamHits+", "+refDivisor+", "+queryDivisor+", "+c1+", "+c2);
131 // return Tools.max(c1, c2);
132 // float kid=idMaxDivisor(), wkid=idMinDivisor();
133 // return kid==0 ? 0 : kid/wkid;
134 }
135
136 public float contamFraction(){
137 return Tools.min(1, contamHits/(float)Tools.max(1, queryDivisor));
138 }
139
140 public float contam2Fraction(){
141 return Tools.min(1, contam2Hits/(float)Tools.max(1, queryDivisor));
142 }
143
144 public float uContamFraction() {
145 int uContamHits=contamHits-multiContamHits;
146 return Tools.min(1, uContamHits/(float)Tools.max(1, queryDivisor));
147 }
148
149
150
151 /*--------------------------------------------------------------*/
152
153 final float wkid(){
154 final int div=minDivisor();
155 return hits/(float)div;
156 }
157 final float kid(){
158 final int div=maxDivisor();
159 return hits/(float)div;
160 }
161 final float aniOld(){
162 float wkid=wkid();
163 // final float ani=wkidToAni(wkid, k1Fraction());
164 final float ani=wkidToAni(wkid);
165 return ani;
166 }
167 final float ani(){
168 if(hits<1){return 0;}
169 final float ani;
170 if(k2>0 && useToValue2){
171 float ani1=ani1();
172 float ani2=ani2();
173 // ani=0.5f*(ani1+ani2);
174 ani=0.5f*(Tools.max(0.9f*ani2, ani1)+Tools.max(0.8f*ani1, ani2));
175 // return (ani1*qSeen1+ani2*qSeen2())/queryDivisor;
176
177 // System.err.println("ani="+ani+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)aniDual()+"\n"
178 //// +"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
179 // + "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
180 // + "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
181 }else{
182 ani=aniOld();
183 }
184 return ani;
185 }
186
187 final float wkid1(){
188 final int div=minDivisor1();
189 return hits1()/(float)div;
190 }
191 final float kid1(){
192 final int div=maxDivisor1();
193 return hits1()/(float)div;
194 }
195 final float ani1(){
196 float wkid=wkid1();
197 final float ani=wkidToAniExact(wkid, k);
198 return ani;
199 }
200
201 final float wkid2(){
202 final int div=minDivisor2();
203 return hits2()/(float)div;
204 }
205 final float kid2(){
206 final int div=maxDivisor2();
207 return hits2()/(float)div;
208 }
209 final float ani2(){
210 assert(k2>0);
211 float wkid=wkid2();
212 final float ani=wkidToAniExact(wkid, k2);
213 return ani;
214 }
215
216 final float aniDual(){
217 assert(k2>0);
218 float wkid1=wkid1();
219 float wkid2=wkid2();
220 float ratio=(wkid1/wkid2);
221 float exp=1f/(k-k2);//TODO - make this initialized
222 double ani=Math.pow(ratio, exp);
223 double gf=wkid2/Math.pow(ani, k2);
224
225 // System.err.println("ani="+ani()+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)ani+"\n"
226 // +"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
227 // + "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
228 // + "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
229
230 return (float)ani;
231 }
232
233 /*--------------------------------------------------------------*/
234
235 int hits1(){return hits1;}
236 int qSeen1(){return qSeen1;}
237 int rSeen1(){return rSeen1;}
238 int minDivisor1(){return Tools.max(1, Tools.min(qSeen1, rSeen1));}
239 int maxDivisor1(){return Tools.max(1, qSeen1, rSeen1);}
240
241 int hits2(){return hits-hits1;}
242 int qSeen2(){return queryDivisor-qSeen1;}
243 int rSeen2(){return refDivisor-rSeen1;}
244 int minDivisor2(){return Tools.max(1, Tools.min(qSeen2(), rSeen2()));}
245 int maxDivisor2(){return Tools.max(1, qSeen2(), rSeen2());}
246
247 /*--------------------------------------------------------------*/
248
249 // public float aniOld(){
250 // if(hits<1){return 0;}
251 //
252 //// double wkid=aniFromWKID ? idMinDivisor() : idMaxDivisor();
253 // double wkid=idMinDivisor();
254 // return wkidToAni(wkid, k1Fraction());
255 //
256 //// final float rID=hits/(float)(refDivisor);
257 //// final float qID=hits/(float)(queryDivisor-contamHits);
258 //// final float wkid2=Tools.max(qID, rID);
259 //// final float ani=wkidToAni(wkid2);
260 ////
261 ////// System.err.println("rid: "+wkidToAni(rID)+", qid: "+wkidToAni(qID)+", qid2: "+wkidToAni(hits/(float)(queryDivisor)));
262 ////
263 //// return ani;
264 // }
265
266 int minDivisor(){return Tools.max(1, Tools.min(refDivisor, queryDivisor));}
267 int maxDivisor(){return Tools.max(1, refDivisor, queryDivisor);}
268
269 private float score0_old(){
270 long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
271 float wkid=wkid();
272 float kid=kid();
273 float complt=completeness();
274 float contam=contamFraction();
275 float refHits=Tools.max(avgRefHits, 1f);
276 float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
277 return (float)((Math.log(hits+2)*.25f+0.5f)*(refHitMult*0.2*Math.log(hits+2+(1.2*uHits()+0.25*unique2()+0.1*unique3()))
278 *Math.sqrt(40*(20000+hits+uHits())*(wkid*kid*Math.pow(est*complt, 0.2)*(1-contam*0.1)))+0.1));
279 }
280
281 private float score0(){
282 final long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
283 final float wkid=wkid();
284 final float kid=kid();
285 final float ani=ani();
286 final float complt=completeness();
287 final float contam=contamFraction();
288 final float refHits=Tools.max(avgRefHits, 1f);
289 final float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
290 final float uhits=uHits();
291 final float uhits2=unique2();
292 final float uhits3=unique3();
293 final float contamMult=(1-contam*0.95f);
294 final float estMult=(float)(Math.pow(est, 0.2)*Math.sqrt(complt));
295 final float aniMult=(float)(ani*Math.sqrt(wkid*kid));
296
297 final float hitsSum=1+hits+uhits+0.5f*uhits2+0.25f*uhits3;
298
299 final float score=(float)(Math.log(Tools.max(1.2f, hits-1))*hitsSum*refHitMult*contamMult*aniMult*estMult);//+(Math.log(hitsSum+1)*wkid*complt));
300 return (float)(8*Math.sqrt(score));
301 }
302
303 private long range(){//TODO Make sure these are calculated correctly; it seems like one divisor might be 1 higher than necessary.
304 long maxA=a.keys[Tools.max(0, queryDivisor-1)];
305 long maxB=b.keys[Tools.max(0, refDivisor-1)];
306 // assert(false) : Tools.max(0, queryDivisor-1)+", "+Tools.max(0, refDivisor-1)+
307 // ", "+a.array[Tools.max(0, queryDivisor-1)]+", "+b.array[Tools.max(0, refDivisor-1)]+", "+Tools.max(maxA, maxB);//+"\n\n"+Arrays.toString(a.array)+"\n\n"+Arrays.toString(b.array);
308 return Tools.min(maxA, maxB);
309 }
310
311 private static double eValue(int hits, int minDiv, int maxDiv, long range){
312 if(hits>=range || maxDiv>=range){return 1.0;}
313 double probHit=maxDiv/(double)range;//Saturation of range
314 // double probNoHit=1-probHit;
315 double eValue=Math.pow(probHit, hits); //This is a simplification, assuming hits are very improbable.
316 //Note that this does not take into account minDiv, the number of attempts... but it should.
317 // System.err.println("hits="+hits+", minDiv="+minDiv+", maxDiv="+maxDiv+", range="+range+", eValue="+eValue);
318 return eValue;
319 }
320
321 public double eValue(){
322 double eValue=eValue1()*eValue2();
323 // System.err.println("eValue="+eValue);
324 return eValue;
325 }
326
327 public double eValue1(){
328 long range0=range();
329 int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k;
330 double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
331 int hits=hits1;
332 int minDiv=Tools.min(qSeen1, rSeen1);
333 int maxDiv=Tools.max(qSeen1, rSeen1);
334 long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
335 return eValue(hits, minDiv, maxDiv, range);
336 }
337
338 public double eValue2(){
339 if(k2<1){return 1.0;}
340
341 long range0=range();
342 int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k2;
343 double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
344 int hits=hits2();
345 int minDiv=Tools.min(qSeen2(), rSeen2());
346 int maxDiv=Tools.max(qSeen2(), rSeen2());
347 long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
348 // assert(false) : missingBits+", "+quantizer+", "+range0+", "+range+", "+eValue(hits, minDiv, maxDiv, range);
349 return eValue(hits, minDiv, maxDiv, range);
350 }
351
352 public String scoreS(){
353 float x=score;
354 return format3(x);
355 }
356
357 public double depth(boolean observedToActual){
358 return observedToActual ? Tools.observedToActualCoverage(depth) : depth;
359 }
360
361 public double depth2(boolean observedToActual){
362 return observedToActual ? Tools.observedToActualCoverage(depth2) : depth2;
363 }
364
365 public String depthS(boolean observedToActual){
366 float x=depth;
367 if(observedToActual){x=(float)Tools.observedToActualCoverage(x);}
368 return format3(x);
369 }
370
371 public float avgRefHits(){
372 return avgRefHits;
373 }
374
375 public String avgRefHitsS(){
376 return format2(avgRefHits);
377 }
378
379 public String depth2S(boolean observedToActual){
380 float x=depth2;
381 if(observedToActual){
382 x=(float)(Tools.observedToActualCoverage(depth)*(depth2/depth));
383 }
384 return format3(x);
385 }
386
387 public String volumeS(){
388 double x=volume()*0.001;
389 return format3(x);
390 }
391
392 static String format3(double x){
393 if(x>=999.95){
394 return(""+(long)Math.round(x));
395 }else if(x>=99.995){
396 return String.format(Locale.ROOT, "%.1f", x);
397 }else if(x>=9.9995){
398 return String.format(Locale.ROOT, "%.2f", x);
399 }
400 return String.format(Locale.ROOT, "%.3f", x);
401 }
402
403 static String format2(double x){
404 if(x>=999.95){
405 return(""+(long)Math.round(x));
406 }else if(x>=99.995){
407 return String.format(Locale.ROOT, "%.1f", x);
408 }else if(x>=9.9995){
409 return String.format(Locale.ROOT, "%.2f", x);
410 }
411 return String.format(Locale.ROOT, "%.2f", x);
412 }
413
414 float volume(){
415 return Tools.max(1f, depth)*hits;
416 }
417
418 @Override
419 public String toString(){
420 return "hits="+hits+", refDivisor="+refDivisor+", queryDivisor="+queryDivisor+", refSize="+refSize+", querySize="+querySize+
421 ", contamHits="+contamHits+", contam2Hits="+contam2Hits+", multiContamHits="+multiContamHits+", depth="+depth+", depth2="+depth2+", volume="+volume()+
422 ", hits="+hits+", multiHits="+multiHits+", unique2="+unique2+", unique3="+unique3+", noHits="+noHits+", taxID="+taxID+", taxName="+taxName;
423 }
424
425 /*--------------------------------------------------------------*/
426 /*---------------- Getters ----------------*/
427 /*--------------------------------------------------------------*/
428
429 // public boolean passesFilter(TaxFilter white, TaxFilter black){
430 // if(white==null && black==null){return true;}
431 // int id=taxID();
432 // String s=name();
433 // return passesFilter(white, id, s) && passesFilter(black, id, s);
434 // }
435 //
436 // private boolean passesFilter(TaxFilter filter, int id, String s){
437 // if(filter==null){return true;}
438 // if(id>0 && !filter.passesFilter(id)){return false;}
439 // if(s!=null && !filter.passesFilterByNameOnly(s)){return false;}
440 // return true;
441 // }
442
443 public String name(){return taxName!=null ? taxName : name0()!=null ? name0() : fname()!=null ? fname() : ""+taxID();}
444 public String taxName(){return taxName;}
445 String name0(){return b.name0();}
446 String fname(){return b.fname();}
447
448 // public int taxID(){return b.taxID<minFakeID ? b.taxID : 0;}
449 public int taxID(){return (taxID<minFakeID && taxID>=0) ? taxID : -1;}
450 long imgID(){return (b.imgID>0 ? b.imgID : -1);}
451
452 long genomeSizeBases(){return b.genomeSizeBases;}
453 long genomeSizeKmers(){return b.genomeSizeKmers;}
454 long genomeSequences(){return b.genomeSequences;}
455 long genomeSizeEstimate(){return b.genomeSizeEstimate();}
456 float gc(){return b.gc();}
457 boolean hasGC(){return b.baseCounts!=null;}
458
459 public boolean hasSSU() {
460 return (a.r16S()!=null && b.r16S()!=null) || (a.r18S()!=null && b.r18S()!=null);
461 }
462 public boolean hasSSUIdentity() {
463 return ssuIdentity>=0;
464 }
465 public boolean needsAlignment() {
466 return hasSSU() && !hasSSUIdentity();
467 }
468 public boolean hasQueryTaxID() {
469 return a.taxID>0 && a.taxID<minFakeID;
470 }
471
472
473 public int uHits() {return hits-multiHits;}
474
475 /** Common ancestor TaxID, if both Sketches have a TaxID */
476 public int commonAncestor() {
477 if(a.taxID<1 || b.taxID<1){return -1;}
478 assert(taxtree!=null);
479 int id=taxtree.commonAncestor(a.taxID, b.taxID);
480 return id;
481 }
482
483 /** Common ancestor node tax level, if both Sketches have a TaxID */
484 public String commonAncestorLevel() {
485 int id=commonAncestor();
486 if(id<1){return ".";}
487 TaxNode tn=taxtree.getNode(id);
488 while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
489 String s=tn.levelStringExtended(false);
490 return s;
491 }
492
493 /** Common ancestor node tax level, if both Sketches have a TaxID */
494 public int commonAncestorLevelInt() {
495 int id=commonAncestor();
496 if(id<1){return 0;}
497 TaxNode tn=taxtree.getNode(id);
498 while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
499 return tn.levelExtended;
500 }
501
502 /*--------------------------------------------------------------*/
503 /*---------------- Comparators ----------------*/
504 /*--------------------------------------------------------------*/
505
506
507
508 static class ScoreComparator implements Comparator<Comparison>{
509
510 @Override
511 public int compare(Comparison a, Comparison b) {
512 {
513 float pa=a.score, pb=b.score;
514 if(pa>pb){
515 return 1;
516 }else if (pa<pb){
517 return -1;
518 }
519 }
520
521 int x=a.hits-b.hits;
522 if(x!=0){return x;}
523 x=b.minDivisor()-a.minDivisor();
524 if(x!=0){return x;}
525 x=b.maxDivisor()-a.maxDivisor();
526 if(x!=0){return x;}
527 x=b.refDivisor-a.refDivisor;
528 if(x!=0){return x;}
529 x=a.taxID()-b.taxID();
530 if(x!=0){return x;}
531 if(a.name0()!=null && b.name0()!=null){
532 return a.name0().compareTo(b.name0());
533 }
534 if(a.taxName()!=null && b.taxName()!=null){
535 return a.taxName().compareTo(b.taxName());
536 }
537 return 0;
538 }
539
540 @Override
541 public String toString(){return "sortByScore";}
542
543 }
544
545 static class DepthComparator implements Comparator<Comparison>{
546
547 @Override
548 public int compare(Comparison a, Comparison b) {
549 final float da=Tools.max(0.1f, a.depth-0.5f), db=Tools.max(0.1f, b.depth-0.5f);
550 final float sa, sb;
551 if(sqrt){
552 sa=da*(float)Math.sqrt(a.score);
553 sb=db*(float)Math.sqrt(b.score);
554 }else{
555 sa=da*a.score;
556 sb=db*b.score;
557 }
558 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
559 }
560
561 @Override
562 public String toString(){return "sortByDepth";}
563
564 }
565
566 static class Depth2Comparator implements Comparator<Comparison>{
567
568 @Override
569 public int compare(Comparison a, Comparison b) {
570 final float da=Tools.max(0.1f, a.depth2-0.8f), db=Tools.max(0.1f, b.depth2-0.8f);
571 final float sa, sb;
572 if(sqrt){
573 sa=da*(float)Math.sqrt(a.score);
574 sb=db*(float)Math.sqrt(b.score);
575 }else{
576 sa=da*a.score;
577 sb=db*b.score;
578 }
579 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
580 }
581
582 @Override
583 public String toString(){return "sortByDepth2";}
584
585 }
586
587 static class VolumeComparator implements Comparator<Comparison>{
588
589 @Override
590 public int compare(Comparison a, Comparison b) {
591 final float da=a.volume(), db=b.volume();
592 final float sa, sb;
593 if(sqrt){
594 sa=da*(float)Math.sqrt(a.score);
595 sb=db*(float)Math.sqrt(b.score);
596 }else{
597 sa=da*a.score;
598 sb=db*b.score;
599 }
600 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
601 }
602
603 @Override
604 public String toString(){return "sortByVolume";}
605
606 }
607
608 static class KIDComparator implements Comparator<Comparison>{
609
610 @Override
611 public int compare(Comparison a, Comparison b) {
612 final float sa=a.kid(), sb=b.kid();
613 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
614 }
615
616 @Override
617 public String toString(){return "sortByKID";}
618
619 }
620
621 static class WKIDComparator implements Comparator<Comparison>{
622
623 @Override
624 public int compare(Comparison a, Comparison b) {
625 final float sa=a.wkid(), sb=b.wkid();
626 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
627 }
628
629 @Override
630 public String toString(){return "sortByWKID";}
631
632 }
633
634 static class SSUComparator implements Comparator<Comparison>{
635
636 @Override
637 public int compare(Comparison a, Comparison b) {
638 // if((a.has16S() || a.has18S()) && !a.hasSSUIdentity()){
639 // synchronized(a){a.ssuIdentity();}
640 // assert(a.hasSSUIdentity());
641 // }
642 // if((b.has16S() || b.has18S()) && !b.hasSSUIdentity()){
643 // synchronized(b){b.ssuIdentity();}
644 // assert(b.hasSSUIdentity());
645 // }
646
647 if(a.hasSSUIdentity() && b.hasSSUIdentity()){
648 final float sa=a.ssuIdentity(), sb=b.ssuIdentity();
649 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
650 }else if(a.hasSSUIdentity()){
651 return 1;
652 }else if(b.hasSSUIdentity()){
653 return -1;
654 }else{
655 return scoreComparator.compare(a, b);
656 }
657
658 }
659
660 @Override
661 public String toString(){return "sortBySSU";}
662
663 }
664
665 static class HitsComparator implements Comparator<Comparison>{
666
667 @Override
668 public int compare(Comparison a, Comparison b) {
669 final float sa=a.hits(), sb=b.hits();
670 return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
671 }
672
673 @Override
674 public String toString(){return "sortByHits";}
675
676 }
677
678 @Override
679 public int compareTo(Comparison b) {
680 assert(false) : "Please use comparators instead.";
681 return scoreComparator.compare(this, b);
682 }
683
684 @Override
685 public int hashCode() {
686 assert(false) : "TODO";
687 return super.hashCode();
688 }
689
690 public static final ScoreComparator scoreComparator=new ScoreComparator();
691 public static final DepthComparator depthComparator=new DepthComparator();
692 public static final Depth2Comparator depth2Comparator=new Depth2Comparator();
693 public static final VolumeComparator volumeComparator=new VolumeComparator();
694 public static final KIDComparator KIDComparator=new KIDComparator();
695 public static final SSUComparator SSUComparator=new SSUComparator();
696 public static final WKIDComparator WKIDComparator=new WKIDComparator();
697 public static final HitsComparator HitsComparator=new HitsComparator();
698 private static final boolean sqrt=false;
699 private static final double aaBitValue=0.86438561897747246957406388589788; //(log(20)/log(2))/5;
700
701 /*--------------------------------------------------------------*/
702
703 int hits(){return hits;}
704 int multiHits(){return multiHits;}
705 int noHits(){return noHits;}
706 int unique2(){return unique2;}
707 int unique3(){return unique3;}
708
709 float depth(){return depth;}
710 float depth2(){return depth2;}
711 float score(){return score;}
712
713 int contamHits(){return contamHits;}
714 int contam2Hits(){return contam2Hits;}
715 int multiContamHits(){return multiContamHits;}
716
717 int queryDivisor(){return queryDivisor;}
718 int refDivisor(){return refDivisor;}
719
720 int querySize(){return querySize;}
721 int refSize(){return refSize;}
722
723 int ssuType(){return has18S() ? 18 : has16S() ? 16 : 0;}
724 int ssuLen(){return has18S() ? b.r18SLen() : has16S() ? b.r16SLen() : 0;}
725 boolean has16S(){return a.r16S()!=null && b.r16S()!=null;}
726 boolean has18S(){return a.r18S()!=null && b.r18S()!=null;}
727
728 float ssuIdentity(){
729 if(ssuIdentity>0){return ssuIdentity;}
730 if(has18S()){ssuIdentity=calcIdentity(a.r18S(), b.r18S());}
731 else if(has16S()){ssuIdentity=calcIdentity(a.r16S(), b.r16S());}
732 return ssuIdentity;
733 }
734
735 private static float calcIdentity(byte[] ssuA, byte[] ssuB){
736 //assert(false);
737 if(ssuA.length>ssuB.length){
738 byte[] c=ssuA;
739 ssuA=ssuB;
740 ssuB=c;
741 }
742 if(useSSA){
743 Aligner ssa=(useSSA3 ? GeneCaller.getSSA3() : GeneCaller.getSSA());
744 int[] max=ssa.fillUnlimited(ssuA, ssuB, 0, ssuB.length-1, 0);
745 if(max==null){return 0;}
746
747 final int rows=max[0];
748 final int maxCol=max[1];
749 final int maxState=max[2];
750
751 //returns {score, bestRefStart, bestRefStop}
752 //padded: {score, bestRefStart, bestRefStop, padLeft, padRight};
753 int[] score=ssa.score(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
754 int rstart=score[1];
755 int rstop=score[2];
756
757 // byte[] match=ssa.traceback(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
758 // float id=Read.identity(match);
759 float id=ssa.tracebackIdentity(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState, null);
760 return id;
761 }else{
762 GlocalAligner ga=new GlocalAligner();
763 return ga.alignForward(ssuA, ssuB);
764 }
765 }
766
767 /*--------------------------------------------------------------*/
768 /*---------------- Fields ----------------*/
769 /*--------------------------------------------------------------*/
770
771 public Sketch a, b;
772
773 String taxName;
774 int taxID;
775
776 private int hits;
777 private int multiHits;
778 private int unique2;
779 private int unique3;
780 private int noHits;
781
782 private float depth;
783 private float depth2;
784 private float avgRefHits;
785 private float score;
786
787 private int contamHits;
788 private int contam2Hits;
789 private int multiContamHits;
790
791 private int refDivisor;
792 private int queryDivisor;
793
794 private int refSize;
795 private int querySize;
796
797 private int hits1;
798 private int qSeen1;
799 private int rSeen1;
800
801 private float ssuIdentity=-1;
802
803 }