jpayne@68: package fun;
jpayne@68: 
jpayne@68: import java.util.HashSet;
jpayne@68: import java.util.Random;
jpayne@68: 
jpayne@68: import dna.AminoAcid;
jpayne@68: 
jpayne@68: public class ProbShared2 {
jpayne@68: 
jpayne@68: 	public static void main(String args[]){
jpayne@68: 		int k=Integer.parseInt(args[0]);
jpayne@68: 		int len1=Integer.parseInt(args[1]);
jpayne@68: 		int len2=Integer.parseInt(args[2]);
jpayne@68: 		int rounds=Integer.parseInt(args[3]);
jpayne@68: 		
jpayne@68: 		System.out.println("Probability:   "+simulate(k, len1, len2, rounds));
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static double simulate(int k, int len1, int len2, int rounds){
jpayne@68: 		int successes=0;
jpayne@68: 		final HashSet<Long> set=new HashSet<Long>();
jpayne@68: 		for(int i=0; i<rounds; i++){
jpayne@68: 			successes+=simulateOnePair(k, len1, len2, set);
jpayne@68: 		}
jpayne@68: 		return successes/(double)rounds;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static int simulateOnePair(int k, int len1, int len2, HashSet<Long> set){
jpayne@68: 		set.clear();
jpayne@68: 		
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: 		
jpayne@68: 		byte[] bases=randomSequence(len2);
jpayne@68: 		
jpayne@68: 		for(int i=0; i<bases.length; i++){
jpayne@68: 			byte b=bases[i];
jpayne@68: 			long x=baseToNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			if(x<0){len=0;}else{len++;}
jpayne@68: 			if(len>=k){
jpayne@68: 				set.add(kmer);
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		bases=randomSequence(len1);
jpayne@68: 		
jpayne@68: 		for(int i=0; i<bases.length; i++){
jpayne@68: 			byte b=bases[i];
jpayne@68: 			long x=baseToNumber[b];
jpayne@68: 			kmer=((kmer<<2)|x)&mask;
jpayne@68: 			if(x<0){len=0;}else{len++;}
jpayne@68: 			if(len>=k){
jpayne@68: 				if(set.contains(kmer)){return 1;}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		return 0;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static byte[] randomSequence(int len){
jpayne@68: 		byte[] array=new byte[len];
jpayne@68: 		for(int i=0; i<len; i++){
jpayne@68: 			int number=randy.nextInt(4);
jpayne@68: 			array[i]=AminoAcid.numberToBase[number];
jpayne@68: 		}
jpayne@68: 		return array;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	static final Random randy=new Random();
jpayne@68: 	static final byte[] numberToBase=AminoAcid.numberToBase;
jpayne@68: 	static final byte[] baseToNumber=AminoAcid.baseToNumber;
jpayne@68: 	
jpayne@68: }