jpayne@68: package kmer;
jpayne@68: 
jpayne@68: import java.io.File;
jpayne@68: import java.util.ArrayList;
jpayne@68: import java.util.Arrays;
jpayne@68: import java.util.BitSet;
jpayne@68: import java.util.concurrent.atomic.AtomicLong;
jpayne@68: 
jpayne@68: import assemble.Contig;
jpayne@68: import bloom.KmerCountAbstract;
jpayne@68: import dna.AminoAcid;
jpayne@68: import fileIO.ByteStreamWriter;
jpayne@68: import fileIO.FileFormat;
jpayne@68: import fileIO.ReadWrite;
jpayne@68: import jgi.BBMerge;
jpayne@68: import shared.Parse;
jpayne@68: import shared.Parser;
jpayne@68: import shared.PreParser;
jpayne@68: import shared.Primes;
jpayne@68: import shared.ReadStats;
jpayne@68: import shared.Shared;
jpayne@68: import shared.Timer;
jpayne@68: import shared.Tools;
jpayne@68: import shared.TrimRead;
jpayne@68: import stream.ConcurrentReadInputStream;
jpayne@68: import stream.FastaReadInputStream;
jpayne@68: import stream.Read;
jpayne@68: import structures.ByteBuilder;
jpayne@68: import structures.IntList;
jpayne@68: import structures.ListNum;
jpayne@68: import structures.LongList;
jpayne@68: import ukmer.Kmer;
jpayne@68: 
jpayne@68: 
jpayne@68: /**
jpayne@68:  * Loads and holds kmers for Tadpole
jpayne@68:  * @author Brian Bushnell
jpayne@68:  * @date Jun 22, 2015
jpayne@68:  *
jpayne@68:  */
jpayne@68: public class KmerTableSet extends AbstractKmerTableSet {
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Code entrance from the command line.
jpayne@68: 	 * @param args Command line arguments
jpayne@68: 	 */
jpayne@68: 	public static void main(String[] args){
jpayne@68: 		Timer t=new Timer(), t2=new Timer();
jpayne@68: 		t.start();
jpayne@68: 		t2.start();
jpayne@68: 		KmerTableSet set=new KmerTableSet(args);
jpayne@68: 		t2.stop();
jpayne@68: 		outstream.println("Initialization Time:      \t"+t2);
jpayne@68: 		
jpayne@68: 		///And run it
jpayne@68: 		set.process(t);
jpayne@68: 		
jpayne@68: 		//Close the print stream if it was redirected
jpayne@68: 		Shared.closeStream(outstream);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Constructor.
jpayne@68: 	 * @param args Command line arguments
jpayne@68: 	 */
jpayne@68: 	private KmerTableSet(String[] args){
jpayne@68: 		this(args, 12);//+5 if using ownership and building contigs
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Constructor.
jpayne@68: 	 * @param args Command line arguments
jpayne@68: 	 */
jpayne@68: 	public KmerTableSet(String[] args, final int bytesPerKmer_){
jpayne@68: 		{//Preparse block for help, config files, and outstream
jpayne@68: 			PreParser pp=new PreParser(args, getClass(), false);
jpayne@68: 			args=pp.args;
jpayne@68: 			outstream=pp.outstream;
jpayne@68: 		}//TODO - no easy way to close outstream
jpayne@68: 		
jpayne@68: 		/* Initialize local variables with defaults */
jpayne@68: 		Parser parser=new Parser();
jpayne@68: 		boolean prealloc_=false;
jpayne@68: 		int k_=31;
jpayne@68: 		int ways_=-1;
jpayne@68: 		int filterMax_=2;
jpayne@68: 		boolean ecco_=false, merge_=false;
jpayne@68: 		boolean rcomp_=true;
jpayne@68: 		double minProb_=defaultMinprob;
jpayne@68: 		int tableType_=AbstractKmerTable.ARRAY1D;
jpayne@68: 		/* Parse arguments */
jpayne@68: 		for(int i=0; i<args.length; i++){
jpayne@68: 
jpayne@68: 			final String arg=args[i];
jpayne@68: 			String[] split=arg.split("=");
jpayne@68: 			String a=split[0].toLowerCase();
jpayne@68: 			String b=split.length>1 ? split[1] : null;
jpayne@68: 			
jpayne@68: 			if(Parser.parseCommonStatic(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(Parser.parseZip(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(Parser.parseQuality(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(Parser.parseFasta(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(parser.parseInterleaved(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(parser.parseTrim(arg, a, b)){
jpayne@68: 				//do nothing
jpayne@68: 			}else if(a.equals("in") || a.equals("in1")){
jpayne@68: 				in1.clear();
jpayne@68: 				if(b!=null){
jpayne@68: 					String[] s=b.split(",");
jpayne@68: 					for(String ss : s){
jpayne@68: 						in1.add(ss);
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("in2")){
jpayne@68: 				in2.clear();
jpayne@68: 				if(b!=null){
jpayne@68: 					String[] s=b.split(",");
jpayne@68: 					for(String ss : s){
jpayne@68: 						in2.add(ss);
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("extra")){
jpayne@68: 				extra.clear();
jpayne@68: 				if(b!=null){
jpayne@68: 					String[] s=b.split(",");
jpayne@68: 					for(String ss : s){
jpayne@68: 						extra.add(ss);
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("append") || a.equals("app")){
jpayne@68: 				append=ReadStats.append=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("overwrite") || a.equals("ow")){
jpayne@68: 				overwrite=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("initialsize")){
jpayne@68: 				initialSize=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("showstats") || a.equals("stats")){
jpayne@68: 				showStats=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("ways")){
jpayne@68: 				ways_=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("buflen") || a.equals("bufflen") || a.equals("bufferlength")){
jpayne@68: 				buflen=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("k")){
jpayne@68: 				assert(b!=null) : "\nk needs an integer value from 1 to 31, such as k=27.  Default is 31.\n";
jpayne@68: 				k_=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("threads") || a.equals("t")){
jpayne@68: 				THREADS=(b==null || b.equalsIgnoreCase("auto") ? Shared.threads() : Integer.parseInt(b));
jpayne@68: 			}else if(a.equals("showspeed") || a.equals("ss")){
jpayne@68: 				showSpeed=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("ecco")){
jpayne@68: 				ecco_=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("merge")){
jpayne@68: 				merge_=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("verbose")){
jpayne@68: //				assert(false) : "Verbose flag is currently static final; must be recompiled to change.";
jpayne@68: 				verbose=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("verbose2")){
jpayne@68: //				assert(false) : "Verbose flag is currently static final; must be recompiled to change.";
jpayne@68: 				verbose2=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("minprob")){
jpayne@68: 				minProb_=Double.parseDouble(b);
jpayne@68: 			}else if(a.equals("minprobprefilter") || a.equals("mpp")){
jpayne@68: 				minProbPrefilter=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("minprobmain") || a.equals("mpm")){
jpayne@68: 				minProbMain=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("reads") || a.startsWith("maxreads")){
jpayne@68: 				maxReads=Parse.parseKMG(b);
jpayne@68: 			}else if(a.equals("prealloc") || a.equals("preallocate")){
jpayne@68: 				if(b==null || b.length()<1 || Character.isLetter(b.charAt(0))){
jpayne@68: 					prealloc_=Parse.parseBoolean(b);
jpayne@68: 				}else{
jpayne@68: 					preallocFraction=Tools.max(0, Double.parseDouble(b));
jpayne@68: 					prealloc_=(preallocFraction>0);
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("prefilter")){
jpayne@68: 				if(b==null || b.length()<1 || !Tools.isDigit(b.charAt(0))){
jpayne@68: 					prefilter=Parse.parseBoolean(b);
jpayne@68: 				}else{
jpayne@68: 					filterMax_=Parse.parseIntKMG(b);
jpayne@68: 					prefilter=filterMax_>0;
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("prefiltersize") || a.equals("prefilterfraction") || a.equals("pff")){
jpayne@68: 				prefilterFraction=Tools.max(0, Double.parseDouble(b));
jpayne@68: 				assert(prefilterFraction<=1) : "prefiltersize must be 0-1, a fraction of total memory.";
jpayne@68: 				prefilter=prefilterFraction>0;
jpayne@68: 			}else if(a.equals("prehashes") || a.equals("hashes")){
jpayne@68: 				prehashes=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("prefilterpasses") || a.equals("prepasses")){
jpayne@68: 				assert(b!=null) : "Bad parameter: "+arg;
jpayne@68: 				if(b.equalsIgnoreCase("auto")){
jpayne@68: 					prepasses=-1;
jpayne@68: 				}else{
jpayne@68: 					prepasses=Parse.parseIntKMG(b);
jpayne@68: 				}
jpayne@68: 			}else if(a.equals("onepass")){
jpayne@68: 				onePass=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("passes")){
jpayne@68: 				int passes=Parse.parseIntKMG(b);
jpayne@68: 				onePass=(passes<2);
jpayne@68: 			}else if(a.equals("rcomp")){
jpayne@68: 				rcomp_=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("tabletype")){
jpayne@68: 				tableType_=Integer.parseInt(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equalsIgnoreCase("filterMemoryOverride") || a.equalsIgnoreCase("filterMemory") || 
jpayne@68: 					a.equalsIgnoreCase("prefilterMemory") || a.equalsIgnoreCase("filtermem")){
jpayne@68: 				filterMemoryOverride=Parse.parseKMG(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(IGNORE_UNKNOWN_ARGS){
jpayne@68: 				//Do nothing
jpayne@68: 			}else{
jpayne@68: 				throw new RuntimeException("Unknown parameter "+args[i]);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		{//Process parser fields
jpayne@68: 			Parser.processQuality();
jpayne@68: 			
jpayne@68: 			qtrimLeft=parser.qtrimLeft;
jpayne@68: 			qtrimRight=parser.qtrimRight;
jpayne@68: 			trimq=parser.trimq;
jpayne@68: 			trimE=parser.trimE();
jpayne@68: 			
jpayne@68: 			minAvgQuality=parser.minAvgQuality;
jpayne@68: 			minAvgQualityBases=parser.minAvgQualityBases;
jpayne@68: 			amino=Shared.AMINO_IN;
jpayne@68: 			if(amino){k_=Tools.min(k_, 12);}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(prepasses==0 || !prefilter){
jpayne@68: 			prepasses=0;
jpayne@68: 			prefilter=false;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 
jpayne@68: //		assert(false) : prepasses+", "+onePass+", "+prefilter;
jpayne@68: 		
jpayne@68: 		{
jpayne@68: 			long memory=Runtime.getRuntime().maxMemory();
jpayne@68: 			double xmsRatio=Shared.xmsRatio();
jpayne@68: //			long tmemory=Runtime.getRuntime().totalMemory();
jpayne@68: 			usableMemory=(long)Tools.max(((memory-96000000)*(xmsRatio>0.97 ? 0.82 : 0.72)), memory*0.45);
jpayne@68: 			if(prepasses==0 || !prefilter){
jpayne@68: 				filterMemory0=filterMemory1=0;
jpayne@68: 			}else if(filterMemoryOverride>0){
jpayne@68: 				filterMemory0=filterMemory1=filterMemoryOverride;
jpayne@68: 			}else{
jpayne@68: 				double low=Tools.min(prefilterFraction, 1-prefilterFraction);
jpayne@68: 				double high=1-low;
jpayne@68: 				if(prepasses<0 || (prepasses&1)==1){//odd passes
jpayne@68: 					filterMemory0=(long)(usableMemory*low);
jpayne@68: 					filterMemory1=(long)(usableMemory*high);
jpayne@68: 				}else{//even passes
jpayne@68: 					filterMemory0=(long)(usableMemory*high);
jpayne@68: 					filterMemory1=(long)(usableMemory*low);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			tableMemory=(long)(usableMemory*.95-filterMemory0);
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		tableType=tableType_;
jpayne@68: 		prealloc=prealloc_;
jpayne@68: 		bytesPerKmer=bytesPerKmer_;
jpayne@68: 		if(ways_<1){
jpayne@68: 			long maxKmers=(2*tableMemory)/bytesPerKmer;
jpayne@68: 			long minWays=Tools.min(10000, maxKmers/Integer.MAX_VALUE);
jpayne@68: 			ways_=(int)Tools.max(31, (int)(Tools.min(96, Shared.threads())*2.5), minWays);
jpayne@68: 			ways_=(int)Primes.primeAtLeast(ways_);
jpayne@68: 			assert(ways_>0);
jpayne@68: //			System.err.println("ways="+ways_);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* Set final variables; post-process and validate argument combinations */
jpayne@68: 		
jpayne@68: 		onePass=onePass&prefilter;
jpayne@68: 		ways=ways_;
jpayne@68: 		filterMax=Tools.min(filterMax_, 0x7FFFFFFF);
jpayne@68: 		ecco=ecco_;
jpayne@68: 		merge=merge_;
jpayne@68: 		minProb=(float)minProb_;
jpayne@68: 		rcomp=rcomp_;
jpayne@68: //		assert(false) : tableMemory+", "+bytesPerKmer+", "+prealloc+", "+preallocFraction;
jpayne@68: 		estimatedKmerCapacity=(long)((tableMemory*1.0/bytesPerKmer)*((prealloc ? preallocFraction : 0.81))*(HashArray.maxLoadFactorFinal*0.97));
jpayne@68: 		
jpayne@68: //		System.err.println("tableMemory="+tableMemory+", bytesPerKmer="+bytesPerKmer+", estimatedKmerCapacity="+estimatedKmerCapacity);
jpayne@68: 		KmerCountAbstract.minProb=(minProbPrefilter ? minProb : 0);
jpayne@68: 		k=k_;
jpayne@68: 		k2=k-1;
jpayne@68: 		int bitsPerBase=(amino ? 5 : 2);
jpayne@68: 		int mask=(amino ? 31 : 3);
jpayne@68: 		coreMask=(AbstractKmerTableSet.MASK_CORE ? ~(((-1L)<<(bitsPerBase*(k_-1)))|mask) : -1L);
jpayne@68: 		
jpayne@68: 		if(k<1 || k>31){throw new RuntimeException("\nk needs an integer value from 1 to 31, such as k=27.  Default is 31.\n");}
jpayne@68: 		
jpayne@68: 		if(initialSize<1){
jpayne@68: 			final long memOverWays=tableMemory/(bytesPerKmer*ways);
jpayne@68: 			final double mem2=(prealloc ? preallocFraction : 1)*tableMemory;
jpayne@68: 			initialSize=(prealloc || memOverWays<initialSizeDefault ? (int)Tools.min(2142000000, (long)(mem2/(bytesPerKmer*ways))) : initialSizeDefault);
jpayne@68: 			if(initialSize!=initialSizeDefault){
jpayne@68: 				System.err.println("Initial size set to "+initialSize);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* Adjust I/O settings and filenames */
jpayne@68: 		
jpayne@68: 		assert(FastaReadInputStream.settingsOK());
jpayne@68: 		 
jpayne@68: 		if(in1.isEmpty()){
jpayne@68: 			//throw new RuntimeException("Error - at least one input file is required.");
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		for(int i=0; i<in1.size(); i++){
jpayne@68: 			String s=in1.get(i);
jpayne@68: 			if(s!=null && s.contains("#") && !new File(s).exists()){
jpayne@68: 				int pound=s.lastIndexOf('#');
jpayne@68: 				String a=s.substring(0, pound);
jpayne@68: 				String b=s.substring(pound+1);
jpayne@68: 				in1.set(i, a+1+b);
jpayne@68: 				in2.add(a+2+b);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(!extra.isEmpty()){
jpayne@68: 			ArrayList<String> temp=(ArrayList<String>) extra.clone();
jpayne@68: 			extra.clear();
jpayne@68: 			for(String s : temp){
jpayne@68: 				if(s!=null && s.contains("#") && !new File(s).exists()){
jpayne@68: 					int pound=s.lastIndexOf('#');
jpayne@68: 					String a=s.substring(0, pound);
jpayne@68: 					String b=s.substring(pound+1);
jpayne@68: 					extra.add(a);
jpayne@68: 					extra.add(b);
jpayne@68: 				}else{
jpayne@68: 					extra.add(s);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		{
jpayne@68: 			boolean allowDuplicates=true;
jpayne@68: 			if(!Tools.testInputFiles(allowDuplicates, true, in1, in2, extra)){
jpayne@68: 				throw new RuntimeException("\nCan't read some input files.\n");  
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		assert(THREADS>0);
jpayne@68: 		
jpayne@68: 		if(DISPLAY_PROGRESS){
jpayne@68: //			assert(false);
jpayne@68: 			outstream.println("Initial:");
jpayne@68: 			outstream.println("Ways="+ways+", initialSize="+initialSize+", prefilter="+(prefilter ? "t" : "f")+", prealloc="+(prealloc ? (""+preallocFraction) : "f"));
jpayne@68: 			Shared.printMemory();
jpayne@68: 			outstream.println();
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Outer Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public void clear(){
jpayne@68: 		tables=null;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public void allocateTables(){
jpayne@68: 		assert(tables==null);
jpayne@68: 		tables=null;
jpayne@68: 		final long coreMask=(MASK_CORE ? ~(((-1L)<<(2*(k-1)))|3) : -1L);
jpayne@68: 		
jpayne@68: 		ScheduleMaker scheduleMaker=new ScheduleMaker(ways, bytesPerKmer, prealloc, 
jpayne@68: 				(prealloc ? preallocFraction : 1.0), -1, (prefilter ? prepasses : 0), prefilterFraction, filterMemoryOverride);
jpayne@68: 		int[] schedule=scheduleMaker.makeSchedule();
jpayne@68: 		//assert(false) : preallocFraction+", "+ways+", "+Arrays.toString(schedule);
jpayne@68: //		System.err.println("DEBUG "+preallocFraction+", "+ways+", "+Arrays.toString(schedule));
jpayne@68: 		tables=AbstractKmerTable.preallocate(ways, tableType, schedule, coreMask);
jpayne@68: 		
jpayne@68: //		tables=AbstractKmerTable.preallocate(ways, tableType, initialSize, coreMask, (!prealloc || preallocFraction<1));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Load reads into tables, using multiple LoadThread.
jpayne@68: 	 */
jpayne@68: 	@Override
jpayne@68: 	public long loadKmers(String fname1, String fname2){
jpayne@68: 		
jpayne@68: 		/* Create read input stream */
jpayne@68: 		final ConcurrentReadInputStream cris;
jpayne@68: 		{
jpayne@68: 			FileFormat ff1=FileFormat.testInput(fname1, FileFormat.FASTQ, null, true, true);
jpayne@68: 			FileFormat ff2=FileFormat.testInput(fname2, FileFormat.FASTQ, null, true, true);
jpayne@68: 			cris=ConcurrentReadInputStream.getReadInputStream(maxReads, false, ff1, ff2);
jpayne@68: 			cris.start(); //4567
jpayne@68: 		}
jpayne@68: 		
jpayne@68: //		/* Optionally skip the first reads, since initial reads may have lower quality */
jpayne@68: //		if(skipreads>0){
jpayne@68: //			long skipped=0;
jpayne@68: //
jpayne@68: //			ListNum<Read> ln=cris.nextList();
jpayne@68: //			ArrayList<Read> reads=(ln!=null ? ln.list : null);
jpayne@68: //
jpayne@68: //			while(skipped<skipreads && reads!=null && reads.size()>0){
jpayne@68: //				skipped+=reads.size();
jpayne@68: //
jpayne@68: //				cris.returnList(ln);
jpayne@68: //				ln=cris.nextList();
jpayne@68: //				reads=(ln!=null ? ln.list : null);
jpayne@68: //			}
jpayne@68: //			cris.returnList(ln);
jpayne@68: //			if(reads==null || reads.isEmpty()){
jpayne@68: //				ReadWrite.closeStreams(cris);
jpayne@68: //				System.err.println("Skipped all of the reads.");
jpayne@68: //				System.exit(0);
jpayne@68: //			}
jpayne@68: //		}
jpayne@68: 		
jpayne@68: 		/* Create ProcessThreads */
jpayne@68: 		ArrayList<LoadThread> alpt=new ArrayList<LoadThread>(THREADS);
jpayne@68: 		for(int i=0; i<THREADS; i++){alpt.add(new LoadThread(cris));}
jpayne@68: 		for(LoadThread pt : alpt){pt.start();}
jpayne@68: 		
jpayne@68: 		long added=0;
jpayne@68: 		
jpayne@68: 		/* Wait for threads to die, and gather statistics */
jpayne@68: 		for(LoadThread pt : alpt){
jpayne@68: 			while(pt.getState()!=Thread.State.TERMINATED){
jpayne@68: 				try {
jpayne@68: 					pt.join();
jpayne@68: 				} catch (InterruptedException e) {
jpayne@68: 					// TODO Auto-generated catch block
jpayne@68: 					e.printStackTrace();
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			added+=pt.added;
jpayne@68: 			
jpayne@68: 			readsIn+=pt.readsInT;
jpayne@68: 			basesIn+=pt.basesInT;
jpayne@68: 			lowqReads+=pt.lowqReadsT;
jpayne@68: 			lowqBases+=pt.lowqBasesT;
jpayne@68: 			readsTrimmed+=pt.readsTrimmedT;
jpayne@68: 			basesTrimmed+=pt.basesTrimmedT;
jpayne@68: 			kmersIn+=pt.kmersInT;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* Shut down I/O streams; capture error status */
jpayne@68: 		errorState|=ReadWrite.closeStreams(cris);
jpayne@68: 		return added;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Classes        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/**
jpayne@68: 	 * Loads kmers.
jpayne@68: 	 */
jpayne@68: 	private class LoadThread extends Thread{
jpayne@68: 		
jpayne@68: 		/**
jpayne@68: 		 * Constructor
jpayne@68: 		 * @param cris_ Read input stream
jpayne@68: 		 */
jpayne@68: 		public LoadThread(ConcurrentReadInputStream cris_){
jpayne@68: 			cris=cris_;
jpayne@68: 			table=new HashBuffer(tables, buflen, k, false, true);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public void run(){
jpayne@68: 			ListNum<Read> ln=cris.nextList();
jpayne@68: 			ArrayList<Read> reads=(ln!=null ? ln.list : null);
jpayne@68: 			
jpayne@68: 			//While there are more reads lists...
jpayne@68: 			while(ln!=null && reads!=null && reads.size()>0){//ln!=null prevents a compiler potential null access warning
jpayne@68: 				
jpayne@68: 				//For each read (or pair) in the list...
jpayne@68: 				for(int i=0; i<reads.size(); i++){
jpayne@68: 					Read r1=reads.get(i);
jpayne@68: 					Read r2=r1.mate;
jpayne@68: 					
jpayne@68: 					if(!r1.validated()){r1.validate(true);}
jpayne@68: 					if(r2!=null && !r2.validated()){r2.validate(true);}
jpayne@68: 					
jpayne@68: 					if(verbose){System.err.println("Considering read "+r1.id+" "+new String(r1.bases));}
jpayne@68: 					
jpayne@68: 					readsInT++;
jpayne@68: 					basesInT+=r1.length();
jpayne@68: 					if(r2!=null){
jpayne@68: 						readsInT++;
jpayne@68: 						basesInT+=r2.length();
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					if(maxNs<Integer.MAX_VALUE){
jpayne@68: 						if(r1!=null && r1.countUndefined()>maxNs){r1.setDiscarded(true);}
jpayne@68: 						if(r2!=null && r2.countUndefined()>maxNs){r2.setDiscarded(true);}
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					//Determine whether to discard the reads based on average quality
jpayne@68: 					if(minAvgQuality>0){
jpayne@68: 						if(r1!=null && r1.quality!=null && r1.avgQuality(false, minAvgQualityBases)<minAvgQuality){r1.setDiscarded(true);}
jpayne@68: 						if(r2!=null && r2.quality!=null && r2.avgQuality(false, minAvgQualityBases)<minAvgQuality){r2.setDiscarded(true);}
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					if(r1!=null){
jpayne@68: 						if(qtrimLeft || qtrimRight){
jpayne@68: 							int x=TrimRead.trimFast(r1, qtrimLeft, qtrimRight, trimq, trimE, 1, false);
jpayne@68: 							basesTrimmedT+=x;
jpayne@68: 							readsTrimmedT+=(x>0 ? 1 : 0);
jpayne@68: 						}
jpayne@68: 						if(r1.length()<k){r1.setDiscarded(true);}
jpayne@68: 					}
jpayne@68: 					if(r2!=null){
jpayne@68: 						if(qtrimLeft || qtrimRight){
jpayne@68: 							int x=TrimRead.trimFast(r2, qtrimLeft, qtrimRight, trimq, trimE, 1, false);
jpayne@68: 							basesTrimmedT+=x;
jpayne@68: 							readsTrimmedT+=(x>0 ? 1 : 0);
jpayne@68: 						}
jpayne@68: 						if(r2.length()<k){r2.setDiscarded(true);}
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					if((ecco || merge) && r1!=null && r2!=null && !r1.discarded() && !r2.discarded()){
jpayne@68: 						if(merge){
jpayne@68: 							final int insert=BBMerge.findOverlapStrict(r1, r2, false);
jpayne@68: 							if(insert>0){
jpayne@68: 								r2.reverseComplement();
jpayne@68: 								r1=r1.joinRead(insert);
jpayne@68: 								r2=null;
jpayne@68: 							}
jpayne@68: 						}else if(ecco){
jpayne@68: 							BBMerge.findOverlapStrict(r1, r2, true);
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					if(minLen>0){
jpayne@68: 						if(r1!=null && r1.length()<minLen){r1.setDiscarded(true);}
jpayne@68: 						if(r2!=null && r2.length()<minLen){r2.setDiscarded(true);}
jpayne@68: 					}
jpayne@68: 
jpayne@68: 					if(r1!=null){
jpayne@68: 						if(r1.discarded()){
jpayne@68: 							lowqBasesT+=r1.length();
jpayne@68: 							lowqReadsT++;
jpayne@68: 						}else{
jpayne@68: 							long temp=addKmersToTable(r1);
jpayne@68: 							added+=temp;
jpayne@68: 							if(verbose){System.err.println("A: Added "+temp);}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 					if(r2!=null){
jpayne@68: 						if(r2.discarded()){
jpayne@68: 							lowqBasesT+=r2.length();
jpayne@68: 							lowqReadsT++;
jpayne@68: 						}else{
jpayne@68: 							long temp=addKmersToTable(r2);
jpayne@68: 							added+=temp;
jpayne@68: 							if(verbose){System.err.println("B: Added "+temp);}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 				
jpayne@68: 				//Fetch a new read list
jpayne@68: 				cris.returnList(ln);
jpayne@68: 				ln=cris.nextList();
jpayne@68: 				reads=(ln!=null ? ln.list : null);
jpayne@68: 			}
jpayne@68: 			cris.returnList(ln);
jpayne@68: 			long temp=table.flush();
jpayne@68: 			if(verbose){System.err.println("Flush: Added "+temp);}
jpayne@68: 			added+=temp;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		private final int addKmersToTableAA(final Read r){
jpayne@68: 			if(r==null || r.bases==null){return 0;}
jpayne@68: 			final float minProb2=(minProbMain ? minProb : 0);
jpayne@68: 			final byte[] bases=r.bases;
jpayne@68: 			final byte[] quals=r.quality;
jpayne@68: 			final int shift=5*k;
jpayne@68: 			final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 			long kmer=0;
jpayne@68: 			int created=0;
jpayne@68: 			int len=0;
jpayne@68: 			
jpayne@68: 			if(bases==null || bases.length<k){return -1;}
jpayne@68: 			
jpayne@68: 			/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 			float prob=1;
jpayne@68: 			for(int i=0; i<bases.length; i++){
jpayne@68: 				final byte b=bases[i];
jpayne@68: 				final long x=AminoAcid.acidToNumber[b];
jpayne@68: 
jpayne@68: 				//Update kmers
jpayne@68: 				kmer=((kmer<<5)|x)&mask;
jpayne@68: 
jpayne@68: 				if(minProb2>0 && quals!=null){//Update probability
jpayne@68: 					prob=prob*PROB_CORRECT[quals[i]];
jpayne@68: 					if(len>k){
jpayne@68: 						byte oldq=quals[i-k];
jpayne@68: 						prob=prob*PROB_CORRECT_INVERSE[oldq];
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 
jpayne@68: 				//Handle Ns
jpayne@68: 				if(x<0){
jpayne@68: 					len=0;
jpayne@68: 					kmer=0;
jpayne@68: 					prob=1;
jpayne@68: 				}else{len++;}
jpayne@68: 
jpayne@68: 				if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 				if(len>=k && prob>=minProb2){
jpayne@68: 					kmersInT++;
jpayne@68: 					final long key=kmer;
jpayne@68: 					if(!prefilter || prefilterArray.read(key)>filterMax2){
jpayne@68: 						int temp=table.incrementAndReturnNumCreated(key, 1);
jpayne@68: 						created+=temp;
jpayne@68: 						if(verbose){System.err.println("C: Added "+temp);}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			return created;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		private final int addKmersToTable(final Read r){
jpayne@68: 			if(amino){return addKmersToTableAA(r);}
jpayne@68: 			if(onePass){return addKmersToTable_onePass(r);}
jpayne@68: 			if(r==null || r.bases==null){return 0;}
jpayne@68: 			final float minProb2=(minProbMain ? minProb : 0);
jpayne@68: 			final byte[] bases=r.bases;
jpayne@68: 			final byte[] quals=r.quality;
jpayne@68: 			final int shift=2*k;
jpayne@68: 			final int shift2=shift-2;
jpayne@68: 			final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 			long kmer=0;
jpayne@68: 			long rkmer=0;
jpayne@68: 			int created=0;
jpayne@68: 			int len=0;
jpayne@68: 			
jpayne@68: 			if(bases==null || bases.length<k){return -1;}
jpayne@68: 			
jpayne@68: 			/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 			float prob=1;
jpayne@68: 			for(int i=0; i<bases.length; i++){
jpayne@68: 				final byte b=bases[i];
jpayne@68: 				final long x=AminoAcid.baseToNumber[b];
jpayne@68: 				final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 
jpayne@68: 				//Update kmers
jpayne@68: 				kmer=((kmer<<2)|x)&mask;
jpayne@68: 				rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 
jpayne@68: 				if(minProb2>0 && quals!=null){//Update probability
jpayne@68: 					prob=prob*PROB_CORRECT[quals[i]];
jpayne@68: 					if(len>k){
jpayne@68: 						byte oldq=quals[i-k];
jpayne@68: 						prob=prob*PROB_CORRECT_INVERSE[oldq];
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 
jpayne@68: 				//Handle Ns
jpayne@68: 				if(x<0){
jpayne@68: 					len=0;
jpayne@68: 					kmer=rkmer=0;
jpayne@68: 					prob=1;
jpayne@68: 				}else{len++;}
jpayne@68: 
jpayne@68: 				if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 				if(len>=k && prob>=minProb2){
jpayne@68: 					kmersInT++;
jpayne@68: 					final long key=toValue(kmer, rkmer);
jpayne@68: 					if(!prefilter || prefilterArray.read(key)>filterMax2){
jpayne@68: 						int temp=table.incrementAndReturnNumCreated(key, 1);
jpayne@68: 						created+=temp;
jpayne@68: 						if(verbose){System.err.println("C: Added "+temp);}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			return created;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		
jpayne@68: 		private final int addKmersToTable_onePass(final Read r){
jpayne@68: 			assert(prefilter);
jpayne@68: 			if(r==null || r.bases==null){return 0;}
jpayne@68: 			final byte[] bases=r.bases;
jpayne@68: 			final byte[] quals=r.quality;
jpayne@68: 			final int shift=2*k;
jpayne@68: 			final int shift2=shift-2;
jpayne@68: 			final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 			long kmer=0;
jpayne@68: 			long rkmer=0;
jpayne@68: 			int created=0;
jpayne@68: 			int len=0;
jpayne@68: 
jpayne@68: 			if(bases==null || bases.length<k){return -1;}
jpayne@68: 			
jpayne@68: 			/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 			float prob=1;
jpayne@68: 			for(int i=0; i<bases.length; i++){
jpayne@68: 				final byte b=bases[i];
jpayne@68: 				final long x=AminoAcid.baseToNumber[b];
jpayne@68: 				final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 
jpayne@68: 				//Update kmers
jpayne@68: 				kmer=((kmer<<2)|x)&mask;
jpayne@68: 				rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 
jpayne@68: 				if(minProb>0 && quals!=null){//Update probability
jpayne@68: 					prob=prob*PROB_CORRECT[quals[i]];
jpayne@68: 					if(len>k){
jpayne@68: 						byte oldq=quals[i-k];
jpayne@68: 						prob=prob*PROB_CORRECT_INVERSE[oldq];
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 
jpayne@68: 				//Handle Ns
jpayne@68: 				if(x<0){
jpayne@68: 					len=0;
jpayne@68: 					kmer=rkmer=0;
jpayne@68: 					prob=1;
jpayne@68: 				}else{len++;}
jpayne@68: 
jpayne@68: 				if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 				if(len>=k){
jpayne@68: 					final long key=toValue(kmer, rkmer);
jpayne@68: 					int count=prefilterArray.incrementAndReturnUnincremented(key, 1);
jpayne@68: 					if(count>=filterMax2){
jpayne@68: 						int temp=table.incrementAndReturnNumCreated(key, 1);
jpayne@68: 						created+=temp;
jpayne@68: 						if(verbose){System.err.println("D: Added "+temp);}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return created;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/*--------------------------------------------------------------*/
jpayne@68: 		
jpayne@68: 		/** Input read stream */
jpayne@68: 		private final ConcurrentReadInputStream cris;
jpayne@68: 		
jpayne@68: 		private final HashBuffer table;
jpayne@68: 		
jpayne@68: 		public long added=0;
jpayne@68: 		
jpayne@68: 		private long readsInT=0;
jpayne@68: 		private long basesInT=0;
jpayne@68: 		private long lowqReadsT=0;
jpayne@68: 		private long lowqBasesT=0;
jpayne@68: 		private long readsTrimmedT=0;
jpayne@68: 		private long basesTrimmedT=0;
jpayne@68: 		private long kmersInT=0;
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------          Convenience         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public void regenerateKmers(byte[] bases, LongList kmers, IntList counts, final int a){
jpayne@68: 		final int loc=a+k;
jpayne@68: 		final int lim=Tools.min(counts.size, a+k+1);
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=kmers.get(a);
jpayne@68: 		long rkmer=rcomp(kmer);
jpayne@68: 		int len=k;
jpayne@68: 		
jpayne@68: //		assert(false) : a+", "+loc+", "+lim;
jpayne@68: 		
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=loc, j=a+1; j<lim; i++, j++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			
jpayne@68: 			if(len>=k){
jpayne@68: 				assert(kmers.get(j)!=kmer);
jpayne@68: 				kmers.set(j, kmer);
jpayne@68: 				int count=getCount(kmer, rkmer);
jpayne@68: 				counts.set(j, count);
jpayne@68: 			}else{
jpayne@68: 				kmers.set(j, -1);
jpayne@68: 				counts.set(j, 0);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public int regenerateCounts(byte[] bases, IntList counts, final Kmer dummy, BitSet changed){
jpayne@68: 		assert(!changed.isEmpty());
jpayne@68: 		final int firstBase=changed.nextSetBit(0), lastBase=changed.length()-1;
jpayne@68: 		final int ca=firstBase-k;
jpayne@68: //		final int b=changed.nextSetBit(0);ca+kbig; //first base changed
jpayne@68: 		final int firstCount=Tools.max(firstBase-k+1, 0), lastCount=Tools.min(counts.size-1, lastBase); //count limit
jpayne@68: //		System.err.println("ca="+ca+", b="+b+", lim="+lim);
jpayne@68: //		System.err.println("Regen from count "+(ca+1)+"-"+lim);
jpayne@68: 		
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		int valid=0;
jpayne@68: 		
jpayne@68: //		System.err.println("ca="+ca+", b="+b+", lim="+lim+", "+counts);
jpayne@68: 		
jpayne@68: 		for(int i=Tools.max(0, firstBase-k+1), lim=Tools.min(lastBase+k-1, bases.length-1); i<=lim; i++){
jpayne@68: 			final byte base=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[base];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[base];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{
jpayne@68: 				len++;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			final int c=i-k+1;
jpayne@68: 			if(i>=firstBase){
jpayne@68: 				if(len>=k){
jpayne@68: 					valid++;
jpayne@68: 					int count=getCount(kmer, rkmer);
jpayne@68: 					counts.set(c, count);
jpayne@68: 				}else if(c>=0){
jpayne@68: 					counts.set(c, 0);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		return valid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public int fillSpecificCounts(byte[] bases, IntList counts, BitSet positions, final Kmer dummy){
jpayne@68: 		final int b=k-1;
jpayne@68: 		final int lim=(positions==null ? bases.length : Tools.min(bases.length, positions.length()+k-1));
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		int valid=0;
jpayne@68: 		
jpayne@68: 		counts.clear();
jpayne@68: 		
jpayne@68: 		for(int i=0, j=-b; i<lim; i++, j++){
jpayne@68: 			final byte base=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[base];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[base];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{
jpayne@68: 				len++;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(j>=0){
jpayne@68: 				if(len>=k && (positions==null || positions.get(j))){
jpayne@68: 					valid++;
jpayne@68: 					int count=getCount(kmer, rkmer);
jpayne@68: 					assert(i-k+1==counts.size()) : "j="+j+", counts="+counts.size()+", b="+(b)+", (i-k+1)="+(i-k+1);
jpayne@68: 					assert(j==counts.size());
jpayne@68: 					counts.add(Tools.max(count, 0));
jpayne@68: //					counts.set(i-k+1, count);
jpayne@68: 				}else{
jpayne@68: 					counts.add(0);
jpayne@68: //					counts.set(i-k+1, 0);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 
jpayne@68: //		assert((counts.size==0 && bases.length<k) || counts.size==bases.length-k+1) : bases.length+", "+k+", "+counts.size;
jpayne@68: 		assert((counts.size==0 && bases.length<k) || counts.size==lim-k+1) : bases.length+", "+k+", "+counts.size;
jpayne@68: 		
jpayne@68: 		return valid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int regenerateCounts(byte[] bases, IntList counts, final int ca){
jpayne@68: 		final int b=ca+k-1;
jpayne@68: 		final int lim=Tools.min(bases.length, b+k+1);
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		int valid=0;
jpayne@68: 		
jpayne@68: 		final int clen=counts.size;
jpayne@68: 		
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts.
jpayne@68: 		 * i is an index in the base array, j is an index in the count array. */
jpayne@68: 		for(int i=b-k+1; i<lim; i++){
jpayne@68: 			final byte base=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[base];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[base];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{
jpayne@68: 				len++;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(i>=b){
jpayne@68: 				if(len>=k){
jpayne@68: 					valid++;
jpayne@68: 					int count=getCount(kmer, rkmer);
jpayne@68: 					counts.set(i-k+1, count);
jpayne@68: 				}else{
jpayne@68: 					counts.set(i-k+1, 0);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		assert((counts.size==0 && bases.length<k) || counts.size==bases.length-k+1) : bases.length+", "+k+", "+counts.size;
jpayne@68: 		assert(clen==counts.size);
jpayne@68: 		
jpayne@68: 		return valid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Returns number of valid kmers */
jpayne@68: 	public int fillKmers(byte[] bases, LongList kmers){
jpayne@68: 		final int blen=bases.length;
jpayne@68: 		if(blen<k){return 0;}
jpayne@68: 		final int min=k-1;
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		int valid=0;
jpayne@68: 
jpayne@68: 		kmers.clear();
jpayne@68: 
jpayne@68: 		/* Loop through the bases, maintaining a forward kmer via bitshifts */
jpayne@68: 		for(int i=0; i<blen; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			assert(b>=0) : Arrays.toString(bases);
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(i>=min){
jpayne@68: 				if(len>=k){
jpayne@68: 					kmers.add(kmer);
jpayne@68: 					valid++;
jpayne@68: 				}else{
jpayne@68: 					kmers.add(-1);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return valid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void fillCounts(LongList kmers, IntList counts){
jpayne@68: 		counts.clear();
jpayne@68: 		for(int i=0; i<kmers.size; i++){
jpayne@68: 			long kmer=kmers.get(i);
jpayne@68: 			if(kmer>=0){
jpayne@68: 				long rkmer=rcomp(kmer);
jpayne@68: 				int count=getCount(kmer, rkmer);
jpayne@68: 				counts.add(count);
jpayne@68: 			}else{
jpayne@68: 				counts.add(0);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Helper Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public long regenerate(final int limit){
jpayne@68: 		long sum=0;
jpayne@68: 		for(AbstractKmerTable akt : tables){
jpayne@68: 			sum+=akt.regenerate(limit);
jpayne@68: 		}
jpayne@68: 		return sum;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public HashArray1D getTableForKey(long key){
jpayne@68: 		return (HashArray1D) tables[kmerToWay(key)];
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public HashArray1D getTable(int tnum){
jpayne@68: 		return (HashArray1D) tables[tnum];
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public long[] fillHistogram(int histMax) {
jpayne@68: 		return HistogramMaker.fillHistogram(tables, histMax);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public void countGC(long[] gcCounts, int max) {
jpayne@68: 		for(AbstractKmerTable set : tables){
jpayne@68: 			set.countGC(gcCounts, max);
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public void initializeOwnership(){
jpayne@68: 		OwnershipThread.initialize(tables);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public void clearOwnership(){
jpayne@68: 		OwnershipThread.clear(tables);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public long rightmostKmer(final ByteBuilder bb){
jpayne@68: 		return rightmostKmer(bb.array, bb.length());
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public long rightmostKmer(final byte[] bases, final int blen){
jpayne@68: 		if(blen<k){return -1;}
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0;
jpayne@68: 		long rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts, to get the rightmost kmer */
jpayne@68: 		{
jpayne@68: 			for(int i=blen-k; i<blen; i++){
jpayne@68: 				final byte b=bases[i];
jpayne@68: 				final long x=AminoAcid.baseToNumber[b];
jpayne@68: 				final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 				kmer=((kmer<<2)|x)&mask;
jpayne@68: 				rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 				if(x<0){
jpayne@68: 					len=0;
jpayne@68: 					kmer=rkmer=0;
jpayne@68: 				}else{len++;}
jpayne@68: 				if(verbose){outstream.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(len<k){return -1;}
jpayne@68: 		else{assert(len==k);}
jpayne@68: 		return kmer;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public long leftmostKmer(final ByteBuilder bb){
jpayne@68: 		return leftmostKmer(bb.array, bb.length());
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public long leftmostKmer(final byte[] bases, final int blen){
jpayne@68: 		if(blen<k){return -1;}
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0;
jpayne@68: 		long rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts, to get the leftmost kmer */
jpayne@68: 		{
jpayne@68: 			for(int i=0; i<k; i++){
jpayne@68: 				final byte b=bases[i];
jpayne@68: 				final long x=AminoAcid.baseToNumber[b];
jpayne@68: 				final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 				kmer=((kmer<<2)|x)&mask;
jpayne@68: 				rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 				if(x<0){
jpayne@68: 					len=0;
jpayne@68: 					kmer=rkmer=0;
jpayne@68: 				}else{len++;}
jpayne@68: 				if(verbose){outstream.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(len<k){return -1;}
jpayne@68: 		else{assert(len==k);}
jpayne@68: 		return kmer;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean doubleClaim(final ByteBuilder bb, final int id){
jpayne@68: 		return doubleClaim(bb.array, bb.length(), id);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Ensures there can be only one owner. */
jpayne@68: 	public boolean doubleClaim(final byte[] bases, final int blength, final int id){
jpayne@68: 		boolean success=claim(bases, blength, id, true);
jpayne@68: 		if(verbose){outstream.println("success1="+success+", id="+id+", blength="+blength);}
jpayne@68: 		if(!success){return false;}
jpayne@68: 		success=claim(bases, blength, id+CLAIM_OFFSET, true);
jpayne@68: 		if(verbose){outstream.println("success2="+success+", id="+id+", blength="+blength);}
jpayne@68: 		return success;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean claim(final ByteBuilder bb, final int id, final boolean exitEarly){
jpayne@68: 		return claim(bb.array, bb.length(), id, exitEarly);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float calcCoverage(final byte[] bases, final int blength){
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		long sum=0, max=0;
jpayne@68: 		int kmers=0;
jpayne@68: 		
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=0; i<blength; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(len>=k){
jpayne@68: 				int count=getCount(kmer, rkmer);
jpayne@68: 				sum+=count;
jpayne@68: 				max=Tools.max(count, max);
jpayne@68: 				kmers++;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return sum==0 ? 0 : sum/(float)kmers;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public float calcCoverage(final Contig contig){
jpayne@68: 		final byte[] bases=contig.bases;
jpayne@68: 		if(bases.length<k){return 0;}
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		long sum=0;
jpayne@68: 		int max=0, min=Integer.MAX_VALUE;
jpayne@68: 		int kmers=0;
jpayne@68: 		
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=0; i<bases.length; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(len>=k){
jpayne@68: 				int count=getCount(kmer, rkmer);
jpayne@68: 				sum+=count;
jpayne@68: 				max=Tools.max(count, max);
jpayne@68: 				min=Tools.min(count, min);
jpayne@68: 				kmers++;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		contig.coverage=sum==0 ? 0 : sum/(float)kmers;
jpayne@68: 		contig.maxCov=max;
jpayne@68: 		contig.minCov=sum==0 ? 0 : min;
jpayne@68: 		return contig.coverage;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean claim(final byte[] bases, final int blength, final int id, boolean exitEarly){
jpayne@68: 		if(verbose){outstream.println("Thread "+id+" claim start.");}
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		boolean success=true;
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=0; i<blength && success; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 			if(len>=k){
jpayne@68: 				success=claim(kmer, rkmer, id/*, rid, i*/);
jpayne@68: 				success=(success || !exitEarly);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return success;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean claim(final long kmer, final long rkmer, final int id/*, final long rid, final int pos*/){
jpayne@68: 		//TODO: rid and pos are just for debugging.
jpayne@68: 		final long key=toValue(kmer, rkmer);
jpayne@68: 		final int way=kmerToWay(key);
jpayne@68: 		final AbstractKmerTable table=tables[way];
jpayne@68: 		final int count=table.getValue(key);
jpayne@68: 		assert(count==-1 || count>0) : count;
jpayne@68: //		if(verbose  /*|| true*/){outstream.println("Count="+count+".");}
jpayne@68: 		if(count<0){return true;}
jpayne@68: 		assert(count>0) : count;
jpayne@68: 		final int owner=table.setOwner(key, id);
jpayne@68: 		if(verbose){outstream.println("owner="+owner+".");}
jpayne@68: //		assert(owner==id) : id+", "+owner+", "+rid+", "+pos;
jpayne@68: 		return owner==id;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void release(ByteBuilder bb, final int id){
jpayne@68: 		release(bb.array, bb.length(), id);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void release(final byte[] bases, final int blength, final int id){
jpayne@68: 		if(verbose  /*|| true*/){outstream.println("*Thread "+id+" release start.");}
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=0; i<blength; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 			if(len>=k){
jpayne@68: 				release(kmer, rkmer, id);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean release(final long kmer, final long rkmer, final int id){
jpayne@68: 		return release(toValue(kmer, rkmer), id);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public boolean release(final long key, final int id){
jpayne@68: 		final int way=kmerToWay(key);
jpayne@68: 		final AbstractKmerTable table=tables[way];
jpayne@68: 		final int count=table.getValue(key);
jpayne@68: //		if(verbose  /*|| true*/){outstream.println("Count="+count+".");}
jpayne@68: 		if(count<1){return true;}
jpayne@68: 		return table.clearOwner(key, id);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int findOwner(ByteBuilder bb, final int id){
jpayne@68: 		return findOwner(bb.array, bb.length(), id);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int findOwner(final byte[] bases, final int blength, final int id){
jpayne@68: 		final int shift=2*k;
jpayne@68: 		final int shift2=shift-2;
jpayne@68: 		final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
jpayne@68: 		long kmer=0, rkmer=0;
jpayne@68: 		int len=0;
jpayne@68: 		int maxOwner=-1;
jpayne@68: 		/* Loop through the bases, maintaining a forward and reverse kmer via bitshifts */
jpayne@68: 		for(int i=0; i<blength; i++){
jpayne@68: 			final byte b=bases[i];
jpayne@68: 			final long x=AminoAcid.baseToNumber[b];
jpayne@68: 			final long x2=AminoAcid.baseToComplementNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			rkmer=((rkmer>>>2)|(x2<<shift2))&mask;
jpayne@68: 			if(x<0){
jpayne@68: 				len=0;
jpayne@68: 				kmer=rkmer=0;
jpayne@68: 			}else{len++;}
jpayne@68: 			if(verbose){System.err.println("Scanning i="+i+", len="+len+", kmer="+kmer+", rkmer="+rkmer+"\t"+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 			if(len>=k){
jpayne@68: 				int owner=findOwner(kmer, rkmer);
jpayne@68: 				maxOwner=Tools.max(owner, maxOwner);
jpayne@68: 				if(maxOwner>id){break;}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return maxOwner;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int findOwner(final long kmer){
jpayne@68: 		return findOwner(kmer, rcomp(kmer));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int findOwner(final long kmer, final long rkmer){
jpayne@68: 		final long key=toValue(kmer, rkmer);
jpayne@68: 		final int way=kmerToWay(key);
jpayne@68: 		final AbstractKmerTable table=tables[way];
jpayne@68: 		final int count=table.getValue(key);
jpayne@68: 		if(count<0){return -1;}
jpayne@68: 		final int owner=table.getOwner(key);
jpayne@68: 		return owner;
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	public int getCount(long kmer, long rkmer){
jpayne@68: 		long key=toValue(kmer, rkmer);
jpayne@68: 		int way=kmerToWay(key);
jpayne@68: 		return tables[way].getValue(key);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int getCount(long key){
jpayne@68: 		int way=kmerToWay(key);
jpayne@68: 		return tables[way].getValue(key);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------          Fill Counts         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public int fillRightCounts(long kmer, long rkmer, int[] counts, long mask, int shift2){
jpayne@68: 		if(FAST_FILL && MASK_CORE && k>2/*((k&1)==1)*/){
jpayne@68: 			return fillRightCounts_fast(kmer, rkmer, counts, mask, shift2);
jpayne@68: 		}else{
jpayne@68: 			return fillRightCounts_safe(kmer, rkmer, counts, mask, shift2);
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int fillRightCounts_fast(final long kmer0, final long rkmer0, int[] counts,
jpayne@68: 			long mask, int shift2){
jpayne@68: //		assert((k&1)==1) : k;
jpayne@68: 		assert(MASK_CORE);
jpayne@68: 		final long kmer=(kmer0<<2)&mask;
jpayne@68: 		final long rkmer=(rkmer0>>>2);
jpayne@68: 		final long a=kmer&coreMask, b=rkmer&coreMask;
jpayne@68: 
jpayne@68: 		int max=-1, maxPos=0;
jpayne@68: 		if(a==b){
jpayne@68: 			return fillRightCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 		}else if(a>b){
jpayne@68: //			return fillRightCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 			final long core=a;
jpayne@68: 			final int way=kmerToWay(core);
jpayne@68: //			final AbstractKmerTable table=tables[way];
jpayne@68: 			final HashArray table=(HashArray) tables[way];
jpayne@68: 			final int cell=table.kmerToCell(kmer);
jpayne@68: 			for(int i=0; i<=3; i++){
jpayne@68: 				final long key=kmer|i;
jpayne@68: 				final int count=Tools.max(0, table.getValue(key, cell));
jpayne@68: //				final int count=Tools.max(0, table.getValue(key));
jpayne@68: 				counts[i]=count;
jpayne@68: 				if(count>max){
jpayne@68: 					max=count;
jpayne@68: 					maxPos=i;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}else{
jpayne@68: //			return fillRightCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 			//Use a higher increment for the left bits
jpayne@68: 			final long core=b;
jpayne@68: 			final int way=kmerToWay(core);
jpayne@68: //			final AbstractKmerTable table=tables[way];
jpayne@68: 			final HashArray table=(HashArray) tables[way];
jpayne@68: 			final int cell=table.kmerToCell(rkmer);
jpayne@68: 			final long incr=1L<<(2*(k-1));
jpayne@68: 			long j=3*incr;
jpayne@68: 			for(int i=0; i<=3; i++, j-=incr){
jpayne@68: 				final long key=rkmer|j;
jpayne@68: 				final int count=Tools.max(0, table.getValue(key, cell));
jpayne@68: //				final int count=Tools.max(0, table.getValue(key));
jpayne@68: 				counts[i]=count;
jpayne@68: 				if(count>max){
jpayne@68: 					max=count;
jpayne@68: 					maxPos=i;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return maxPos;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	//TODO: Change this to take advantage of coreMask
jpayne@68: 	//Requires special handling of core palindromes.
jpayne@68: 	//Thus it would be easiest to just handle odd kmers, and K is normally 31 anyway.
jpayne@68: 	public int fillRightCounts_safe(long kmer, long rkmer, int[] counts, long mask, int shift2){
jpayne@68: 		assert(kmer==rcomp(rkmer));
jpayne@68: 		if(verbose){outstream.println("fillRightCounts:   "+toText(kmer)+",   "+toText(rkmer));}
jpayne@68: 		kmer=(kmer<<2)&mask;
jpayne@68: 		rkmer=(rkmer>>>2);
jpayne@68: 		int max=-1, maxPos=0;
jpayne@68: 		
jpayne@68: 		for(int i=0; i<=3; i++){
jpayne@68: 			long kmer2=kmer|((long)i);
jpayne@68: 			long rkmer2=rkmer|(((long)AminoAcid.numberToComplement[i])<<shift2);
jpayne@68: 			if(verbose){outstream.println("kmer:               "+toText(kmer2)+", "+toText(rkmer2));}
jpayne@68: 			assert(kmer2==(kmer2&mask));
jpayne@68: 			assert(rkmer2==(rkmer2&mask));
jpayne@68: 			assert(kmer2==rcomp(rkmer2));
jpayne@68: 			long key=toValue(kmer2, rkmer2);
jpayne@68: 			int way=kmerToWay(key);
jpayne@68: 			int count=tables[way].getValue(key);
jpayne@68: 			assert(count==NOT_PRESENT || count>=0);
jpayne@68: 			count=Tools.max(count, 0);
jpayne@68: 			counts[i]=count;
jpayne@68: 			if(count>max){
jpayne@68: 				max=count;
jpayne@68: 				maxPos=i;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return maxPos;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** For KmerCompressor. */
jpayne@68: 	public int fillRightCountsRcompOnly(long kmer, long rkmer, int[] counts, long mask, int shift2){
jpayne@68: 		assert(kmer==rcomp(rkmer));
jpayne@68: 		if(verbose){outstream.println("fillRightCounts:   "+toText(kmer)+",   "+toText(rkmer));}
jpayne@68: 		kmer=(kmer<<2)&mask;
jpayne@68: 		rkmer=(rkmer>>>2);
jpayne@68: 		int max=-1, maxPos=0;
jpayne@68: 		
jpayne@68: 		for(int i=0; i<=3; i++){
jpayne@68: 			long kmer2=kmer|((long)i);
jpayne@68: 			long rkmer2=rkmer|(((long)AminoAcid.numberToComplement[i])<<shift2);
jpayne@68: 			if(verbose){outstream.println("kmer:               "+toText(kmer2)+", "+toText(rkmer2));}
jpayne@68: 			assert(kmer2==(kmer2&mask));
jpayne@68: 			assert(rkmer2==(rkmer2&mask));
jpayne@68: 			assert(kmer2==rcomp(rkmer2));
jpayne@68: 			long key=rkmer2;
jpayne@68: 			int way=kmerToWay(key);
jpayne@68: 			int count=tables[way].getValue(key);
jpayne@68: 			assert(count==NOT_PRESENT || count>=0);
jpayne@68: 			count=Tools.max(count, 0);
jpayne@68: 			counts[i]=count;
jpayne@68: 			if(count>max){
jpayne@68: 				max=count;
jpayne@68: 				maxPos=i;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return maxPos;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int fillLeftCounts(long kmer, long rkmer, int[] counts, long mask, int shift2){
jpayne@68: 		if(FAST_FILL && MASK_CORE && k>2/*((k&1)==1)*/){
jpayne@68: 			return fillLeftCounts_fast(kmer, rkmer, counts, mask, shift2);
jpayne@68: 		}else{
jpayne@68: 			return fillLeftCounts_safe(kmer, rkmer, counts, mask, shift2);
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int fillLeftCounts_fast(final long kmer0, final long rkmer0, int[] counts,
jpayne@68: 			long mask, int shift2){
jpayne@68: //		assert((k&1)==1) : k;
jpayne@68: 		assert(MASK_CORE);
jpayne@68: 		final long rkmer=(rkmer0<<2)&mask;
jpayne@68: 		final long kmer=(kmer0>>>2);
jpayne@68: 		final long a=kmer&coreMask, b=rkmer&coreMask;
jpayne@68: 		
jpayne@68: 		int max=-1, maxPos=0;
jpayne@68: 		if(a==b){
jpayne@68: 			return fillLeftCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 		}else if(a>b){
jpayne@68: //			return fillLeftCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 			
jpayne@68: 			final long core=a;
jpayne@68: 			final int way=kmerToWay(core);
jpayne@68: //			final AbstractKmerTable table=tables[way];
jpayne@68: 			final HashArray table=(HashArray) tables[way];
jpayne@68: 			final int cell=table.kmerToCell(kmer);
jpayne@68: 			final long incr=1L<<(2*(k-1));
jpayne@68: 			long j=0;//Must be declared as a long, outside of loop
jpayne@68: 			for(int i=0; i<=3; i++, j+=incr){
jpayne@68: //				final long rkmer2=rkmer|((long)AminoAcid.numberToComplement[i]);
jpayne@68: //				final long kmer2=kmer|(((long)i)<<shift2);
jpayne@68: //				final long key2=toValue(rkmer2, kmer2);
jpayne@68: //				int way2=kmerToWay(key2);
jpayne@68: //				int count2=tables[way2].getValue(key2);
jpayne@68: //				count2=Tools.max(count2, 0);
jpayne@68: 				
jpayne@68: 				final long key=kmer|j;
jpayne@68: 				final int count=Tools.max(0, table.getValue(key, cell));
jpayne@68: //				int count=Tools.max(0, table.getValue(key));
jpayne@68: 
jpayne@68: //				assert(way==way2) : i+", "+way+", "+way2;
jpayne@68: //				assert(key==key2) : i+", "+Long.toBinaryString(kmer)+", "+
jpayne@68: //					Long.toBinaryString(key)+", "+Long.toBinaryString(key2)+
jpayne@68: //					", "+Long.toBinaryString(j)+", "+Long.toBinaryString(incr);
jpayne@68: //				assert(count==count2) : i+", "+count+", "+count2;
jpayne@68: 				
jpayne@68: 				counts[i]=count;
jpayne@68: 				if(count>max){
jpayne@68: 					max=count;
jpayne@68: 					maxPos=i;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return maxPos;
jpayne@68: 		}else{
jpayne@68: //			return fillLeftCounts_safe(kmer0, rkmer0, counts, mask, shift2);
jpayne@68: 			final long core=b;
jpayne@68: 			final int way=kmerToWay(core);
jpayne@68: //			final AbstractKmerTable table=tables[way];
jpayne@68: 			final HashArray table=(HashArray) tables[way];
jpayne@68: 			final int cell=table.kmerToCell(rkmer);
jpayne@68: 			for(int i=0, j=3; i<=3; i++, j--){
jpayne@68: 				final long key=rkmer|j;
jpayne@68: 				final int count=Tools.max(0, table.getValue(key, cell));
jpayne@68: //				final int count=Tools.max(0, table.getValue(key));
jpayne@68: 				
jpayne@68: 				counts[i]=count;
jpayne@68: 				if(count>max){
jpayne@68: 					max=count;
jpayne@68: 					maxPos=i;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return maxPos;
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public int fillLeftCounts_safe(final long kmer0, final long rkmer0, int[] counts, long mask, int shift2){
jpayne@68: 		assert(kmer0==rcomp(rkmer0));
jpayne@68: 		if(verbose){outstream.println("fillLeftCounts:    "+toText(kmer0)+",   "+toText(rkmer0));}
jpayne@68: 		long rkmer=(rkmer0<<2)&mask;
jpayne@68: 		long kmer=(kmer0>>>2);
jpayne@68: 		int max=-1, maxPos=0;
jpayne@68: //		assert(false) : shift2+", "+k;
jpayne@68: 		for(int i=0; i<=3; i++){
jpayne@68: 			final long rkmer2=rkmer|((long)AminoAcid.numberToComplement[i]);
jpayne@68: 			final long kmer2=kmer|(((long)i)<<shift2);
jpayne@68: 			if(verbose){outstream.println("kmer:             "+toText(kmer2)+",     "+toText(rkmer2));}
jpayne@68: 			assert(kmer2==(kmer2&mask));
jpayne@68: 			assert(rkmer2==(rkmer2&mask));
jpayne@68: 			assert(kmer2==rcomp(rkmer2)) : "\n"+"kmer:      \t"+toText(rcomp(rkmer2))+", "+toText(rcomp(kmer2));
jpayne@68: 			long key=toValue(rkmer2, kmer2);
jpayne@68: 			int way=kmerToWay(key);
jpayne@68: 			int count=tables[way].getValue(key);
jpayne@68: 			assert(count==NOT_PRESENT || count>=0);
jpayne@68: 			count=Tools.max(count, 0);
jpayne@68: 			counts[i]=count;
jpayne@68: 			if(count>max){
jpayne@68: 				max=count;
jpayne@68: 				maxPos=i;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return maxPos;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------       Printing Methods       ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public boolean dumpKmersAsBytes(String fname, int mincount, int maxcount, boolean printTime, AtomicLong remaining){
jpayne@68: 		if(fname==null){return false;}
jpayne@68: 		Timer t=new Timer();
jpayne@68: 		
jpayne@68: 		ByteStreamWriter bsw=new ByteStreamWriter(fname, overwrite, append, true);
jpayne@68: 		bsw.start();
jpayne@68: 		for(AbstractKmerTable set : tables){
jpayne@68: 			set.dumpKmersAsBytes(bsw, k, mincount, maxcount, remaining);
jpayne@68: 		}
jpayne@68: 		bsw.poisonAndWait();
jpayne@68: 		
jpayne@68: 		t.stop();
jpayne@68: 		if(printTime){outstream.println("Kmer Dump Time:             \t"+t);}
jpayne@68: 		return bsw.errorState;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	@Override
jpayne@68: 	public boolean dumpKmersAsBytes_MT(String fname, int mincount, int maxcount, boolean printTime, AtomicLong remaining){
jpayne@68: 		final int threads=Tools.min(Shared.threads(), tables.length);
jpayne@68: 		if(threads<3 || DumpThread.NUM_THREADS==1){return dumpKmersAsBytes(fname, mincount, maxcount, printTime, remaining);}
jpayne@68: 		
jpayne@68: 		if(fname==null){return false;}
jpayne@68: 		Timer t=new Timer();
jpayne@68: 		
jpayne@68: 		ByteStreamWriter bsw=new ByteStreamWriter(fname, overwrite, append, true);
jpayne@68: 		bsw.start();
jpayne@68: 		DumpThread.dump(k, mincount, maxcount, tables, bsw, remaining);
jpayne@68: 		bsw.poisonAndWait();
jpayne@68: 		
jpayne@68: 		t.stop();
jpayne@68: 		if(printTime){outstream.println("Kmer Dump Time:             \t"+t);}
jpayne@68: 		return bsw.errorState;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Recall Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	private final long rcomp(long kmer){return AminoAcid.reverseComplementBinaryFast(kmer, k);}
jpayne@68: 	private final StringBuilder toText(long kmer){return AbstractKmerTable.toText(kmer, k);}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Static Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Transforms a kmer into a canonical value stored in the table.  Expected to be inlined.
jpayne@68: 	 * @param kmer Forward kmer
jpayne@68: 	 * @param rkmer Reverse kmer
jpayne@68: 	 * @return Canonical value
jpayne@68: 	 */
jpayne@68: 	public final long toValue(long kmer, long rkmer){
jpayne@68: //		long value=(rcomp ? Tools.max(kmer, rkmer) : kmer);
jpayne@68: //		return value;
jpayne@68: 		if(!rcomp){return kmer;}
jpayne@68: 		final long a=kmer&coreMask, b=rkmer&coreMask;
jpayne@68: 		return a>b ? kmer : a<b ? rkmer : Tools.max(kmer, rkmer);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------       Final Primitives       ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	@Override
jpayne@68: 	public int kbig(){return k;}
jpayne@68: 	@Override
jpayne@68: 	public long filterMemory(int pass){return ((pass&1)==0) ? filterMemory0 : filterMemory1;}
jpayne@68: 	@Override
jpayne@68: 	public boolean ecco(){return ecco;}
jpayne@68: 	@Override
jpayne@68: 	public boolean qtrimLeft(){return qtrimLeft;}
jpayne@68: 	@Override
jpayne@68: 	public boolean qtrimRight(){return qtrimRight;}
jpayne@68: 	@Override
jpayne@68: 	public float minAvgQuality(){return minAvgQuality;}
jpayne@68: 	@Override
jpayne@68: 	public long tableMemory(){return tableMemory;}
jpayne@68: 	@Override
jpayne@68: 	public long estimatedKmerCapacity(){return estimatedKmerCapacity;}
jpayne@68: 	@Override
jpayne@68: 	public int ways(){return ways;}
jpayne@68: 	@Override
jpayne@68: 	public boolean rcomp(){return rcomp;}
jpayne@68: 	
jpayne@68: 	public final int kmerToWay(final long kmer){
jpayne@68: 		final int way=(int)((kmer&coreMask)%ways);
jpayne@68: 		return way;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Hold kmers.  A kmer X such that X%WAYS=Y will be stored in tables[Y] */
jpayne@68: 	private AbstractKmerTable[] tables;
jpayne@68: 	
jpayne@68: 	public AbstractKmerTable[] tables(){return tables;}
jpayne@68: 	
jpayne@68: 	public long filterMemoryOverride=0;
jpayne@68: 	
jpayne@68: 	public final int tableType; //AbstractKmerTable.ARRAY1D;
jpayne@68: 	
jpayne@68: 	private final int bytesPerKmer;
jpayne@68: 
jpayne@68: 	private final long usableMemory;
jpayne@68: 	private final long filterMemory0;
jpayne@68: 	private final long filterMemory1;
jpayne@68: 	private final long tableMemory;
jpayne@68: 	private final long estimatedKmerCapacity;
jpayne@68: 	
jpayne@68: 	/** Number of tables (and threads, during loading) */
jpayne@68: 	private final boolean prealloc;
jpayne@68: 	
jpayne@68: 	/** Number of tables (and threads, during loading) */
jpayne@68: 	public final int ways;
jpayne@68: 	
jpayne@68: 	/** Normal kmer length */
jpayne@68: 	public final int k;
jpayne@68: 	/** k-1; used in some expressions */
jpayne@68: 	public final int k2;
jpayne@68: 	
jpayne@68: 	public final long coreMask;
jpayne@68: 	
jpayne@68: 	/** Look for reverse-complements as well as forward kmers.  Default: true */
jpayne@68: 	private final boolean rcomp;
jpayne@68: 	
jpayne@68: 	/** Quality-trim the left side */
jpayne@68: 	public final boolean qtrimLeft;
jpayne@68: 	/** Quality-trim the right side */
jpayne@68: 	public final boolean qtrimRight;
jpayne@68: 	/** Trim bases at this quality or below.  Default: 4 */
jpayne@68: 	public final float trimq;
jpayne@68: 	/** Error rate for trimming (derived from trimq) */
jpayne@68: 	private final float trimE;
jpayne@68: 	
jpayne@68: 	/** Throw away reads below this average quality after trimming.  Default: 0 */
jpayne@68: 	public final float minAvgQuality;
jpayne@68: 	/** If positive, calculate average quality from the first X bases only.  Default: 0 */
jpayne@68: 	public final int minAvgQualityBases;
jpayne@68: 	
jpayne@68: 	/** Ignore kmers with probability of correctness less than this */
jpayne@68: 	public final float minProb;
jpayne@68: 	
jpayne@68: 	/** Correct via overlap */
jpayne@68: 	private final boolean ecco;
jpayne@68: 	
jpayne@68: 	/** Attempt to merge via overlap prior to counting kmers */
jpayne@68: 	private final boolean merge;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Walker            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public Walker walk(){
jpayne@68: 		return new WalkerKST();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public class WalkerKST extends Walker {
jpayne@68: 		
jpayne@68: 		WalkerKST(){
jpayne@68: 			w=tables[0].walk();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** 
jpayne@68: 		 * Fills this object with the next key and value.
jpayne@68: 		 * @return True if successful.
jpayne@68: 		 */
jpayne@68: 		public boolean next(){
jpayne@68: 			if(w==null){return false;}
jpayne@68: 			if(w.next()){return true;}
jpayne@68: 			if(tnum<tables.length){tnum++;}
jpayne@68: 			w=(tnum<tables.length ? tables[tnum].walk() : null);
jpayne@68: 			return w==null ? false : w.next();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		public long kmer(){return w.kmer();}
jpayne@68: 		public int value(){return w.value();}
jpayne@68: 		
jpayne@68: 		private Walker w=null;
jpayne@68: 
jpayne@68: 		/** current table number */
jpayne@68: 		private int tnum=0;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: }