jpayne@68: package sketch;
jpayne@68: 
jpayne@68: import java.util.Comparator;
jpayne@68: import java.util.Locale;
jpayne@68: 
jpayne@68: import aligner.Aligner;
jpayne@68: import prok.GeneCaller;
jpayne@68: import shared.Tools;
jpayne@68: import tax.TaxNode;
jpayne@68: 
jpayne@68: public final class Comparison extends SketchObject implements Comparable<Comparison> {
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Constructors         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: //	public Comparison(CompareBuffer buffer){
jpayne@68: //		this(buffer, null, null);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: //	public Comparison(Sketch a_, Sketch b_){
jpayne@68: //		this(null, a_, b_);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	public Comparison(CompareBuffer buffer, Sketch a_, Sketch b_){
jpayne@68: 		
jpayne@68: 		a=a_;
jpayne@68: 		b=b_;
jpayne@68: 		
jpayne@68: 		if(buffer!=null){setFrom(buffer);}
jpayne@68: 		
jpayne@68: 		if(b!=null){
jpayne@68: 			taxName=b.taxName();
jpayne@68: 			taxID=b.taxID;
jpayne@68: 		}
jpayne@68: 
jpayne@68: //		System.err.println(this);
jpayne@68: //		System.err.println(b.present);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------           Mutators           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public void setFrom(CompareBuffer buffer){
jpayne@68: 		hits=buffer.hits();
jpayne@68: 		multiHits=buffer.multiHits();
jpayne@68: 		unique2=buffer.unique2();
jpayne@68: 		unique3=buffer.unique3();
jpayne@68: 		noHits=buffer.noHits();
jpayne@68: 
jpayne@68: 		contamHits=buffer.contamHits();
jpayne@68: 		contam2Hits=buffer.contam2Hits();
jpayne@68: 		multiContamHits=buffer.multiContamHits();
jpayne@68: 		
jpayne@68: 		refDivisor=buffer.refDivisor();
jpayne@68: 		queryDivisor=buffer.queryDivisor();
jpayne@68: 		
jpayne@68: 		refSize=buffer.refSize();
jpayne@68: 		querySize=buffer.querySize();
jpayne@68: 
jpayne@68: 		depth=buffer.depth();
jpayne@68: 		depth2=buffer.depth2();
jpayne@68: 		float x=buffer.avgRefHits();
jpayne@68: 		if(x>0){avgRefHits=x;}
jpayne@68: //		volume=volume0();
jpayne@68: 		score=score0();
jpayne@68: 
jpayne@68: 		hits1=buffer.hits1();
jpayne@68: 		qSeen1=buffer.qSeen1();
jpayne@68: 		rSeen1=buffer.rSeen1();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void recompare(CompareBuffer buffer, int[][] taxHits, int contamLevel){
jpayne@68: 		
jpayne@68: //		for(int[] row : taxHits){
jpayne@68: //			if(row!=null){
jpayne@68: //				System.err.println(Arrays.toString(row));
jpayne@68: //			}
jpayne@68: //		}//Tested; correctly indicates most rows have octopus but some have other things.
jpayne@68: 		
jpayne@68: 		assert(a.merged());
jpayne@68: //		int oldContam2=contam2Hits;
jpayne@68: 		int x=a.countMatches(b, buffer, a.compareBitSet(), false, taxHits, contamLevel);
jpayne@68: 		assert(x==hits);
jpayne@68: 		setFrom(buffer);
jpayne@68: //		contam2Hits=oldContam2;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Methods           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public boolean equals(Object b){
jpayne@68: 		if(b==null || b.getClass()!=this.getClass()){return false;}
jpayne@68: 		return scoreComparator.compare(this, (Comparison)b)==0;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: //	//WKID
jpayne@68: //	public float wkid(){return idMinDivisor();}
jpayne@68: //	public float idMinDivisor(){
jpayne@68: //		return hits/(float)minDivisor();
jpayne@68: //	}
jpayne@68: 	
jpayne@68: //	public float k1Fraction(){
jpayne@68: //		return a.k1Fraction();
jpayne@68: //	}
jpayne@68: 	
jpayne@68: //	//KID
jpayne@68: //	public float kid(){return idMaxDivisor();}
jpayne@68: //	public float idMaxDivisor(){
jpayne@68: //		return hits/(float)maxDivisor();
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	public float idQueryDivisor(){
jpayne@68: 		return hits/(float)(Tools.max(1, refDivisor));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float idRefDivisor(){
jpayne@68: 		return hits/(float)(Tools.max(1, refDivisor));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float completeness(){
jpayne@68: 		float complt=Tools.min(1, (queryDivisor-contamHits)/(float)Tools.max(1, refDivisor));
jpayne@68: 		return complt;
jpayne@68: //		float c2=hits/(float)refDivisor;
jpayne@68: //		assert(queryDivisor-contamHits>=hits);
jpayne@68: //		assert(c1>=c2);
jpayne@68: //		System.err.println(hits+", "+contamHits+", "+refDivisor+", "+queryDivisor+", "+c1+", "+c2);
jpayne@68: //		return Tools.max(c1, c2);
jpayne@68: //		float kid=idMaxDivisor(), wkid=idMinDivisor();
jpayne@68: //		return kid==0 ? 0 : kid/wkid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float contamFraction(){
jpayne@68: 		return Tools.min(1, contamHits/(float)Tools.max(1, queryDivisor));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float contam2Fraction(){
jpayne@68: 		return Tools.min(1, contam2Hits/(float)Tools.max(1, queryDivisor));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float uContamFraction() {
jpayne@68: 		int uContamHits=contamHits-multiContamHits;
jpayne@68: 		return Tools.min(1, uContamHits/(float)Tools.max(1, queryDivisor));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	final float wkid(){
jpayne@68: 		final int div=minDivisor();
jpayne@68: 		return hits/(float)div;
jpayne@68: 	}
jpayne@68: 	final float kid(){
jpayne@68: 		final int div=maxDivisor();
jpayne@68: 		return hits/(float)div;
jpayne@68: 	}
jpayne@68: 	final float aniOld(){
jpayne@68: 		float wkid=wkid();
jpayne@68: //		final float ani=wkidToAni(wkid, k1Fraction());
jpayne@68: 		final float ani=wkidToAni(wkid);
jpayne@68: 		return ani;
jpayne@68: 	}
jpayne@68: 	final float ani(){
jpayne@68: 		if(hits<1){return 0;}
jpayne@68: 		final float ani;
jpayne@68: 		if(k2>0 && useToValue2){
jpayne@68: 			float ani1=ani1();
jpayne@68: 			float ani2=ani2();
jpayne@68: //			ani=0.5f*(ani1+ani2);
jpayne@68: 			ani=0.5f*(Tools.max(0.9f*ani2, ani1)+Tools.max(0.8f*ani1, ani2));
jpayne@68: //			return (ani1*qSeen1+ani2*qSeen2())/queryDivisor;
jpayne@68: 			
jpayne@68: //			System.err.println("ani="+ani+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)aniDual()+"\n"
jpayne@68: ////					+"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
jpayne@68: //							+ "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
jpayne@68: //									+ "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
jpayne@68: 		}else{
jpayne@68: 			ani=aniOld();
jpayne@68: 		}
jpayne@68: 		return ani;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	final float wkid1(){
jpayne@68: 		final int div=minDivisor1();
jpayne@68: 		return hits1()/(float)div;
jpayne@68: 	}
jpayne@68: 	final float kid1(){
jpayne@68: 		final int div=maxDivisor1();
jpayne@68: 		return hits1()/(float)div;
jpayne@68: 	}
jpayne@68: 	final float ani1(){
jpayne@68: 		float wkid=wkid1();
jpayne@68: 		final float ani=wkidToAniExact(wkid, k);
jpayne@68: 		return ani;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	final float wkid2(){
jpayne@68: 		final int div=minDivisor2();
jpayne@68: 		return hits2()/(float)div;
jpayne@68: 	}
jpayne@68: 	final float kid2(){
jpayne@68: 		final int div=maxDivisor2();
jpayne@68: 		return hits2()/(float)div;
jpayne@68: 	}
jpayne@68: 	final float ani2(){
jpayne@68: 		assert(k2>0);
jpayne@68: 		float wkid=wkid2();
jpayne@68: 		final float ani=wkidToAniExact(wkid, k2);
jpayne@68: 		return ani;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	final float aniDual(){
jpayne@68: 		assert(k2>0);
jpayne@68: 		float wkid1=wkid1();
jpayne@68: 		float wkid2=wkid2();
jpayne@68: 		float ratio=(wkid1/wkid2);
jpayne@68: 		float exp=1f/(k-k2);//TODO - make this initialized
jpayne@68: 		double ani=Math.pow(ratio, exp);
jpayne@68: 		double gf=wkid2/Math.pow(ani, k2);
jpayne@68: 		
jpayne@68: //		System.err.println("ani="+ani()+"aniOld="+aniOld()+", ani1="+ani1()+", ani2="+ani2()+", anid="+(float)ani+"\n"
jpayne@68: //				+"gf="+(float)gf+", wkid1="+wkid1+", wkid2="+wkid2+"\n"
jpayne@68: //						+ "k1f="+k1Fraction()+", hits="+hits+", hits1="+hits1+", hits2="+hits2()+", qSeen1()="+qSeen1()+", rSeen1()="+rSeen1()+"\n"
jpayne@68: //								+ "qSeen2()="+qSeen2()+", rSeen2()="+rSeen2()+", minDivisor1()="+minDivisor1()+", minDivisor2()="+minDivisor2()+"\n");
jpayne@68: 		
jpayne@68: 		return (float)ani;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	int hits1(){return hits1;}
jpayne@68: 	int qSeen1(){return qSeen1;}
jpayne@68: 	int rSeen1(){return rSeen1;}
jpayne@68: 	int minDivisor1(){return Tools.max(1, Tools.min(qSeen1, rSeen1));}
jpayne@68: 	int maxDivisor1(){return Tools.max(1, qSeen1, rSeen1);}
jpayne@68: 
jpayne@68: 	int hits2(){return hits-hits1;}
jpayne@68: 	int qSeen2(){return queryDivisor-qSeen1;}
jpayne@68: 	int rSeen2(){return refDivisor-rSeen1;}
jpayne@68: 	int minDivisor2(){return Tools.max(1, Tools.min(qSeen2(), rSeen2()));}
jpayne@68: 	int maxDivisor2(){return Tools.max(1, qSeen2(), rSeen2());}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: //	public float aniOld(){
jpayne@68: //		if(hits<1){return 0;}
jpayne@68: //
jpayne@68: ////		double wkid=aniFromWKID ? idMinDivisor() : idMaxDivisor();
jpayne@68: //		double wkid=idMinDivisor();
jpayne@68: //		return wkidToAni(wkid, k1Fraction());
jpayne@68: //
jpayne@68: ////		final float rID=hits/(float)(refDivisor);
jpayne@68: ////		final float qID=hits/(float)(queryDivisor-contamHits);
jpayne@68: ////		final float wkid2=Tools.max(qID, rID);
jpayne@68: ////		final float ani=wkidToAni(wkid2);
jpayne@68: ////
jpayne@68: //////		System.err.println("rid: "+wkidToAni(rID)+", qid: "+wkidToAni(qID)+", qid2: "+wkidToAni(hits/(float)(queryDivisor)));
jpayne@68: ////
jpayne@68: ////		return ani;
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	int minDivisor(){return Tools.max(1, Tools.min(refDivisor, queryDivisor));}
jpayne@68: 	int maxDivisor(){return Tools.max(1, refDivisor, queryDivisor);}
jpayne@68: 	
jpayne@68: 	private float score0_old(){
jpayne@68: 		long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
jpayne@68: 		float wkid=wkid();
jpayne@68: 		float kid=kid();
jpayne@68: 		float complt=completeness();
jpayne@68: 		float contam=contamFraction();
jpayne@68: 		float refHits=Tools.max(avgRefHits, 1f);
jpayne@68: 		float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
jpayne@68: 		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()))
jpayne@68: 				*Math.sqrt(40*(20000+hits+uHits())*(wkid*kid*Math.pow(est*complt, 0.2)*(1-contam*0.1)))+0.1));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private float score0(){
jpayne@68: 		final long est=useSizeEstimate ? genomeSizeEstimate() : genomeSizeKmers();
jpayne@68: 		final float wkid=wkid();
jpayne@68: 		final float kid=kid();
jpayne@68: 		final float ani=ani();
jpayne@68: 		final float complt=completeness();
jpayne@68: 		final float contam=contamFraction();
jpayne@68: 		final float refHits=Tools.max(avgRefHits, 1f);
jpayne@68: 		final float refHitMult=1f+(0.6f/(float)Math.sqrt(refHits+1));
jpayne@68: 		final float uhits=uHits();
jpayne@68: 		final float uhits2=unique2();
jpayne@68: 		final float uhits3=unique3();
jpayne@68: 		final float contamMult=(1-contam*0.95f);
jpayne@68: 		final float estMult=(float)(Math.pow(est, 0.2)*Math.sqrt(complt));
jpayne@68: 		final float aniMult=(float)(ani*Math.sqrt(wkid*kid));
jpayne@68: 		
jpayne@68: 		final float hitsSum=1+hits+uhits+0.5f*uhits2+0.25f*uhits3;
jpayne@68: 		
jpayne@68: 		final float score=(float)(Math.log(Tools.max(1.2f, hits-1))*hitsSum*refHitMult*contamMult*aniMult*estMult);//+(Math.log(hitsSum+1)*wkid*complt));
jpayne@68: 		return (float)(8*Math.sqrt(score));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private long range(){//TODO Make sure these are calculated correctly; it seems like one divisor might be 1 higher than necessary.
jpayne@68: 		long maxA=a.keys[Tools.max(0, queryDivisor-1)];
jpayne@68: 		long maxB=b.keys[Tools.max(0, refDivisor-1)];
jpayne@68: //		assert(false) : Tools.max(0, queryDivisor-1)+", "+Tools.max(0, refDivisor-1)+
jpayne@68: //			", "+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);
jpayne@68: 		return Tools.min(maxA, maxB);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private static double eValue(int hits, int minDiv, int maxDiv, long range){
jpayne@68: 		if(hits>=range || maxDiv>=range){return 1.0;}
jpayne@68: 		double probHit=maxDiv/(double)range;//Saturation of range
jpayne@68: //		double probNoHit=1-probHit;
jpayne@68: 		double eValue=Math.pow(probHit, hits);  //This is a simplification, assuming hits are very improbable.
jpayne@68: 		//Note that this does not take into account minDiv, the number of attempts...  but it should.
jpayne@68: //		System.err.println("hits="+hits+", minDiv="+minDiv+", maxDiv="+maxDiv+", range="+range+", eValue="+eValue);
jpayne@68: 		return eValue;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public double eValue(){
jpayne@68: 		double eValue=eValue1()*eValue2();
jpayne@68: //		System.err.println("eValue="+eValue);
jpayne@68: 		return eValue;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public double eValue1(){
jpayne@68: 		long range0=range();
jpayne@68: 		int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k;
jpayne@68: 		double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
jpayne@68: 		int hits=hits1;
jpayne@68: 		int minDiv=Tools.min(qSeen1, rSeen1);
jpayne@68: 		int maxDiv=Tools.max(qSeen1, rSeen1);
jpayne@68: 		long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
jpayne@68: 		return eValue(hits, minDiv, maxDiv, range);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public double eValue2(){
jpayne@68: 		if(k2<1){return 1.0;}
jpayne@68: 		
jpayne@68: 		long range0=range();
jpayne@68: 		int missingBits=64-(aminoOrTranslate() ? 5 : 2)*k2;
jpayne@68: 		double quantizer=1.0/(aminoOrTranslate() ? Math.pow(2, missingBits*aaBitValue) : 1L<<missingBits);
jpayne@68: 		int hits=hits2();
jpayne@68: 		int minDiv=Tools.min(qSeen2(), rSeen2());
jpayne@68: 		int maxDiv=Tools.max(qSeen2(), rSeen2());
jpayne@68: 		long range=Tools.max((long)Math.ceil(range0*quantizer), maxDiv);
jpayne@68: //		assert(false) : missingBits+", "+quantizer+", "+range0+", "+range+", "+eValue(hits, minDiv, maxDiv, range);
jpayne@68: 		return eValue(hits, minDiv, maxDiv, range);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public String scoreS(){
jpayne@68: 		float x=score;
jpayne@68: 		return format3(x);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public double depth(boolean observedToActual){
jpayne@68: 		return observedToActual ? Tools.observedToActualCoverage(depth) : depth;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public double depth2(boolean observedToActual){
jpayne@68: 		return observedToActual ? Tools.observedToActualCoverage(depth2) : depth2;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public String depthS(boolean observedToActual){
jpayne@68: 		float x=depth;
jpayne@68: 		if(observedToActual){x=(float)Tools.observedToActualCoverage(x);}
jpayne@68: 		return format3(x);
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public float avgRefHits(){
jpayne@68: 		return avgRefHits;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public String avgRefHitsS(){
jpayne@68: 		return format2(avgRefHits);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public String depth2S(boolean observedToActual){
jpayne@68: 		float x=depth2;
jpayne@68: 		if(observedToActual){
jpayne@68: 			x=(float)(Tools.observedToActualCoverage(depth)*(depth2/depth));
jpayne@68: 		}
jpayne@68: 		return format3(x);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public String volumeS(){
jpayne@68: 		double x=volume()*0.001;
jpayne@68: 		return format3(x);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static String format3(double x){
jpayne@68: 		if(x>=999.95){
jpayne@68: 			return(""+(long)Math.round(x));
jpayne@68: 		}else if(x>=99.995){
jpayne@68: 			return String.format(Locale.ROOT, "%.1f", x);
jpayne@68: 		}else if(x>=9.9995){
jpayne@68: 			return String.format(Locale.ROOT, "%.2f", x);
jpayne@68: 		}
jpayne@68: 		return String.format(Locale.ROOT, "%.3f", x);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static String format2(double x){
jpayne@68: 		if(x>=999.95){
jpayne@68: 			return(""+(long)Math.round(x));
jpayne@68: 		}else if(x>=99.995){
jpayne@68: 			return String.format(Locale.ROOT, "%.1f", x);
jpayne@68: 		}else if(x>=9.9995){
jpayne@68: 			return String.format(Locale.ROOT, "%.2f", x);
jpayne@68: 		}
jpayne@68: 		return String.format(Locale.ROOT, "%.2f", x);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	float volume(){
jpayne@68: 		return Tools.max(1f, depth)*hits;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public String toString(){
jpayne@68: 		return "hits="+hits+", refDivisor="+refDivisor+", queryDivisor="+queryDivisor+", refSize="+refSize+", querySize="+querySize+
jpayne@68: 				", contamHits="+contamHits+", contam2Hits="+contam2Hits+", multiContamHits="+multiContamHits+", depth="+depth+", depth2="+depth2+", volume="+volume()+
jpayne@68: 				", hits="+hits+", multiHits="+multiHits+", unique2="+unique2+", unique3="+unique3+", noHits="+noHits+", taxID="+taxID+", taxName="+taxName;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Getters           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: //	public boolean passesFilter(TaxFilter white, TaxFilter black){
jpayne@68: //		if(white==null && black==null){return true;}
jpayne@68: //		int id=taxID();
jpayne@68: //		String s=name();
jpayne@68: //		return passesFilter(white, id, s) && passesFilter(black, id, s);
jpayne@68: //	}
jpayne@68: //	
jpayne@68: //	private boolean passesFilter(TaxFilter filter, int id, String s){
jpayne@68: //		if(filter==null){return true;}
jpayne@68: //		if(id>0 && !filter.passesFilter(id)){return false;}
jpayne@68: //		if(s!=null && !filter.passesFilterByNameOnly(s)){return false;}
jpayne@68: //		return true;
jpayne@68: //	}
jpayne@68: 
jpayne@68: 	public String name(){return taxName!=null ? taxName : name0()!=null ? name0() : fname()!=null ? fname() : ""+taxID();}
jpayne@68: 	public String taxName(){return taxName;}
jpayne@68: 	String name0(){return b.name0();}
jpayne@68: 	String fname(){return b.fname();}
jpayne@68: 
jpayne@68: //	public int taxID(){return b.taxID<minFakeID ? b.taxID : 0;}
jpayne@68: 	public int taxID(){return (taxID<minFakeID && taxID>=0) ? taxID : -1;}
jpayne@68: 	long imgID(){return (b.imgID>0 ? b.imgID : -1);}
jpayne@68: 	
jpayne@68: 	long genomeSizeBases(){return b.genomeSizeBases;}
jpayne@68: 	long genomeSizeKmers(){return b.genomeSizeKmers;}
jpayne@68: 	long genomeSequences(){return b.genomeSequences;}
jpayne@68: 	long genomeSizeEstimate(){return b.genomeSizeEstimate();}
jpayne@68: 	float gc(){return b.gc();}
jpayne@68: 	boolean hasGC(){return b.baseCounts!=null;}
jpayne@68: 
jpayne@68: 	public boolean hasSSU() {
jpayne@68: 		return (a.r16S()!=null && b.r16S()!=null) || (a.r18S()!=null && b.r18S()!=null);
jpayne@68: 	}
jpayne@68: 	public boolean hasSSUIdentity() {
jpayne@68: 		return ssuIdentity>=0;
jpayne@68: 	}
jpayne@68: 	public boolean needsAlignment() {
jpayne@68: 		return hasSSU() && !hasSSUIdentity();
jpayne@68: 	}
jpayne@68: 	public boolean hasQueryTaxID() {
jpayne@68: 		return a.taxID>0 && a.taxID<minFakeID;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 
jpayne@68: 	public int uHits() {return hits-multiHits;}
jpayne@68: 
jpayne@68: 	/** Common ancestor TaxID, if both Sketches have a TaxID */
jpayne@68: 	public int commonAncestor() {
jpayne@68: 		if(a.taxID<1 || b.taxID<1){return -1;}
jpayne@68: 		assert(taxtree!=null);
jpayne@68: 		int id=taxtree.commonAncestor(a.taxID, b.taxID);
jpayne@68: 		return id;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/** Common ancestor node tax level, if both Sketches have a TaxID */
jpayne@68: 	public String commonAncestorLevel() {
jpayne@68: 		int id=commonAncestor();
jpayne@68: 		if(id<1){return ".";}
jpayne@68: 		TaxNode tn=taxtree.getNode(id);
jpayne@68: 		while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
jpayne@68: 		String s=tn.levelStringExtended(false);
jpayne@68: 		return s;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/** Common ancestor node tax level, if both Sketches have a TaxID */
jpayne@68: 	public int commonAncestorLevelInt() {
jpayne@68: 		int id=commonAncestor();
jpayne@68: 		if(id<1){return 0;}
jpayne@68: 		TaxNode tn=taxtree.getNode(id);
jpayne@68: 		while(!tn.isRanked() && tn.pid!=tn.id){tn=taxtree.getNode(tn.pid);}
jpayne@68: 		return tn.levelExtended;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Comparators          ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	static class ScoreComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			{
jpayne@68: 				float pa=a.score, pb=b.score;
jpayne@68: 				if(pa>pb){
jpayne@68: 					return 1;
jpayne@68: 				}else if (pa<pb){
jpayne@68: 					return -1;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			int x=a.hits-b.hits;
jpayne@68: 			if(x!=0){return x;}
jpayne@68: 			x=b.minDivisor()-a.minDivisor();
jpayne@68: 			if(x!=0){return x;}
jpayne@68: 			x=b.maxDivisor()-a.maxDivisor();
jpayne@68: 			if(x!=0){return x;}
jpayne@68: 			x=b.refDivisor-a.refDivisor;
jpayne@68: 			if(x!=0){return x;}
jpayne@68: 			x=a.taxID()-b.taxID();
jpayne@68: 			if(x!=0){return x;}
jpayne@68: 			if(a.name0()!=null && b.name0()!=null){
jpayne@68: 				return a.name0().compareTo(b.name0());
jpayne@68: 			}
jpayne@68: 			if(a.taxName()!=null && b.taxName()!=null){
jpayne@68: 				return a.taxName().compareTo(b.taxName());
jpayne@68: 			}
jpayne@68: 			return 0;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByScore";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class DepthComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float da=Tools.max(0.1f, a.depth-0.5f), db=Tools.max(0.1f, b.depth-0.5f);
jpayne@68: 			final float sa, sb;
jpayne@68: 			if(sqrt){
jpayne@68: 				sa=da*(float)Math.sqrt(a.score);
jpayne@68: 				sb=db*(float)Math.sqrt(b.score);
jpayne@68: 			}else{
jpayne@68: 				sa=da*a.score;
jpayne@68: 				sb=db*b.score;
jpayne@68: 			}
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByDepth";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class Depth2Comparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float da=Tools.max(0.1f, a.depth2-0.8f), db=Tools.max(0.1f, b.depth2-0.8f);
jpayne@68: 			final float sa, sb;
jpayne@68: 			if(sqrt){
jpayne@68: 				sa=da*(float)Math.sqrt(a.score);
jpayne@68: 				sb=db*(float)Math.sqrt(b.score);
jpayne@68: 			}else{
jpayne@68: 				sa=da*a.score;
jpayne@68: 				sb=db*b.score;
jpayne@68: 			}
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByDepth2";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class VolumeComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float da=a.volume(), db=b.volume();
jpayne@68: 			final float sa, sb;
jpayne@68: 			if(sqrt){
jpayne@68: 				sa=da*(float)Math.sqrt(a.score);
jpayne@68: 				sb=db*(float)Math.sqrt(b.score);
jpayne@68: 			}else{
jpayne@68: 				sa=da*a.score;
jpayne@68: 				sb=db*b.score;
jpayne@68: 			}
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByVolume";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class KIDComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float sa=a.kid(), sb=b.kid();
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByKID";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class WKIDComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float sa=a.wkid(), sb=b.wkid();
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByWKID";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class SSUComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: //			if((a.has16S() || a.has18S()) && !a.hasSSUIdentity()){
jpayne@68: //				synchronized(a){a.ssuIdentity();}
jpayne@68: //				assert(a.hasSSUIdentity());
jpayne@68: //			}
jpayne@68: //			if((b.has16S() || b.has18S()) && !b.hasSSUIdentity()){
jpayne@68: //				synchronized(b){b.ssuIdentity();}
jpayne@68: //				assert(b.hasSSUIdentity());
jpayne@68: //			}
jpayne@68: 			
jpayne@68: 			if(a.hasSSUIdentity() && b.hasSSUIdentity()){
jpayne@68: 				final float sa=a.ssuIdentity(), sb=b.ssuIdentity();
jpayne@68: 				return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 			}else if(a.hasSSUIdentity()){
jpayne@68: 				return 1;
jpayne@68: 			}else if(b.hasSSUIdentity()){
jpayne@68: 				return -1;
jpayne@68: 			}else{
jpayne@68: 				return scoreComparator.compare(a, b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortBySSU";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static class HitsComparator implements Comparator<Comparison>{
jpayne@68: 
jpayne@68: 		@Override
jpayne@68: 		public int compare(Comparison a, Comparison b) {
jpayne@68: 			final float sa=a.hits(), sb=b.hits();
jpayne@68: 			return sa>sb ? 1 : sa<sb ? -1 : scoreComparator.compare(a, b);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public String toString(){return "sortByHits";}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public int compareTo(Comparison b) {
jpayne@68: 		assert(false) : "Please use comparators instead.";
jpayne@68: 		return scoreComparator.compare(this, b);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public int hashCode() {
jpayne@68: 		assert(false) : "TODO";
jpayne@68: 		return super.hashCode();
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public static final ScoreComparator scoreComparator=new ScoreComparator();
jpayne@68: 	public static final DepthComparator depthComparator=new DepthComparator();
jpayne@68: 	public static final Depth2Comparator depth2Comparator=new Depth2Comparator();
jpayne@68: 	public static final VolumeComparator volumeComparator=new VolumeComparator();
jpayne@68: 	public static final KIDComparator KIDComparator=new KIDComparator();
jpayne@68: 	public static final SSUComparator SSUComparator=new SSUComparator();
jpayne@68: 	public static final WKIDComparator WKIDComparator=new WKIDComparator();
jpayne@68: 	public static final HitsComparator HitsComparator=new HitsComparator();
jpayne@68: 	private static final boolean sqrt=false;
jpayne@68: 	private static final double aaBitValue=0.86438561897747246957406388589788; //(log(20)/log(2))/5;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	int hits(){return hits;}
jpayne@68: 	int multiHits(){return multiHits;}
jpayne@68: 	int noHits(){return noHits;}
jpayne@68: 	int unique2(){return unique2;}
jpayne@68: 	int unique3(){return unique3;}
jpayne@68: 
jpayne@68: 	float depth(){return depth;}
jpayne@68: 	float depth2(){return depth2;}
jpayne@68: 	float score(){return score;}
jpayne@68: 
jpayne@68: 	int contamHits(){return contamHits;}
jpayne@68: 	int contam2Hits(){return contam2Hits;}
jpayne@68: 	int multiContamHits(){return multiContamHits;}
jpayne@68: 	
jpayne@68: 	int queryDivisor(){return queryDivisor;}
jpayne@68: 	int refDivisor(){return refDivisor;}
jpayne@68: 	
jpayne@68: 	int querySize(){return querySize;}
jpayne@68: 	int refSize(){return refSize;}
jpayne@68: 
jpayne@68: 	int ssuType(){return has18S() ? 18 : has16S() ? 16 : 0;}
jpayne@68: 	int ssuLen(){return has18S() ? b.r18SLen() : has16S() ? b.r16SLen() : 0;}
jpayne@68: 	boolean has16S(){return a.r16S()!=null && b.r16S()!=null;}
jpayne@68: 	boolean has18S(){return a.r18S()!=null && b.r18S()!=null;}
jpayne@68: 	
jpayne@68: 	float ssuIdentity(){
jpayne@68: 		if(ssuIdentity>0){return ssuIdentity;}
jpayne@68: 		if(has18S()){ssuIdentity=calcIdentity(a.r18S(), b.r18S());}
jpayne@68: 		else if(has16S()){ssuIdentity=calcIdentity(a.r16S(), b.r16S());}
jpayne@68: 		return ssuIdentity;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private static float calcIdentity(byte[] ssuA, byte[] ssuB){
jpayne@68: 		//assert(false);
jpayne@68: 		if(ssuA.length>ssuB.length){
jpayne@68: 			byte[] c=ssuA;
jpayne@68: 			ssuA=ssuB;
jpayne@68: 			ssuB=c;
jpayne@68: 		}
jpayne@68: 		if(useSSA){
jpayne@68: 			Aligner ssa=(useSSA3 ? GeneCaller.getSSA3() : GeneCaller.getSSA());
jpayne@68: 			int[] max=ssa.fillUnlimited(ssuA, ssuB, 0, ssuB.length-1, 0);
jpayne@68: 			if(max==null){return 0;}
jpayne@68: 			
jpayne@68: 			final int rows=max[0];
jpayne@68: 			final int maxCol=max[1];
jpayne@68: 			final int maxState=max[2];
jpayne@68: 			
jpayne@68: 			//returns {score, bestRefStart, bestRefStop} 
jpayne@68: 			//padded: {score, bestRefStart, bestRefStop, padLeft, padRight};
jpayne@68: 			int[] score=ssa.score(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
jpayne@68: 			int rstart=score[1];
jpayne@68: 			int rstop=score[2];
jpayne@68: 			
jpayne@68: //			byte[] match=ssa.traceback(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState);
jpayne@68: //			float id=Read.identity(match);
jpayne@68: 			float id=ssa.tracebackIdentity(ssuA, ssuB, 0, ssuB.length-1, rows, maxCol, maxState, null);
jpayne@68: 			return id;
jpayne@68: 		}else{
jpayne@68: 			GlocalAligner ga=new GlocalAligner();
jpayne@68: 			return ga.alignForward(ssuA, ssuB);
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Fields            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public Sketch a, b;
jpayne@68: 
jpayne@68: 	String taxName;
jpayne@68: 	int taxID;
jpayne@68: 	
jpayne@68: 	private int hits;
jpayne@68: 	private int multiHits;
jpayne@68: 	private int unique2;
jpayne@68: 	private int unique3;
jpayne@68: 	private int noHits;
jpayne@68: 
jpayne@68: 	private float depth;
jpayne@68: 	private float depth2;
jpayne@68: 	private float avgRefHits;
jpayne@68: 	private float score;
jpayne@68: 
jpayne@68: 	private int contamHits;
jpayne@68: 	private int contam2Hits;
jpayne@68: 	private int multiContamHits;
jpayne@68: 	
jpayne@68: 	private int refDivisor;
jpayne@68: 	private int queryDivisor;
jpayne@68: 	
jpayne@68: 	private int refSize;
jpayne@68: 	private int querySize;
jpayne@68: 
jpayne@68: 	private int hits1;
jpayne@68: 	private int qSeen1;
jpayne@68: 	private int rSeen1;
jpayne@68: 	
jpayne@68: 	private float ssuIdentity=-1;
jpayne@68: 	
jpayne@68: }