<commit_before>package edu.psu.compbio.seqcode.projects.naomi.multiscalesignal;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;
import java.util.Set;
import java.util.SortedSet;
import java.util.TreeMap;
import java.util.TreeSet;

import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.distribution.NormalDistribution;

import edu.psu.compbio.seqcode.deepseq.experiments.ExptConfig;
import edu.psu.compbio.seqcode.genome.GenomeConfig;
import edu.psu.compbio.seqcode.genome.location.Region;
import edu.psu.compbio.seqcode.projects.naomi.utilities.MapUtility;
import edu.psu.compbio.seqcode.projects.seed.SEEDConfig;

/**
 * HyperstackSegmentation
 *
 * Methods for MultiScaleSignalRepresentation
 * Probabilistic Multiscale Image Segmentation, Vincken et al. IEEE (1997)
 * Multiscale representation of genomic signals, Knijnenburg et al. Nature Methods (2014)
 * 
 * @author <NAME>
 * 
 * copied from SegmentationTree to take more than one experiment target
 *
 **/

public class HyperstackSegmentation {
	
	protected GenomeConfig gconfig;
	protected ExptConfig econfig;
	protected SEEDConfig sconfig;
	protected int numScale;
	
	/*********************
	 * Gaussian scale space and window parameters	
	 */
	final static double DELTA_TAU = 0.5*Math.log(2);
	final static double MINIMUM_VALUE = Math.pow(10, -100); //arbitrary minimum value; I cannot use Double.MIN_VALUE because it becomes zero after several multiplication with Double.MIN_VALUE
//	final static double P_MIN = Math.pow(10,-3); // Knijnenburg paper uses this value
	final static double P_MIN = Math.pow(10,-2); //needs to be carefully determined because it will substantially affect Gaussian window size
	final static double P_BIN = 0.4995; // to determine window size of growing bin
	final static double K_MIN = 1/Math.sqrt(1-Math.exp(-2*DELTA_TAU));	
	final static double K_N = Math.ceil(K_MIN);

	/*********************
	 * Linkage weights
	 */
	final static double WEIGHT_I = 1.00;
	final static double WEIGHT_G = 0.0000001;
	final static double WEIGHT_M = 1000;
		
	public HyperstackSegmentation(GenomeConfig gcon, ExptConfig econ, SEEDConfig scon, int scale){	
		gconfig = gcon;
		econfig = econ;
		sconfig = scon;
		numScale = scale; 
	}	

	protected Map <Integer, Set<Integer>> buildTree (float[][][] gaussianBlur, Map <Integer, Integer> linkageMap, float[] DImax, int trailingZero, int zeroEnd){
		
		int currchromBinSize = gaussianBlur.length;
		int numConditions = gaussianBlur[0][0].length;
		float fchromBinSize = currchromBinSize;
		
		//cumGaussianBlur stores cumulative signal counts; this to be used to get the mean signal intensity
		float[][][] cumGaussianBlur = new float [gaussianBlur.length][2][numConditions];
		for (int i = 0; i<currchromBinSize; i++){
			for (int j = 0; j<2; j++){
				for (int k = 0 ; k < numConditions ; k++)
					gaussianBlur[i][j][k] = 0;
			}
		}
		for (int k = 0 ; k < numConditions ; k++){
			cumGaussianBlur[0][1][k] = gaussianBlur[0][1][k];
			for (int i = 1; i < currchromBinSize; i++){
				cumGaussianBlur[i][1][k] = cumGaussianBlur[i-1][1][k] + gaussianBlur[i][1][k];
			}
		}
		
		Map<Integer,Set<Integer>> segmentationTree =new HashMap<Integer,Set<Integer>>();
		segmentationTree.put(0, linkageMap.keySet());
		System.out.println("curr Scale 0 size, printing from segmentationTree "+segmentationTree.get(0).size());
		Set<Integer> startingNodes = new TreeSet<Integer>(segmentationTree.get(0));
		
//		Map<Region, HashMap<Integer,Set<Integer>>> segmentationTree = new HashMap<Region, HashMap<Integer, Set<Integer>>>();
		
		/*********************
		 * Matrices parameters
		 */	
		double[] sigma = new double[numScale];
		double[] radius = new double[numScale];
		for (int i = 0; i<numScale;i++){
			sigma[i] = 1;
			radius[i] = 1;
		}
		
		for (int n = 1; n<numScale; n++){
			
			final long gaussianStartTime = System.currentTimeMillis();
			
			//sigma calculation
			sigma[n] = Math.exp(n*DELTA_TAU);
			NormalDistribution normDistribution = new NormalDistribution(0.00,sigma[n]);
			//take inverse CDF based on the normal distribution using probability
			double inverseCDF = normDistribution.inverseCumulativeProbability(P_MIN);	
			double binInverseCDF = normDistribution.inverseCumulativeProbability(P_BIN);			
			int windowSize = (int) (-Math.round(inverseCDF)*2+1);	
			float binWindowSize = -Math.round(binInverseCDF)*2+1;
			
			//window calculation based on Gaussian(normal) distribution with sigma, mean=zero,x=X[i]			
			double window[] = new double[windowSize];
			double windowSum = 0;
			for (int i = 0;i<windowSize;i++){
				window[i] = normDistribution.density(Math.round(inverseCDF)+i);
				windowSum = windowSum+window[i];
			}
			double normalizedWindow[]=new double[windowSize];
			for (int i = 0;i<windowSize;i++)
				normalizedWindow[i] = window[i]/windowSum;				
			
			// smoothing signals for numConditions times
			for (int k = 0 ; k < numConditions; k++){
				
				double polyCoeffi[] = new double [(int) Math.ceil(fchromBinSize/binWindowSize)];
				System.out.println("binWindowSize is "+binWindowSize+"\tpolyCoeffi length is "+(int) Math.ceil(fchromBinSize/binWindowSize));
				//copy from column[1] to column[0];this procedure need to be repeated for each iteration of scale
				// copy from column[1] to array to store polynomial coefficient
				for (int i = 0 ; i<currchromBinSize; i++){
					gaussianBlur[i][0][k]=gaussianBlur[i][1][k];
					polyCoeffi[(int) Math.floor(((float) i)/binWindowSize)] += gaussianBlur[i][1][k]/binWindowSize;
				}
				
				for (int i = 0; i < Math.ceil(fchromBinSize/binWindowSize); i++){
					if (polyCoeffi[i] == 0)
						polyCoeffi[i] = MINIMUM_VALUE;
				}

				PolynomialFunction poly1 = new PolynomialFunction(polyCoeffi);
				PolynomialFunction poly2 = new PolynomialFunction(normalizedWindow);
				PolynomialFunction polyMultiplication=poly1.multiply(poly2);
				double coefficients[]= polyMultiplication.getCoefficients();
		
				//taking mid point of polynomial coefficients			
				int coeffiMid = (int) Math.floor(((float) coefficients.length)/ 2.0);
		
				System.out.println("currchromBin Size is : "+currchromBinSize+"\twindowSize is: "+windowSize+
						"\tpolyCoeffi length is "+polyCoeffi.length+"\tcoefficients length is: "+coefficients.length);
				
			//copy Gaussian blur results to the column[1] without increasing bin size
//			for (int i = 0; i<currchromBinSize;i++){
//				if (currchromBinSize % 2 ==0 && coefficients.length % 2 == 1){
//					gaussianBlur[i][1]=(float) coefficients[polyMid-currchromBinSize/2+i+1];
//				}else{
//					gaussianBlur[i][1]=(float) coefficients[polyMid-currchromBinSize/2+i];
//				}
			
				// copy Gaussian blur results to the column[1] with increasing bin size
				for (int i = 0; i<currchromBinSize;i++){
					if (polyCoeffi.length % 2 == 0 && coefficients.length % 2 == 1)
						gaussianBlur[i][1][k] = (float) coefficients[(int) (coeffiMid-Math.floor((fchromBinSize/2-i)/binWindowSize))+1];
					else
						gaussianBlur[i][1][k] = (float) coefficients[(int) (coeffiMid-Math.floor((fchromBinSize/2-i)/binWindowSize))];
				}			
			}
			
			// printing region is sensitive to bin size; current bin size is 100 bp
			for (int k = 0 ; k < numConditions ; k ++){
				for (int i = 0; i< 1151698;i += 20000)
					System.out.println(gaussianBlur[i][0][k]+" : "+gaussianBlur[i][1][k]);
			}
			
			final long gaussianEndTime = System.currentTimeMillis();
			
			System.out.println("scale "+n+" Gausisian blur execusion time "+ (gaussianEndTime-gaussianStartTime));
		
			/***************
			 * Search Volume
			 */ 	<commit_msg>Fixed a minor bug in HyperstackSegmentation<commit_after>package edu.psu.compbio.seqcode.projects.naomi.multiscalesignal;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.Map;
import java.util.Set;
import java.util.SortedSet;
import java.util.TreeMap;
import java.util.TreeSet;

import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math3.distribution.NormalDistribution;

import edu.psu.compbio.seqcode.deepseq.experiments.ExptConfig;
import edu.psu.compbio.seqcode.genome.GenomeConfig;
import edu.psu.compbio.seqcode.genome.location.Region;
import edu.psu.compbio.seqcode.projects.naomi.utilities.MapUtility;
import edu.psu.compbio.seqcode.projects.seed.SEEDConfig;

/**
 * HyperstackSegmentation
 *
 * Methods for MultiScaleSignalRepresentation
 * Probabilistic Multiscale Image Segmentation, Vincken et al. IEEE (1997)
 * Multiscale representation of genomic signals, Knijnenburg et al. Nature Methods (2014)
 * 
 * @author <NAME>
 * 
 * copied from SegmentationTree to take more than one experiment target
 *
 **/

public class HyperstackSegmentation {
	
	protected GenomeConfig gconfig;
	protected ExptConfig econfig;
	protected SEEDConfig sconfig;
	protected int numScale;
	
	/*********************
	 * Gaussian scale space and window parameters	
	 */
	final static double DELTA_TAU = 0.5*Math.log(2);
	final static double MINIMUM_VALUE = Math.pow(10, -100); //arbitrary minimum value; I cannot use Double.MIN_VALUE because it becomes zero after several multiplication with Double.MIN_VALUE
//	final static double P_MIN = Math.pow(10,-3); // Knijnenburg paper uses this value
	final static double P_MIN = Math.pow(10,-2); //needs to be carefully determined because it will substantially affect Gaussian window size
	final static double P_BIN = 0.4995; // to determine window size of growing bin
	final static double K_MIN = 1/Math.sqrt(1-Math.exp(-2*DELTA_TAU));	
	final static double K_N = Math.ceil(K_MIN);

	/*********************
	 * Linkage weights
	 */
	final static double WEIGHT_I = 1.00;
	final static double WEIGHT_G = 0.0000001;
	final static double WEIGHT_M = 1000;
		
	public HyperstackSegmentation(GenomeConfig gcon, ExptConfig econ, SEEDConfig scon, int scale){	
		gconfig = gcon;
		econfig = econ;
		sconfig = scon;
		numScale = scale; 
	}	

	protected Map <Integer, Set<Integer>> buildTree (float[][][] gaussianBlur, Map <Integer, Integer> linkageMap, float[] DImax, int trailingZero, int zeroEnd){
		
		int currchromBinSize = gaussianBlur.length;
		int numConditions = gaussianBlur[0][0].length;
		float fchromBinSize = currchromBinSize;
		
		//cumGaussianBlur stores cumulative signal counts; this to be used to get the mean signal intensity
		float[][][] cumGaussianBlur = new float [gaussianBlur.length][2][numConditions];
		for (int i = 0; i<currchromBinSize; i++){
			for (int j = 0; j<2; j++){
				for (int k = 0 ; k < numConditions ; k++)
					gaussianBlur[i][j][k] = 0;
			}
		}
		for (int k = 0 ; k < numConditions ; k++){
			cumGaussianBlur[0][1][k] = gaussianBlur[0][1][k];
			for (int i = 1; i < currchromBinSize; i++){
				cumGaussianBlur[i][1][k] = cumGaussianBlur[i-1][1][k] + gaussianBlur[i][1][k];
			}
		}
		
		Map<Integer,Set<Integer>> segmentationTree =new HashMap<Integer,Set<Integer>>();
		segmentationTree.put(0, linkageMap.keySet());
		System.out.println("curr Scale 0 size, printing from segmentationTree "+segmentationTree.get(0).size());
		Set<Integer> startingNodes = new TreeSet<Integer>(segmentationTree.get(0));
		
//		Map<Region, HashMap<Integer,Set<Integer>>> segmentationTree = new HashMap<Region, HashMap<Integer, Set<Integer>>>();
		
		/*********************
		 * Matrices parameters
		 */	
		double[] sigma = new double[numScale];
		double[] radius = new double[numScale];
		for (int i = 0; i<numScale;i++){
			sigma[i] = 1;
			radius[i] = 1;
		}
		
		for (int n = 1; n<numScale; n++){
			
			final long gaussianStartTime = System.currentTimeMillis();
			
			//sigma calculation
			sigma[n] = Math.exp(n*DELTA_TAU);
			NormalDistribution normDistribution = new NormalDistribution(0.00,sigma[n]);
			//take inverse CDF based on the normal distribution using probability
			double inverseCDF = normDistribution.inverseCumulativeProbability(P_MIN);	
			double binInverseCDF = normDistribution.inverseCumulativeProbability(P_BIN);			
			int windowSize = (int) (-Math.round(inverseCDF)*2+1);	
			float binWindowSize = -Math.round(binInverseCDF)*2+1;
			
			//window calculation based on Gaussian(normal) distribution with sigma, mean=zero,x=X[i]			
			double window[] = new double[windowSize];
			double windowSum = 0;
			for (int i = 0;i<windowSize;i++){
				window[i] = normDistribution.density(Math.round(inverseCDF)+i);
				windowSum = windowSum+window[i];
			}
			double normalizedWindow[]=new double[windowSize];
			for (int i = 0;i<windowSize;i++)
				normalizedWindow[i] = window[i]/windowSum;				
			
			// smoothing signals for numConditions times
			for (int k = 0 ; k < numConditions; k++){
				
				double polyCoeffi[] = new double [(int) Math.ceil(fchromBinSize/binWindowSize)];
				System.out.println("binWindowSize is "+binWindowSize+"\tpolyCoeffi length is "+(int) Math.ceil(fchromBinSize/binWindowSize));
				//copy from column[1] to column[0];this procedure need to be repeated for each iteration of scale
				// copy from column[1] to array to store polynomial coefficient
				for (int i = 0 ; i<currchromBinSize; i++){
					gaussianBlur[i][0][k]=gaussianBlur[i][1][k];
					polyCoeffi[(int) Math.floor(((float) i)/binWindowSize)] += gaussianBlur[i][1][k]/binWindowSize;
				}
				
				for (int i = 0; i < Math.ceil(fchromBinSize/binWindowSize); i++){
					if (polyCoeffi[i] == 0)
						polyCoeffi[i] = MINIMUM_VALUE;
				}

				PolynomialFunction poly1 = new PolynomialFunction(polyCoeffi);
				PolynomialFunction poly2 = new PolynomialFunction(normalizedWindow);
				PolynomialFunction polyMultiplication=poly1.multiply(poly2);
				double coefficients[]= polyMultiplication.getCoefficients();
		
				//taking mid point of polynomial coefficients			
				int coeffiMid = (int) Math.floor(((float) coefficients.length)/ 2.0);
		
				System.out.println("currchromBin Size is : "+currchromBinSize+"\twindowSize is: "+windowSize+
						"\tpolyCoeffi length is "+polyCoeffi.length+"\tcoefficients length is: "+coefficients.length);
				
			//copy Gaussian blur results to the column[1] without increasing bin size
//			for (int i = 0; i<currchromBinSize;i++){
//				if (currchromBinSize % 2 ==0 && coefficients.length % 2 == 1){
//					gaussianBlur[i][1]=(float) coefficients[polyMid-currchromBinSize/2+i+1];
//				}else{
//					gaussianBlur[i][1]=(float) coefficients[polyMid-currchromBinSize/2+i];
//				}
			
				// copy Gaussian blur results to the column[1] with increasing bin size
				for (int i = 0; i<currchromBinSize;i++){
					if (polyCoeffi.length % 2 == 0 && coefficients.length % 2 == 1)
						gaussianBlur[i][1][k] = (float) coefficients[(int) (coeffiMid-Math.floor((fchromBinSize/2-i)/binWindowSize))+1];
					else
						gaussianBlur[i][1][k] = (float) coefficients[(int) (coeffiMid-Math.floor((fchromBinSize/2-i)/binWindowSize))];
				}			
			}
			
			// printing region is sensitive to bin size; current bin size is 100 bp
			for (int k = 0 ; k < numConditions ; k ++){
				System.out.println("condition "+k);
				for (int i = 0; i< 1151698;i += 20000)
					System.out.println(gaussianBlur[i][0][k]+" : "+gaussianBlur[i][1][k]);
			}
			
			final long gaussianEndTime = System.currentTimeMillis();
			
			System.out.println("scale "+n+" Gausisian blur execusion time "+ (gaussianEndTime-gaussianStartTime));
		
			/***************
			 * Search Volume
			 */ 	
{'fraction_non_alphanumeric': 0.13824524863316845, 'fraction_numerical': 0.01822442072376985, 'mean_word_length': 4.251196172248804, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3877195', 'n_tokens_mistral': 5897, 'n_tokens_neox': 4955, 'n_words': 1237}