jpayne@68: package kmer;
jpayne@68: 
jpayne@68: import java.io.PrintStream;
jpayne@68: import java.util.ArrayList;
jpayne@68: import java.util.Arrays;
jpayne@68: import java.util.concurrent.ArrayBlockingQueue;
jpayne@68: 
jpayne@68: import dna.AminoAcid;
jpayne@68: import fileIO.FileFormat;
jpayne@68: import fileIO.ReadWrite;
jpayne@68: import fileIO.TextStreamWriter;
jpayne@68: import jgi.BBMerge;
jpayne@68: import shared.PreParser;
jpayne@68: import shared.Shared;
jpayne@68: import shared.Timer;
jpayne@68: import shared.Tools;
jpayne@68: import stream.ConcurrentReadInputStream;
jpayne@68: import stream.Read;
jpayne@68: import structures.IntList;
jpayne@68: import structures.ListNum;
jpayne@68: 
jpayne@68: /**
jpayne@68:  * @author Brian Bushnell
jpayne@68:  * @date Mar 4, 2015
jpayne@68:  *
jpayne@68:  */
jpayne@68: public class TableLoaderLockFree {
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Initialization        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
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: 
jpayne@68: 		{//Preparse block for help, config files, and outstream
jpayne@68: 			PreParser pp=new PreParser(args, null, false);
jpayne@68: 			args=pp.args;
jpayne@68: 			outstream=pp.outstream;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		Timer t=new Timer();
jpayne@68: 		
jpayne@68: 		AbstractKmerTable[] tables=makeTables(AbstractKmerTable.ARRAY1D, 12, -1L, false, 1.0);
jpayne@68: 		
jpayne@68: 		int k=31;
jpayne@68: 		int mink=0;
jpayne@68: 		int speed=0;
jpayne@68: 		int hdist=0;
jpayne@68: 		int edist=0;
jpayne@68: 		boolean rcomp=true;
jpayne@68: 		boolean maskMiddle=false;
jpayne@68: 		
jpayne@68: 		//Create a new Loader instance
jpayne@68: 		TableLoaderLockFree loader=new TableLoaderLockFree(tables, k, mink, speed, hdist, edist, rcomp, maskMiddle);
jpayne@68: 		loader.setRefSkip(0);
jpayne@68: 		loader.hammingDistance2=0;
jpayne@68: 		loader.editDistance2=0;
jpayne@68: 		loader.storeMode(SET_IF_NOT_PRESENT);
jpayne@68: 		
jpayne@68: 		///And run it
jpayne@68: 		String[] refs=args;
jpayne@68: 		String[] literals=null;
jpayne@68: 		boolean keepNames=false;
jpayne@68: 		boolean useRefNames=false;
jpayne@68: 		long kmers=loader.processData(refs, literals, keepNames, useRefNames, false);
jpayne@68: 		t.stop();
jpayne@68: 
jpayne@68: 		outstream.println("Time:     \t"+t);
jpayne@68: 		outstream.println("Return:   \t"+kmers);
jpayne@68: 		outstream.println("refKmers: \t"+loader.refKmers);
jpayne@68: 		outstream.println("refBases: \t"+loader.refBases);
jpayne@68: 		outstream.println("refReads: \t"+loader.refReads);
jpayne@68: 		
jpayne@68: 		//Close the print stream if it was redirected
jpayne@68: 		Shared.closeStream(outstream);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public TableLoaderLockFree(AbstractKmerTable[] tables_, int k_){
jpayne@68: 		this(tables_, k_, 0, 0, 0, 0, true, false);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public TableLoaderLockFree(AbstractKmerTable[] tables_, int k_, int mink_, int speed_, int hdist_, int edist_, boolean rcomp_, boolean maskMiddle_){
jpayne@68: 		tables=tables_;
jpayne@68: 		k=k_;
jpayne@68: 		k2=k-1;
jpayne@68: 		mink=mink_;
jpayne@68: 		rcomp=rcomp_;
jpayne@68: 		useShortKmers=(mink>0 && mink<k);
jpayne@68: 		speed=speed_;
jpayne@68: 		hammingDistance=hdist_;
jpayne@68: 		editDistance=edist_;
jpayne@68: 		middleMask=maskMiddle ? ~(3L<<(2*(k/2))) : -1L;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Outer Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	
jpayne@68: //	public static AbstractKmerTable[] makeTables(int tableType, int initialSize, long coreMask, boolean growable){
jpayne@68: //		return AbstractKmerTable.preallocate(WAYS, tableType, initialSize, coreMask, growable);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: //	public static AbstractKmerTable[] makeTables(int tableType, int[] schedule, long coreMask){
jpayne@68: //		return AbstractKmerTable.preallocate(WAYS, tableType, schedule, coreMask);
jpayne@68: //	}
jpayne@68: 	
jpayne@68: 	public static AbstractKmerTable[] makeTables(int tableType, int bytesPerKmer, long coreMask, 
jpayne@68: 			boolean prealloc, double memRatio){
jpayne@68: 		ScheduleMaker scheduleMaker=new ScheduleMaker(WAYS, bytesPerKmer, prealloc, memRatio);
jpayne@68: 		int[] schedule=scheduleMaker.makeSchedule();
jpayne@68: 		return AbstractKmerTable.preallocate(WAYS, tableType, schedule, coreMask);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	ScheduleMaker scheduleMaker=new ScheduleMaker(WAYS, 12, false, 0.8);
jpayne@68: 	int[] schedule=scheduleMaker.makeSchedule();
jpayne@68: 	
jpayne@68: 	public long processData(String[] ref, String[] literal, boolean keepNames, boolean useRefNames, boolean ecc_){
jpayne@68: 		
jpayne@68: 		scaffoldNames=null;
jpayne@68: 		refNames=null;
jpayne@68: 		refScafCounts=null;
jpayne@68: 		scaffoldLengths=null;
jpayne@68: 		ecc=ecc_;
jpayne@68: 		
jpayne@68: 		if(keepNames){
jpayne@68: 			scaffoldNames=new ArrayList<String>();
jpayne@68: 			refNames=new ArrayList<String>();
jpayne@68: 			scaffoldLengths=new IntList();
jpayne@68: 			
jpayne@68: 			if(ref!=null){
jpayne@68: 				for(String s : ref){refNames.add(s);}
jpayne@68: 			}
jpayne@68: 			if(literal!=null){refNames.add("literal");}
jpayne@68: 			refScafCounts=new int[refNames.size()];
jpayne@68: 			
jpayne@68: 			if(useRefNames){toRefNames();}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		return spawnLoadThreads(ref, literal);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void setRefSkip(int x){setRefSkip(x, x);}
jpayne@68: 	
jpayne@68: 	public void setRefSkip(int min, int max){
jpayne@68: 		max=Tools.max(min, max);
jpayne@68: 		if(min==max){
jpayne@68: 			minRefSkip=maxRefSkip=min;
jpayne@68: 		}else{
jpayne@68: 			minRefSkip=min;
jpayne@68: 			maxRefSkip=max;
jpayne@68: 		}
jpayne@68: 		variableRefSkip=(minRefSkip!=maxRefSkip);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	public void storeMode(final int x){
jpayne@68: 		assert(x==SET_IF_NOT_PRESENT || x==SET_ALWAYS || x==INCREMENT);
jpayne@68: 		storeMode=x;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 
jpayne@68: 	/**
jpayne@68: 	 * Fills tables with kmers from references, using multiple LoadThread.
jpayne@68: 	 * @return Number of kmers stored.
jpayne@68: 	 */
jpayne@68: 	private long spawnLoadThreads(String[] ref, String[] literal){
jpayne@68: 		Timer t=new Timer();
jpayne@68: 		if((ref==null || ref.length<1) && (literal==null || literal.length<1)){return 0;}
jpayne@68: 		long added=0;
jpayne@68: 		
jpayne@68: 		/* Create load threads */
jpayne@68: 		LoadThread[] loaders=new LoadThread[WAYS];
jpayne@68: 		for(int i=0; i<loaders.length; i++){
jpayne@68: 			loaders[i]=new LoadThread(i);
jpayne@68: 			loaders[i].start();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* For each reference file... */
jpayne@68: 		int refNum=0;
jpayne@68: 		if(ref!=null){
jpayne@68: 			for(String refname : ref){
jpayne@68: 
jpayne@68: 				/* Start an input stream */
jpayne@68: 				FileFormat ff=FileFormat.testInput(refname, FileFormat.FASTA, null, false, true);
jpayne@68: 				ConcurrentReadInputStream cris=ConcurrentReadInputStream.getReadInputStream(-1L, false, ff, null, null, null, Shared.USE_MPI, true);
jpayne@68: 				cris.start(); //4567
jpayne@68: 				ListNum<Read> ln=cris.nextList();
jpayne@68: 				ArrayList<Read> reads=(ln!=null ? ln.list : null);
jpayne@68: 				
jpayne@68: 				/* Iterate through read lists from the input stream */
jpayne@68: 				while(ln!=null && reads!=null && reads.size()>0){//ln!=null prevents a compiler potential null access warning
jpayne@68: 					{
jpayne@68: 						/* Assign a unique ID number to each scaffold */
jpayne@68: 						ArrayList<Read> reads2=new ArrayList<Read>(reads);
jpayne@68: 						if(scaffoldNames!=null){
jpayne@68: 							for(Read r1 : reads2){
jpayne@68: 								final Read r2=r1.mate;
jpayne@68: 								final Integer id=scaffoldNames.size();
jpayne@68: 								if(ecc && r1!=null && r2!=null){BBMerge.findOverlapStrict(r1, r2, true);}
jpayne@68: 								refScafCounts[refNum]++;
jpayne@68: 								scaffoldNames.add(r1.id==null ? id.toString() : r1.id);
jpayne@68: 								int len=r1.length();
jpayne@68: 								r1.obj=id;
jpayne@68: 								if(r2!=null){
jpayne@68: 									r2.obj=id;
jpayne@68: 									len+=r2.length();
jpayne@68: 								}
jpayne@68: 								scaffoldLengths.add(len);
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 						
jpayne@68: 						if(REPLICATE_AMBIGUOUS){
jpayne@68: 							reads2=Tools.replicateAmbiguous(reads2, Tools.min(k, mink));
jpayne@68: 						}
jpayne@68: 
jpayne@68: 						/* Send a pointer to the read list to each LoadThread */
jpayne@68: 						for(LoadThread lt : loaders){
jpayne@68: 							boolean b=true;
jpayne@68: 							while(b){
jpayne@68: 								try {
jpayne@68: 									lt.queue.put(reads2);
jpayne@68: 									b=false;
jpayne@68: 								} catch (InterruptedException e) {
jpayne@68: 									//TODO:  This will hang due to still-running threads.
jpayne@68: 									throw new RuntimeException(e);
jpayne@68: 								}
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 
jpayne@68: 					/* Dispose of the old list and fetch a new one */
jpayne@68: 					cris.returnList(ln);
jpayne@68: 					ln=cris.nextList();
jpayne@68: 					reads=(ln!=null ? ln.list : null);
jpayne@68: 				}
jpayne@68: 				/* Cleanup */
jpayne@68: 				cris.returnList(ln);
jpayne@68: 				errorState|=ReadWrite.closeStream(cris);
jpayne@68: 				refNum++;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		/* If there are literal sequences to use as references */
jpayne@68: 		if(literal!=null){
jpayne@68: 			ArrayList<Read> list=new ArrayList<Read>(literal.length);
jpayne@68: 			if(verbose){outstream.println("Adding literals "+Arrays.toString(literal));}
jpayne@68: 
jpayne@68: 			/* Assign a unique ID number to each literal sequence */
jpayne@68: 			for(int i=0; i<literal.length; i++){
jpayne@68: 				if(scaffoldNames!=null){
jpayne@68: 					final Integer id=scaffoldNames.size();
jpayne@68: 					final Read r=new Read(literal[i].getBytes(), null, id);
jpayne@68: 					refScafCounts[refNum]++;
jpayne@68: 					scaffoldNames.add(id.toString());
jpayne@68: 					scaffoldLengths.add(r.length());
jpayne@68: 					r.obj=id;
jpayne@68: 				}else{
jpayne@68: 					final Read r=new Read(literal[i].getBytes(), null, i);
jpayne@68: 					list.add(r);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(REPLICATE_AMBIGUOUS){
jpayne@68: 				list=Tools.replicateAmbiguous(list, Tools.min(k, mink));
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			/* Send a pointer to the read list to each LoadThread */
jpayne@68: 			for(LoadThread lt : loaders){
jpayne@68: 				boolean b=true;
jpayne@68: 				while(b){
jpayne@68: 					try {
jpayne@68: 						lt.queue.put(list);
jpayne@68: 						b=false;
jpayne@68: 					} catch (InterruptedException e) {
jpayne@68: 						//TODO:  This will hang due to still-running threads.
jpayne@68: 						throw new RuntimeException(e);
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* Signal loaders to terminate */
jpayne@68: 		for(LoadThread lt : loaders){
jpayne@68: 			boolean b=true;
jpayne@68: 			while(b){
jpayne@68: 				try {
jpayne@68: 					lt.queue.put(POISON);
jpayne@68: 					b=false;
jpayne@68: 				} catch (InterruptedException e) {
jpayne@68: 					//TODO:  This will hang due to still-running threads.
jpayne@68: 					throw new RuntimeException(e);
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/* Wait for loaders to die, and gather statistics */
jpayne@68: 		for(LoadThread lt : loaders){
jpayne@68: 			while(lt.getState()!=Thread.State.TERMINATED){
jpayne@68: 				try {
jpayne@68: 					lt.join();
jpayne@68: 				} catch (InterruptedException e) {
jpayne@68: 					// TODO Auto-generated catch block
jpayne@68: 					e.printStackTrace();
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			added+=lt.addedT;
jpayne@68: 			refKmers+=lt.refKmersT;
jpayne@68: 			refBases+=lt.refBasesT;
jpayne@68: 			refReads+=lt.refReadsT;
jpayne@68: 		}
jpayne@68: 		//Correct statistics for number of threads, since each thread processes all reference data
jpayne@68: 		refKmers/=WAYS;
jpayne@68: 		refBases/=WAYS;
jpayne@68: 		refReads/=WAYS;
jpayne@68: 
jpayne@68: 		t.stop();
jpayne@68: 		if(DISPLAY_PROGRESS){
jpayne@68: 			outstream.println("Added "+added+" kmers; time: \t"+t);
jpayne@68: 			Shared.printMemory();
jpayne@68: 			outstream.println();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(verbose){
jpayne@68: 			TextStreamWriter tsw=new TextStreamWriter("stdout", false, false, false, FileFormat.TEXT);
jpayne@68: 			tsw.start();
jpayne@68: 			for(AbstractKmerTable table : tables){
jpayne@68: 				table.dumpKmersAsText(tsw, k, 1, Integer.MAX_VALUE);
jpayne@68: 			}
jpayne@68: 			tsw.poisonAndWait();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		return added;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Fills the scaffold names array with reference names.
jpayne@68: 	 */
jpayne@68: 	public void toRefNames(){
jpayne@68: 		final int numRefs=refNames.size();
jpayne@68: 		for(int r=0, s=1; r<numRefs; r++){
jpayne@68: 			final int scafs=refScafCounts[r];
jpayne@68: 			final int lim=s+scafs;
jpayne@68: 			final String name=ReadWrite.stripToCore(refNames.get(r));
jpayne@68: //			outstream.println("r="+r+", s="+s+", scafs="+scafs+", lim="+lim+", name="+name);
jpayne@68: 			while(s<lim){
jpayne@68: //				outstream.println(r+", "+s+". Setting "+scaffoldNames.get(s)+" -> "+name);
jpayne@68: 				scaffoldNames.set(s, name);
jpayne@68: 				s++;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Classes        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/**
jpayne@68: 	 * Loads kmers into a table.  Each thread handles all kmers X such that X%WAYS==tnum.
jpayne@68: 	 */
jpayne@68: 	private class LoadThread extends Thread{
jpayne@68: 		
jpayne@68: 		public LoadThread(final int tnum_){
jpayne@68: 			tnum=tnum_;
jpayne@68: 			map=tables[tnum];
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/**
jpayne@68: 		 * Get the next list of reads (or scaffolds) from the queue.
jpayne@68: 		 * @return List of reads
jpayne@68: 		 */
jpayne@68: 		private ArrayList<Read> fetch(){
jpayne@68: 			ArrayList<Read> list=null;
jpayne@68: 			while(list==null){
jpayne@68: 				try {
jpayne@68: 					list=queue.take();
jpayne@68: 				} catch (InterruptedException e) {
jpayne@68: 					// TODO Auto-generated catch block
jpayne@68: 					e.printStackTrace();
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return list;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		@Override
jpayne@68: 		public void run(){
jpayne@68: 			ArrayList<Read> reads=fetch();
jpayne@68: 			while(reads!=POISON){
jpayne@68: 				for(Read r1 : reads){
jpayne@68: 					assert(r1.pairnum()==0);
jpayne@68: 					final Read r2=r1.mate;
jpayne@68: 					
jpayne@68: 					addedT+=addToMap(r1, minRefSkip);
jpayne@68: 					if(r2!=null){addedT+=addToMap(r2, minRefSkip);}
jpayne@68: 				}
jpayne@68: 				reads=fetch();
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(map.canRebalance() && map.size()>2L*map.arrayLength()){
jpayne@68: 				map.rebalance();
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		/**
jpayne@68: 		 * @param r The current read to process
jpayne@68: 		 * @param skip Number of bases to skip between kmers
jpayne@68: 		 * @return Number of kmers stored
jpayne@68: 		 */
jpayne@68: 		private long addToMap(Read r, int skip){
jpayne@68: 			if(variableRefSkip){
jpayne@68: 				int rblen=r.length();
jpayne@68: 				skip=(rblen>20000000 ? k : rblen>5000000 ? 11 : rblen>500000 ? 2 : 0);
jpayne@68: 				skip=Tools.mid(minRefSkip, maxRefSkip, skip);
jpayne@68: 			}
jpayne@68: 			final byte[] bases=r.bases;
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: 			final long kmask=kMasks[k];
jpayne@68: 			long kmer=0;
jpayne@68: 			long rkmer=0;
jpayne@68: 			long added=0;
jpayne@68: 			int len=0;
jpayne@68: 			
jpayne@68: 			if(bases!=null){
jpayne@68: 				refReadsT++;
jpayne@68: 				refBasesT+=bases.length;
jpayne@68: 			}
jpayne@68: 			if(bases==null || bases.length<k){return 0;}
jpayne@68: 			
jpayne@68: 			final int id=(r.obj==null ? 1 : ((Integer)r.obj).intValue());
jpayne@68: 			
jpayne@68: 			if(skip>1){ //Process while skipping some kmers
jpayne@68: 				for(int i=0; i<bases.length; i++){
jpayne@68: 					byte b=bases[i];
jpayne@68: 					long x=AminoAcid.baseToNumber[b];
jpayne@68: 					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){len=0; rkmer=0;}else{len++;}
jpayne@68: 					if(verbose){outstream.println("Scanning1 i="+i+", kmer="+kmer+", rkmer="+rkmer+", bases="+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 					if(len>=k){
jpayne@68: 						refKmersT++;
jpayne@68: 						if(len%skip==0){
jpayne@68: 							final long extraBase=(i>=bases.length-1 ? -1 : AminoAcid.baseToNumber[bases[i+1]]);
jpayne@68: 							added+=addToMap(kmer, rkmer, k, extraBase, id, kmask, hammingDistance, editDistance);
jpayne@68: 							if(useShortKmers){
jpayne@68: 								if(i==k2){added+=addToMapRightShift(kmer, rkmer, id);}
jpayne@68: 								if(i==bases.length-1){added+=addToMapLeftShift(kmer, rkmer, extraBase, id);}
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}else{ //Process all kmers
jpayne@68: 				for(int i=0; i<bases.length; i++){
jpayne@68: 					byte b=bases[i];
jpayne@68: 					long x=AminoAcid.baseToNumber[b];
jpayne@68: 					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){len=0; rkmer=0;}else{len++;}
jpayne@68: 					if(verbose){outstream.println("Scanning2 i="+i+", kmer="+kmer+", rkmer="+rkmer+", bases="+new String(bases, Tools.max(0, i-k2), Tools.min(i+1, k)));}
jpayne@68: 					if(len>=k){
jpayne@68: 						refKmersT++;
jpayne@68: 						final long extraBase=(i>=bases.length-1 ? -1 : AminoAcid.baseToNumber[bases[i+1]]);
jpayne@68: 						final long atm=addToMap(kmer, rkmer, k, extraBase, id, kmask, hammingDistance, editDistance);
jpayne@68: 						added+=atm;
jpayne@68: //						assert(false) : atm+", "+map.contains(toValue(kmer, rkmer, kmask));
jpayne@68: 						if(useShortKmers){
jpayne@68: 							if(i==k2){added+=addToMapRightShift(kmer, rkmer, id);}
jpayne@68: 							if(i==bases.length-1){added+=addToMapLeftShift(kmer, rkmer, extraBase, id);}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return added;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 
jpayne@68: 		/**
jpayne@68: 		 * Adds short kmers on the left end of the read.
jpayne@68: 		 * @param kmer Forward kmer
jpayne@68: 		 * @param rkmer Reverse kmer
jpayne@68: 		 * @param extraBase Base added to end in case of deletions
jpayne@68: 		 * @param id Scaffold number
jpayne@68: 		 * @return Number of kmers stored
jpayne@68: 		 */
jpayne@68: 		private long addToMapLeftShift(long kmer, long rkmer, final long extraBase, final int id){
jpayne@68: 			if(verbose){outstream.println("addToMapLeftShift");}
jpayne@68: 			long added=0;
jpayne@68: 			for(int i=k-1; i>=mink; i--){
jpayne@68: 				kmer=kmer&rightMasks[i];
jpayne@68: 				rkmer=rkmer>>>2;
jpayne@68: 				long x=addToMap(kmer, rkmer, i, extraBase, id, kMasks[i], hammingDistance2, editDistance2);
jpayne@68: 				added+=x;
jpayne@68: 				if(verbose){
jpayne@68: 					if((toValue(kmer, rkmer, kMasks[i]))%WAYS==tnum){
jpayne@68: 						outstream.println("added="+x+"; i="+i+"; tnum="+tnum+"; Added left-shift kmer "+AminoAcid.kmerToString(kmer&~kMasks[i], i)+"; value="+(toValue(kmer, rkmer, kMasks[i]))+"; kmer="+kmer+"; rkmer="+rkmer+"; kmask="+kMasks[i]+"; rightMasks[i+1]="+rightMasks[i+1]);
jpayne@68: 						outstream.println("i="+i+"; tnum="+tnum+"; Looking for left-shift kmer "+AminoAcid.kmerToString(kmer&~kMasks[i], i));
jpayne@68: 						final long value=toValue(kmer, rkmer, kMasks[i]);
jpayne@68: 						if(map.contains(value)){outstream.println("Found "+value);}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return added;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 
jpayne@68: 		/**
jpayne@68: 		 * Adds short kmers on the right end of the read.
jpayne@68: 		 * @param kmer Forward kmer
jpayne@68: 		 * @param rkmer Reverse kmer
jpayne@68: 		 * @param id Scaffold number
jpayne@68: 		 * @return Number of kmers stored
jpayne@68: 		 */
jpayne@68: 		private long addToMapRightShift(long kmer, long rkmer, final int id){
jpayne@68: 			if(verbose){outstream.println("addToMapRightShift");}
jpayne@68: 			long added=0;
jpayne@68: 			for(int i=k-1; i>=mink; i--){
jpayne@68: 				long extraBase=kmer&3L;
jpayne@68: 				kmer=kmer>>>2;
jpayne@68: 				rkmer=rkmer&rightMasks[i];
jpayne@68: //				assert(Long.numberOfLeadingZeros(kmer)>=2*(32-i)) : Long.numberOfLeadingZeros(kmer)+", "+i+", "+kmer+", "+kMasks[i];
jpayne@68: //				assert(Long.numberOfLeadingZeros(rkmer)>=2*(32-i)) : Long.numberOfLeadingZeros(rkmer)+", "+i+", "+rkmer+", "+kMasks[i];
jpayne@68: 				long x=addToMap(kmer, rkmer, i, extraBase, id, kMasks[i], hammingDistance2, editDistance2);
jpayne@68: 				added+=x;
jpayne@68: 				if(verbose){
jpayne@68: 					if((toValue(kmer, rkmer, kMasks[i]))%WAYS==tnum){
jpayne@68: 						outstream.println("added="+x+"; i="+i+"; tnum="+tnum+"; Added right-shift kmer "+AminoAcid.kmerToString(kmer&~kMasks[i], i)+"; value="+(toValue(kmer, rkmer, kMasks[i]))+"; kmer="+kmer+"; rkmer="+rkmer+"; kmask="+kMasks[i]+"; rightMasks[i+1]="+rightMasks[i+1]);
jpayne@68: 						outstream.println("i="+i+"; tnum="+tnum+"; Looking for right-shift kmer "+AminoAcid.kmerToString(kmer&~kMasks[i], i));
jpayne@68: 						final long value=toValue(kmer, rkmer, kMasks[i]);
jpayne@68: 						if(map.contains(value)){outstream.println("Found "+value);}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			return added;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		
jpayne@68: 		/**
jpayne@68: 		 * Adds this kmer to the table, including any mutations implied by editDistance or hammingDistance.
jpayne@68: 		 * @param kmer Forward kmer
jpayne@68: 		 * @param rkmer Reverse kmer
jpayne@68: 		 * @param len Kmer length
jpayne@68: 		 * @param extraBase Base added to end in case of deletions
jpayne@68: 		 * @param id Scaffold number
jpayne@68: 		 * @param kmask0
jpayne@68: 		 * @return Number of kmers stored
jpayne@68: 		 */
jpayne@68: 		private long addToMap(final long kmer, final long rkmer, final int len, final long extraBase, final int id, final long kmask0, final int hdist, final int edist){
jpayne@68: 			
jpayne@68: 			assert(kmask0==kMasks[len]) : kmask0+", "+len+", "+kMasks[len]+", "+Long.numberOfTrailingZeros(kmask0)+", "+Long.numberOfTrailingZeros(kMasks[len]);
jpayne@68: 			
jpayne@68: 			if(verbose){outstream.println("addToMap_A; len="+len+"; kMasks[len]="+kMasks[len]);}
jpayne@68: 			assert((kmer&kmask0)==0);
jpayne@68: 			final long added;
jpayne@68: 			if(hdist==0){
jpayne@68: 				final long key=toValue(kmer, rkmer, kmask0);
jpayne@68: 				if(failsSpeed(key)){return 0;}
jpayne@68: 				if(key%WAYS!=tnum){return 0;}
jpayne@68: 				if(verbose){outstream.println("addToMap_B: "+AminoAcid.kmerToString(kmer&~kMasks[len], len)+" = "+key);}
jpayne@68: 				if(storeMode==SET_IF_NOT_PRESENT){
jpayne@68: 					added=map.setIfNotPresent(key, id);
jpayne@68: 				}else if(storeMode==SET_ALWAYS){
jpayne@68: 					added=map.set(key, id);
jpayne@68: 				}else{
jpayne@68: 					assert(storeMode==INCREMENT);
jpayne@68: 					added=map.increment(key, 1);
jpayne@68: 				}
jpayne@68: 			}else if(edist>0){
jpayne@68: //				long extraBase=(i>=bases.length-1 ? -1 : AminoAcid.baseToNumber[bases[i+1]]);
jpayne@68: 				added=mutate(kmer, rkmer, len, id, edist, extraBase);
jpayne@68: 			}else{
jpayne@68: 				added=mutate(kmer, rkmer, len, id, hdist, -1);
jpayne@68: 			}
jpayne@68: 			if(verbose){outstream.println("addToMap added "+added+" keys.");}
jpayne@68: 			return added;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/**
jpayne@68: 		 * Mutate and store this kmer through 'dist' recursions.
jpayne@68: 		 * @param kmer Forward kmer
jpayne@68: 		 * @param rkmer Reverse kmer
jpayne@68: 		 * @param id Scaffold number
jpayne@68: 		 * @param dist Number of mutations
jpayne@68: 		 * @param extraBase Base added to end in case of deletions
jpayne@68: 		 * @return Number of kmers stored
jpayne@68: 		 */
jpayne@68: 		private long mutate(final long kmer, final long rkmer, final int len, final int id, final int dist, final long extraBase){
jpayne@68: 			long added=0;
jpayne@68: 			
jpayne@68: 			final long key=toValue(kmer, rkmer, kMasks[len]);
jpayne@68: 			
jpayne@68: 			if(verbose){outstream.println("mutate_A; len="+len+"; kmer="+kmer+"; rkmer="+rkmer+"; kMasks[len]="+kMasks[len]);}
jpayne@68: 			if(key%WAYS==tnum){
jpayne@68: 				if(verbose){outstream.println("mutate_B: "+AminoAcid.kmerToString(kmer&~kMasks[len], len)+" = "+key);}
jpayne@68: 				int x;
jpayne@68: 				if(storeMode==SET_IF_NOT_PRESENT){
jpayne@68: 					x=map.setIfNotPresent(key, id);
jpayne@68: 				}else if(storeMode==SET_ALWAYS){
jpayne@68: 					x=map.set(key, id);
jpayne@68: 				}else{
jpayne@68: 					assert(storeMode==INCREMENT);
jpayne@68: 					x=map.increment(key, 1);
jpayne@68: 					x=(x==1 ? 1 : 0);
jpayne@68: 				}
jpayne@68: 				if(verbose){outstream.println("mutate_B added "+x+" keys.");}
jpayne@68: 				added+=x;
jpayne@68: 				assert(map.contains(key));
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(dist>0){
jpayne@68: 				final int dist2=dist-1;
jpayne@68: 				
jpayne@68: 				//Sub
jpayne@68: 				for(int j=0; j<4; j++){
jpayne@68: 					for(int i=0; i<len; i++){
jpayne@68: 						final long temp=(kmer&clearMasks[i])|setMasks[j][i];
jpayne@68: 						if(temp!=kmer){
jpayne@68: 							long rtemp=AminoAcid.reverseComplementBinaryFast(temp, len);
jpayne@68: 							added+=mutate(temp, rtemp, len, id, dist2, extraBase);
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 				
jpayne@68: 				if(editDistance>0){
jpayne@68: 					//Del
jpayne@68: 					if(extraBase>=0 && extraBase<=3){
jpayne@68: 						for(int i=1; i<len; i++){
jpayne@68: 							final long temp=(kmer&leftMasks[i])|((kmer<<2)&rightMasks[i])|extraBase;
jpayne@68: 							if(temp!=kmer){
jpayne@68: 								long rtemp=AminoAcid.reverseComplementBinaryFast(temp, len);
jpayne@68: 								added+=mutate(temp, rtemp, len, id, dist2, -1);
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 
jpayne@68: 					//Ins
jpayne@68: 					final long eb2=kmer&3;
jpayne@68: 					for(int i=1; i<len; i++){
jpayne@68: 						final long temp0=(kmer&leftMasks[i])|((kmer&rightMasks[i])>>2);
jpayne@68: 						for(int j=0; j<4; j++){
jpayne@68: 							final long temp=temp0|setMasks[j][i-1];
jpayne@68: 							if(temp!=kmer){
jpayne@68: 								long rtemp=AminoAcid.reverseComplementBinaryFast(temp, len);
jpayne@68: 								added+=mutate(temp, rtemp, len, id, dist2, eb2);
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			return added;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/*--------------------------------------------------------------*/
jpayne@68: 		
jpayne@68: 		/** Number of kmers stored by this thread */
jpayne@68: 		public long addedT=0;
jpayne@68: 		/** Number of items encountered by this thread */
jpayne@68: 		public long refKmersT=0, refReadsT=0, refBasesT=0;
jpayne@68: 		/** Thread number; used to determine which kmers to store */
jpayne@68: 		public final int tnum;
jpayne@68: 		/** Buffer of input read lists */
jpayne@68: 		public final ArrayBlockingQueue<ArrayList<Read>> queue=new ArrayBlockingQueue<ArrayList<Read>>(32);
jpayne@68: 		
jpayne@68: 		/** Destination for storing kmers */
jpayne@68: 		private final AbstractKmerTable map;
jpayne@68: 		
jpayne@68: 	}
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: 	 * @param lengthMask Bitmask with single '1' set to left of kmer
jpayne@68: 	 * @return Canonical value
jpayne@68: 	 */
jpayne@68: 	private final long toValue(long kmer, long rkmer, long lengthMask){
jpayne@68: 		assert(lengthMask==0 || (kmer<lengthMask && rkmer<lengthMask)) : lengthMask+", "+kmer+", "+rkmer;
jpayne@68: 		long value=(rcomp ? Tools.max(kmer, rkmer) : kmer);
jpayne@68: 		return (value&middleMask)|lengthMask;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	final boolean passesSpeed(long key){
jpayne@68: 		return speed<1 || ((key&Long.MAX_VALUE)%17)>=speed;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	final boolean failsSpeed(long key){
jpayne@68: 		return speed>0 && ((key&Long.MAX_VALUE)%17)<speed;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Fields            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Has this class encountered errors while processing? */
jpayne@68: 	public boolean errorState=false;
jpayne@68: 	
jpayne@68: 	/** How to associate values with kmers */
jpayne@68: 	private int storeMode=SET_IF_NOT_PRESENT;
jpayne@68: 	
jpayne@68: 	/** Hold kmers.  A kmer X such that X%WAYS=Y will be stored in keySets[Y] */
jpayne@68: 	public AbstractKmerTable[] tables;
jpayne@68: 	/** A scaffold's name is stored at scaffoldNames.get(id).
jpayne@68: 	 * scaffoldNames[0] is reserved, so the first id is 1. */
jpayne@68: 	public ArrayList<String> scaffoldNames;
jpayne@68: 	/** Names of reference files (refNames[0] is valid). */
jpayne@68: 	public ArrayList<String> refNames;
jpayne@68: 	/** Number of scaffolds per reference. */
jpayne@68: 	public int[] refScafCounts;
jpayne@68: 	/** scaffoldLengths[id] stores the length of that scaffold */
jpayne@68: 	public IntList scaffoldLengths=new IntList();
jpayne@68: 	
jpayne@68: 	/** Make the middle base in a kmer a wildcard to improve sensitivity */
jpayne@68: 	public final boolean maskMiddle=false;
jpayne@68: 	/** Correct errors via read overlap */
jpayne@68: 	public boolean ecc=false;
jpayne@68: 	
jpayne@68: 	/** Store reference kmers with up to this many substitutions */
jpayne@68: 	public final int hammingDistance;
jpayne@68: 	/** Store reference kmers with up to this many edits (including indels) */
jpayne@68: 	public final int editDistance;
jpayne@68: 	/** Store short reference kmers with up to this many substitutions */
jpayne@68: 	public int hammingDistance2=-1;
jpayne@68: 	/** Store short reference kmers with up to this many edits (including indels) */
jpayne@68: 	public int editDistance2=-1;
jpayne@68: 	/** Always skip at least this many consecutive kmers when hashing reference.
jpayne@68: 	 * 1 means every kmer is used, 2 means every other, etc. */
jpayne@68: 	private int minRefSkip=0;
jpayne@68: 	/** Never skip more than this many consecutive kmers when hashing reference. */
jpayne@68: 	private int maxRefSkip=0;
jpayne@68: 	
jpayne@68: 	private boolean variableRefSkip=false;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------          Statistics          ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	long refReads=0;
jpayne@68: 	long refBases=0;
jpayne@68: 	long refKmers=0;
jpayne@68: 	
jpayne@68: 	long storedKmers=0;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------       Final Primitives       ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Look for reverse-complements as well as forward kmers.  Default: true */
jpayne@68: 	private final boolean rcomp;
jpayne@68: 	/** AND bitmask with 0's at the middle base */
jpayne@68: 	private final long middleMask;
jpayne@68: 	
jpayne@68: 	/** Normal kmer length */
jpayne@68: 	private final int k;
jpayne@68: 	/** k-1; used in some expressions */
jpayne@68: 	private final int k2;
jpayne@68: 	/** Shortest kmer to use for trimming */
jpayne@68: 	private final int mink;
jpayne@68: 	/** Attempt to match kmers shorter than normal k on read ends when doing kTrimming. */
jpayne@68: 	private final boolean useShortKmers;
jpayne@68: 	
jpayne@68: 	/** Fraction of kmers to skip, 0 to 16 out of 17 */
jpayne@68: 	private final int speed;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Static Fields        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Default initial size of data structures */
jpayne@68: 	private static final int initialSizeDefault=128000;
jpayne@68: 	
jpayne@68: 	/** Number of tables (and threads, during loading) */
jpayne@68: 	private static final int WAYS=7; //123
jpayne@68: 	/** Verbose messages */
jpayne@68: 	public static final boolean verbose=false; //123
jpayne@68: 	
jpayne@68: 	/** Print messages to this stream */
jpayne@68: 	private static PrintStream outstream=System.err;
jpayne@68: 	/** Display progress messages such as memory usage */
jpayne@68: 	public static boolean DISPLAY_PROGRESS=true;
jpayne@68: 	/** Indicates end of input stream */
jpayne@68: 	private static final ArrayList<Read> POISON=new ArrayList<Read>(0);
jpayne@68: 	/** Make unambiguous copies of ref sequences with ambiguous bases */
jpayne@68: 	public static boolean REPLICATE_AMBIGUOUS=false;
jpayne@68: 	
jpayne@68: 	/** x&clearMasks[i] will clear base i */
jpayne@68: 	private static final long[] clearMasks;
jpayne@68: 	/** x|setMasks[i][j] will set base i to j */
jpayne@68: 	private static final long[][] setMasks;
jpayne@68: 	/** x&leftMasks[i] will clear all bases to the right of i (exclusive) */
jpayne@68: 	private static final long[] leftMasks;
jpayne@68: 	/** x&rightMasks[i] will clear all bases to the left of i (inclusive) */
jpayne@68: 	private static final long[] rightMasks;
jpayne@68: 	/** x|kMasks[i] will set the bit to the left of the leftmost base */
jpayne@68: 	private static final long[] kMasks;
jpayne@68: 	
jpayne@68: 	public static final int SET_IF_NOT_PRESENT=1, SET_ALWAYS=2, INCREMENT=3;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------      Static Initializers     ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	static{
jpayne@68: 		clearMasks=new long[32];
jpayne@68: 		leftMasks=new long[32];
jpayne@68: 		rightMasks=new long[32];
jpayne@68: 		kMasks=new long[32];
jpayne@68: 		setMasks=new long[4][32];
jpayne@68: 		for(int i=0; i<32; i++){
jpayne@68: 			clearMasks[i]=~(3L<<(2*i));
jpayne@68: 		}
jpayne@68: 		for(int i=0; i<32; i++){
jpayne@68: 			leftMasks[i]=((-1L)<<(2*i));
jpayne@68: 		}
jpayne@68: 		for(int i=0; i<32; i++){
jpayne@68: 			rightMasks[i]=~((-1L)<<(2*i));
jpayne@68: 		}
jpayne@68: 		for(int i=0; i<32; i++){
jpayne@68: 			kMasks[i]=((1L)<<(2*i));
jpayne@68: 		}
jpayne@68: 		for(int i=0; i<32; i++){
jpayne@68: 			for(long j=0; j<4; j++){
jpayne@68: 				setMasks[(int)j][i]=(j<<(2*i));
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: }