Mercurial > repos > rliterman > csp2
diff CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/fun/ProbShared2.java @ 68:5028fdace37b
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author | jpayne |
---|---|
date | Tue, 18 Mar 2025 16:23:26 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/fun/ProbShared2.java Tue Mar 18 16:23:26 2025 -0400 @@ -0,0 +1,75 @@ +package fun; + +import java.util.HashSet; +import java.util.Random; + +import dna.AminoAcid; + +public class ProbShared2 { + + public static void main(String args[]){ + int k=Integer.parseInt(args[0]); + int len1=Integer.parseInt(args[1]); + int len2=Integer.parseInt(args[2]); + int rounds=Integer.parseInt(args[3]); + + System.out.println("Probability: "+simulate(k, len1, len2, rounds)); + } + + static double simulate(int k, int len1, int len2, int rounds){ + int successes=0; + final HashSet<Long> set=new HashSet<Long>(); + for(int i=0; i<rounds; i++){ + successes+=simulateOnePair(k, len1, len2, set); + } + return successes/(double)rounds; + } + + static int simulateOnePair(int k, int len1, int len2, HashSet<Long> set){ + set.clear(); + + final int shift=2*k; + final long mask=(shift>63 ? -1L : ~((-1L)<<shift)); + long kmer=0; + int len=0; + + byte[] bases=randomSequence(len2); + + for(int i=0; i<bases.length; i++){ + byte b=bases[i]; + long x=baseToNumber[b]; + kmer=((kmer<<2)|x)&mask; + if(x<0){len=0;}else{len++;} + if(len>=k){ + set.add(kmer); + } + } + + bases=randomSequence(len1); + + for(int i=0; i<bases.length; i++){ + byte b=bases[i]; + long x=baseToNumber[b]; + kmer=((kmer<<2)|x)&mask; + if(x<0){len=0;}else{len++;} + if(len>=k){ + if(set.contains(kmer)){return 1;} + } + } + return 0; + } + + static byte[] randomSequence(int len){ + byte[] array=new byte[len]; + for(int i=0; i<len; i++){ + int number=randy.nextInt(4); + array[i]=AminoAcid.numberToBase[number]; + } + return array; + } + + static final Random randy=new Random(); + static final byte[] numberToBase=AminoAcid.numberToBase; + static final byte[] baseToNumber=AminoAcid.baseToNumber; + +}