import cern.jet.random.Poisson; //导入依赖的package包/类
/**
   * Creates a new exponential number generator.
   * 
   * @param mean the mean for the distribution.
   */
  public PoissonDistr(double mean) {
engine = new MersenneTwister(new Date());
poisson = new Poisson(mean, engine);

//always sleep for some milliseconds in order not to have same seed for iterative PoissonDistr contruction
try {
	TimeUnit.MILLISECONDS.sleep(10);
} catch (InterruptedException e) {
   	SimLogger.printLine("impossible is occured! Poisson random number cannot be created!");
	e.printStackTrace();
   	System.exit(0);
}
  }
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * This resets the generators
 * @param seed
 */
public static void setSeed(int seed) {
    uniform = new Random(seed);
    poisson = new Poisson(1.0, new MersenneTwister(uniform.nextInt()));
    gamma = new Gamma(1.0, 5.0, new MersenneTwister(uniform.nextInt()));
    gaussian = new Normal(0.0, 1.0, new MersenneTwister(uniform.nextInt()));
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
public EnrichmentSignificance(EventsConfig con, ExperimentManager exptman, BindingManager bman, double minFoldChange, double genomeLength){
	this.config = con;
	this.manager = exptman;
	this.bindingManager = bman;
	this.minFoldChange = minFoldChange;
	this.genomeLength = genomeLength;
	binomial = new Binomial(100, .5, new DRand());
	poisson = new Poisson(1, new DRand());
	chisquare = new ChiSquare(1, new DRand());
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * NegativeBinomial (r p parameterization)
 * Parameterized as the distribution of the number of failures (X) before the rth success in independent trials
 * @param r : Number of successes 
 * @param p : Probability of success in each trial
 */
public NegativeBinomialDistrib(double r, double p){
	this.r = r;
	this.p = p;
	gamma = new Gamma(1, 1.0, new DRand());
	poisson = new Poisson(10, new DRand());
	nb = new NegativeBinomialDist(r, p);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
protected int calcCountThreshold(){
	int countThres=0;
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
	lambda = (totalReads*binWidth)/(regionLength*mappableRegion); 
	P.setMean(lambda);
	double l=1;
	for(int b=1; l>confThreshold; b++){
		l=1-P.cdf(b);
		countThres=b;
	}
	return(Math.max(1,countThres));
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
public PointEnrichmentTester(String base, GenomeConfig gcon,List<Point> g, List<Region> r, int numTest){
	outbase=base;
	gconfig=gcon;
	gff=g;
	regions=r;	
	poisson = new Poisson(1, new DRand());
	numItr = numTest;
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * NegativeBinomial (r p parameterization)
 * @param r : Number of successes 
 * @param p : Probability of success in each trial
 */
public NegativeBinomialDistrib(double r, double p){
	this.r = r;
	this.p = p;
	gamma = new Gamma(1, 1.0, new DRand());
	poisson = new Poisson(10, new DRand());
	nb = new NegativeBinomialDist(r, p);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
public EnrichmentSignificance(Config con, ExperimentSet exptSet, List<BindingEvent> features, double minFoldChange, double genomeLength){
	this.config = con;
	this.exptSet = exptSet;
	this.features = features;
	this.minFoldChange = minFoldChange;
	this.genomeLength = genomeLength;
	binomial = new Binomial(100, .5, new DRand());
	poisson = new Poisson(1, new DRand());
	chisquare = new ChiSquare(1, new DRand());
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
protected int calcCountThreshold(){
	int countThres=0;
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
	lambda = (totalReads*(readLength/binStep + binWidth/binStep))/(regionLength*mappableRegion/binStep); 
	P.setMean(lambda);
	double l=1;
	for(int b=1; l>confThreshold; b++){
		l=1-P.cdf(b);
		countThres=b;
	}
	return(Math.max(1,countThres));
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
private int calcHitCount_per_BP(double totalReads, double threshold){
	int countThres=0;
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
	double lambda = totalReads/config.mappable_genome_length;
	P.setMean(lambda);
	double l=1;
	for(int b=1; l>threshold; b++){
		l=1-P.cdf(b);	//p-value as the tail of Poisson
		countThres=b;
	}
	return Math.max(1,countThres);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
protected int calcHitCount_per_BP(double totalReads, double threshold){
	int countThres=0;
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
	double lambda = totalReads/mappable_genome_length;
	P.setMean(lambda);
	double l=1;
	for(int b=1; l>threshold; b++){
		l=1-P.cdf(b);	//p-value as the tail of Poisson
		countThres=b;
	}
	return Math.max(1,countThres);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
private int calcExpectedHitCount(double totalReads, double threshold, int regionWidth){
	int countThres=0;
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
	double lambda = totalReads*regionWidth /config.mappable_genome_length;
	P.setMean(lambda);
	double l=1;
	for(int b=1; l>threshold; b++){
		l=1-P.cdf(b);	//p-value as the tail of Poisson
		countThres=b;
	}
	return Math.max(1,countThres);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * Reset duplicate reads that pass Poisson threshold. 
 * The Poisson lambda parameter is calculated by an Gaussian average
 * that puts more weight for nearby bases (same chrom, same strand)
 */
private void capPerBaseCountWithPoissonGaussianFilter(double threshold, int width){
       double g[] = new double[width*4+1];
	NormalDistribution gaussianDist = new NormalDistribution(0, width*width);
	for (int i=0;i<g.length;i++)
		g[i]=gaussianDist.calcProbability((double)i);
	
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
		
	for(int i = 0; i < fivePrimeCounts.length; i++)
		for(int j = 0; j < fivePrimeCounts[i].length; j++){
			float counts[] = fivePrimeCounts[i][j];
			int pos[] = fivePrimePos[i][j]; 
			if(counts!=null){
				for(int k = 0; k < counts.length; k++){
					int posK = pos[k]; 
					double sum = 0;
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k+x>=counts.length|| pos[k+x]-posK>width*4)
							break;
						sum += counts[k+x]*g[pos[k+x]-posK];
					}
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k-x<0 || posK-pos[k-x]>width*4)
							break;
						sum += counts[k-x]*g[posK-pos[k-x]];
					}
					sum = sum/(1-g[0]);				// exclude this position for evaluation
					
					double countThres=0;
					P.setMean(sum);
					double pvalue=1;
					for(int b=1; pvalue>threshold; b++){
						pvalue=1-P.cdf(b);	//p-value as the tail of Poisson
						countThres=b;
					}
					if (counts[k] > Math.max(1,countThres))
						counts[k] = (float) Math.max(1,countThres);					
				}
			}
		}
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * Fit a truncated Poisson to the contents of a histogram.
 * Returns the background proportion in the sample
 *
 */
public double fitPoisson(RealValuedHistogram h, Sample samp){
	DRand re = new DRand();
	//Heuristic to find the upper bound for the truncated Poisson
	int pUpper = poissUpperBound;
	double uniformMean = samp.getHitCount()/(gen.getGenomeLength()/binWidth);
	if(cdfPercOfUniform>0 && cdfPercOfUniform<=1){
		Poisson uniPoiss = new Poisson(uniformMean, re);
		double tmpProp=0;
		int i=0;
		while(tmpProp<cdfPercOfUniform){
			tmpProp = uniPoiss.cdf(i);
			i++;
		}
		pUpper=Math.max(i,poissUpperBoundMin);
	}
	System.out.println("Truncated Poisson Upper Bound:\t"+pUpper);
	
	//Fit the Poisson
	int left=0, right=pUpper;
	double xsum=0, xcount=0;
	for(double i=left; i<=right; i++){
		xsum += i*h.getBin( h.getBinContainingVal(i));
		xcount += h.getBin( h.getBinContainingVal(i));
	}
	double xavg = xsum/xcount;
	UnivariateFunction func = new truncPoisson(xavg, left, right);
	double relativeAccuracy = 1.0e-6;
	double absoluteAccuracy = 1.0e-4;
	UnivariateOptimizer solver = new BrentOptimizer(relativeAccuracy, absoluteAccuracy);
	UnivariatePointValuePair pvp = solver.optimize(100, func, GoalType.MINIMIZE, 0.001, 50.0, xavg);
	double lambda = pvp.getPoint();
	System.out.println("xavg: "+ xavg+"\tlambda: "+lambda);
	
	//Calculate the background proportion
	Poisson poiss = new Poisson(lambda, re);
	double backsize = xsum / (poiss.cdf(right) - poiss.cdf(left - 1));
	double backprop = Math.min(backsize / sampleTotals[samp.getIndex()], 1.0);
	System.out.println("Background= "+ backsize+" / "+sampleTotals[samp.getIndex()]+" =\t"+backprop);
	
	return backprop;
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
public double K(double L, int left, int right, boolean out){
	Poisson poiss = new Poisson(L, re); 
    if(out)
        System.out.println("K: "+poiss.cdf(right)+" "+poiss.cdf(left - 1)+" "+ L+" "+ left+" "+ right);
    return poiss.cdf(right) - poiss.cdf(left - 1);
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 * Reset duplicate reads that pass Poisson threshold. 
 * The Poisson lambda parameter is calculated by an Gaussian average
 * that puts more weight for nearby bases (same chrom, same strand)
 */
public void capPerBaseCountWithPoissonGaussianFilter(double threshold, int width){
       double g[] = new double[width*4+1];
	NormalDistribution gaussianDist = new NormalDistribution(0, width*width);
	for (int i=0;i<g.length;i++)
		g[i]=gaussianDist.calcProbability((double)i);
	
	DRand re = new DRand();
	Poisson P = new Poisson(0, re);
		
	for(int i = 0; i < fivePrimeCounts.length; i++)
		for(int j = 0; j < fivePrimeCounts[i].length; j++){
			float counts[] = fivePrimeCounts[i][j];
			int pos[] = fivePrimePos[i][j]; 
			if(counts!=null){
				for(int k = 0; k < counts.length; k++){
					int posK = pos[k]; 
					double sum = 0;
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k+x>=counts.length|| pos[k+x]-posK>width*4)
							break;
						sum += counts[k+x]*g[pos[k+x]-posK];
					}
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k-x<0 || posK-pos[k-x]>width*4)
							break;
						sum += counts[k-x]*g[posK-pos[k-x]];
					}
					sum = sum/(1-g[0]);				// exclude this position for evaluation
					
					double countThres=0;
					P.setMean(sum);
					double pvalue=1;
					for(int b=1; pvalue>threshold; b++){
						pvalue=1-P.cdf(b);	//p-value as the tail of Poisson
						countThres=b;
					}
					if (counts[k] > Math.max(1,countThres))
						counts[k] = (float) Math.max(1,countThres);					
				}
			}
		}
	updateTotalHits();
}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
	 * Reset duplicate reads that pass Poisson threshold. 
	 * The Poisson lambda parameter is calculated by an Gaussian average
	 * that puts more weight for nearby bases (same chrom, same strand)
	 */
	public void applyPoissonGaussianFilter(double threshold, int width){
        //init the Guassian kernel prob. for smoothing the read profile of called events
		double g[] = new double[width*4+1];
		NormalDistribution gaussianDist = new NormalDistribution(0, width*width);
		for (int i=0;i<g.length;i++)
			g[i]=gaussianDist.calcProbability((double)i);
		
		DRand re = new DRand();
		Poisson P = new Poisson(0, re);
			
		for(int i = 0; i < hitCounts.length; i++)
			for(int j = 0; j < hitCounts[i].length; j++){
				float counts[] = hitCounts[i][j];
				int pos[] = fivePrimes[i][j];
				for(int k = 0; k < counts.length; k++){
					int posK = pos[k]; 
					double sum = 0;
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k+x>=counts.length|| pos[k+x]-posK>width*4)
							break;
						sum += counts[k+x]*g[pos[k+x]-posK];
					}
					for (int x=1;x<=width*4;x++){		// at most extend out 250 idx
						if (k-x<0 || posK-pos[k-x]>width*4)
							break;
						sum += counts[k-x]*g[posK-pos[k-x]];
					}
					sum = sum/(1-g[0]);				// exclude this position for evaluation
					
					double countThres=0;
					P.setMean(sum);
					double pvalue=1;
					for(int b=1; pvalue>threshold; b++){
						pvalue=1-P.cdf(b);	//p-value as the tail of Poisson
						countThres=b;
					}
										
//					System.out.println(String.format("%s%d\t%.1f\t%.1f\t%.3f", (j==0?"+":"-"), posK, counts[k], countThres, sum));
					if (counts[k] > Math.max(1,countThres))
						counts[k] = (float) Math.max(1,countThres);					
				}
			}
		
		updateTotalHits();
	}
 

import cern.jet.random.Poisson; //导入依赖的package包/类
public DNASeqEnrichmentCaller() throws IOException, ClientException, SQLException {
      loader = new SeqDataLoader();
      DRand re = new DRand();
poisson = new Poisson(0, re);
      reads = new HMMReads();
  }
 

import cern.jet.random.Poisson; //导入依赖的package包/类
/**
 *
 * @return Poisson RNG
 */
public static Poisson getPoissonGenerator() {
    return poisson;
}