library(componentSkeleton)
library(spp)

execute <- function(cf){
	# Inputs
	chip <- get.input(cf, 'treatment')
	mock <- get.input(cf, 'control' )

	# Params
	srange_max <- get.parameter(cf, 'srange_max', type = 'int')
	srange_min  <- get.parameter(cf, 'srange_min', type = 'int')
	bin  <- get.parameter(cf, 'bin', type = 'int')
	fdr  <- get.parameter(cf, 'fdr', type = 'float' )

	chip.data <- read.bam.tags(chip)
	input.data <- read.bam.tags(mock)

	cluster <- NULL

	# get binding info from cross-correlation profile
	# srange gives the possible range for the size of the protected region;
	# srange should be higher than tag length; making the upper boundary too high will increase calculation time
	#
	# bin - bin tags within the specified number of basepairs to speed up calculation;
	# increasing bin size decreases the accuracy of the determined parameters
	binding.characteristics <- get.binding.characteristics(chip.data,srange=c(50,500),bin=5,cluster=cluster);

	# print out binding peak separation distance
	#print(paste("binding peak separation distance =",binding.characteristics$peak$x))

	# plot cross-correlation profile
	#pdf(file="example.crosscorrelation.pdf",width=5,height=5)
	#par(mar = c(3.5,3.5,1.0,0.5), mgp = c(2,0.65,0), cex = 0.8);
	#plot(binding.characteristics$cross.correlation,type='l',xlab="strand shift",ylab="cross-correlation");
	#abline(v=binding.characteristics$peak$x,lty=2,col=2)
	#dev.off();

	# select informative tags based on the binding characteristics
	chip.data <- select.informative.tags(chip.data,binding.characteristics);
	input.data <- select.informative.tags(input.data,binding.characteristics);


	# binding detection parameters
	# desired FDR (1%). Alternatively, an E-value can be supplied to the method calls below instead of the fdr parameter
	# the binding.characteristics contains the optimized half-size for binding detection window
	detection.window.halfsize <- binding.characteristics$whs;
	  
	# determine binding positions using wtd method
	bp <- find.binding.positions(signal.data=chip.data,control.data=input.data,fdr=fdr,whs=detection.window.halfsize,cluster=cluster)

	print(paste("detected",sum(unlist(lapply(bp$npl,function(d) length(d$x)))),"peaks"));
	  
	# output detected binding positions
	output.binding.results(bp, get.output(cf, 'peaks'))
	#CSV.write(get.output(cf, 'peaks'), as.data.frame(bp))

	# Alternative:
	# bp <- find.binding.positions(signal.data=chip.data,control.data=input.data,fdr=fdr,method=tag.lwcc,whs=detection.window.halfsize,cluster=cluster)

}

main(execute)