Mercurial > repos > rliterman > csp2
view 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 source
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; }