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;
	
}