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