# TODO: Add comment
# 
# Author: E.Korsching Sep 11, 2013
###############################################################################


coverage.vec <- function(x, window.x){
	# calculates the coverage of a physical position range
	# extracts from a simple  scanBam() object  one list  with  pos and qwidth  and  optional score  object
	# note: depending on the read length there might be boundary effects
	# window.x: c(min,max), the view is downstream
	
	# check x
	if(length(x)==2){ flag <- 0 }
	if(length(x)==3){ flag <- 1 }
	# normalize window
	window.width <- window.x[2]-window.x[1]+1		#+1 because including first
	# matrix col1: position, col2: coverage
	erg <- matrix(0,window.width,2)
	erg[,1] <- seq(window.x[1],window.x[2],1)
	# extract pos and qwidth objects
	pos <- x[["pos"]]
	qwidth <- x[["qwidth"]]
	if(flag==1){ score <- x[["score"]] }
	# to be sure order
	pos.o <- order(pos)
	pos <- pos[pos.o]
	qwidth <- qwidth[pos.o]
	if(flag==1){ score <- score[pos.o] }
	# how many positions 
	lpos <- length(pos)		# greater 28000 might be wrong (read length 35)
	lqwidth <- length(qwidth)
	if(lpos!=lqwidth){ cat("\n length of pos and qwidth differ"); return(matrix(0,1,2)) }	###
	# control qwidth
	qwidth[qwidth<1] <- 1		# RLE
	# normalize pos
	pos <- pos - (window.x[1] -1)		#-1 because starting at 1 not at 0
	pos.qwidth <- pos + (qwidth-1)		# RLE
	# control deviations / clip
	pos[pos<1] <- 1
	if(sum(pos>window.width)==lpos){ cat("\n check for error: matrix: [0,0] ; read number: ",lpos); return(matrix(0,1,2)) }	###
	pos.qwidth[pos.qwidth>window.width] <- window.width
	# coverage   and clip if necessary
	for(i in 1:lpos){
		tmpA <- pos[i]
		tmpB <- pos.qwidth[i]
#		assign("aa", list(pos,pos.qwidth,lpos), envir = .GlobalEnv)
#		cat("\n i ",i," ",tmpA," : ",tmpB)
		if(flag==0){
			erg[tmpA:tmpB, 2] <- erg[tmpA:tmpB, 2] +1		#add count: +1
		}else{
			erg[tmpA:tmpB, 2] <- erg[tmpA:tmpB, 2] +score[i]		#add score value
		}
	}
	return(erg)
}



#ac <- coverage.vec( list( pos=as.integer(c(18,25,25,25,38)), qwidth=as.integer(c(10,10,10,10,10)) ) , c(20,50) )
#ac <- coverage.vec( list( pos=as.integer(c(21,25,25,25,46)), qwidth=as.integer(c(10,10,10,10,10)) ) , c(20,50) )
#ac <- coverage.vec(a.tmp[[1]], c(3556531,3557531))