Mercurial > repos > rliterman > csp2

comparison CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/fun/ProbShared2.java @ 68:5028fdace37b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
comparison
equal deleted inserted replaced
67:0e9998148a16 68:5028fdace37b
1 package fun;
2
3 import java.util.HashSet;
4 import java.util.Random;
5
6 import dna.AminoAcid;
7
8 public class ProbShared2 {
9
10 public static void main(String args[]){
11 int k=Integer.parseInt(args[0]);
12 int len1=Integer.parseInt(args[1]);
13 int len2=Integer.parseInt(args[2]);
14 int rounds=Integer.parseInt(args[3]);
15
16 System.out.println("Probability: "+simulate(k, len1, len2, rounds));
17 }
18
19 static double simulate(int k, int len1, int len2, int rounds){
20 int successes=0;
21 final HashSet<Long> set=new HashSet<Long>();
22 for(int i=0; i<rounds; i++){
23 successes+=simulateOnePair(k, len1, len2, set);
24 }
25 return successes/(double)rounds;
26 }
27
28 static int simulateOnePair(int k, int len1, int len2, HashSet<Long> set){
29 set.clear();
30
31 final int shift=2*k;
32 final long mask=(shift>63 ? -1L : ~((-1L)<<shift));
33 long kmer=0;
34 int len=0;
35
36 byte[] bases=randomSequence(len2);
37
38 for(int i=0; i<bases.length; i++){
39 byte b=bases[i];
40 long x=baseToNumber[b];
41 kmer=((kmer<<2)|x)&mask;
42 if(x<0){len=0;}else{len++;}
43 if(len>=k){
44 set.add(kmer);
45 }
46 }
47
48 bases=randomSequence(len1);
49
50 for(int i=0; i<bases.length; i++){
51 byte b=bases[i];
52 long x=baseToNumber[b];
53 kmer=((kmer<<2)|x)&mask;
54 if(x<0){len=0;}else{len++;}
55 if(len>=k){
56 if(set.contains(kmer)){return 1;}
57 }
58 }
59 return 0;
60 }
61
62 static byte[] randomSequence(int len){
63 byte[] array=new byte[len];
64 for(int i=0; i<len; i++){
65 int number=randy.nextInt(4);
66 array[i]=AminoAcid.numberToBase[number];
67 }
68 return array;
69 }
70
71 static final Random randy=new Random();
72 static final byte[] numberToBase=AminoAcid.numberToBase;
73 static final byte[] baseToNumber=AminoAcid.baseToNumber;
74
75 }