# TODO: Add comment
# 
# Author: E.Korsching Dec 7, 2012
###############################################################################


roc.ek <- function(x, pos.level=NULL, line=T, title="", c.groups=0, add=F){
	# ROC algorithm for generating (single or multiple) ROC points/curves
	# Lit.: ROC Graphs: Notes and Practical Considerations for Researchers, Tom Fawcett, 2004
	# Require: P > 0 and N > 0
	# Inputs: x: list: containing matrix/data.frame(s) with two columns:
	#  col 1: L(i), the set of test examples, True Class: e.g. outcome (dead or alive)
	#  col 2: f(L(i)), the probabilistic classifier’s estimate that example i is positive, Hypothesized Class: e.g. immunohistochemistry test score
	#  pos.level: number or character denoting a positive level
	#  line: T: line, F: points
	#  title: plot title or vector of titles, cgroups: number of color groups, max. 3
	#  lp.col: line point color, add: add graph to existing plot? (more than one roc curve in a figure)
	# P (TP+FN) and N (FP+TN), the number of positive and negative examples in L
	# AUC: average value, neither optimistic nor pessimistic
	# Outputs: RocP, a matrix of ROC points increasing by fp rate and AUC value
	
	#ini
	if(is.null(pos.level)){ cat("\nSupply positive level (number or character) - stop"); return() }
	l.len <- length(x)
	erg <- list(l.len)		#create list for results
	AUC.v <- vector(mode="numeric",length=l.len)		#create vector for all AUC values
	#create color groups
	if(c.groups<0 | c.groups>2){ cat("\n Number of color groups [0..2]"); return() }
	if(c.groups==0){ lp.col <- rep("black", l.len) }
	if(c.groups==1){ lp.col <- colorRampPalette(c("green", "orange"))( l.len ) }
	if(c.groups==2){ lp.col <- colorRampPalette(c("green", "grey", "orange"))( l.len ) }
	
	
	# process each roc data set
	for(i in 1:l.len){
		y <- x[[i]]
		Llevel <- unique(y[,1])
		Llevel.len <- length(Llevel)
		if(Llevel.len>2){ cat("\nOnly two levels allowed - stop in data set : ",i); return() }
		
		Ysorted <- y[order(y[,2],y[,1], decreasing=T),]		# L sorted decreasing by f scores
		Ysorted.len <- nrow(y)
		
		FP <- 0; TP <- 0
		P <- sum(Ysorted[,1]==pos.level)
		N <- sum(Ysorted[,1]!=pos.level)
		if((P+N)!=Ysorted.len){ cat("\nSum of P+N != number of rows - check data set : ",i); return() }
		RocP <- matrix(0,Ysorted.len+1,2)		# +1, last point
		
		fprev <- -Inf		# start value lower left corner
		j <- 1
		
		#
		while(j<=Ysorted.len){
			if(Ysorted[j,2]!=fprev){
				RocP[j,1] <- FP/N
				RocP[j,2] <- TP/P
				fprev <- Ysorted[j,2]
			}
			if(Ysorted[j,1]==pos.level){		# j is a positive example
				TP <- TP + 1
			}else{						# j is a negative example
				FP <- FP + 1
			}
			j <- j+1
		}
		RocP[j,1] <- FP/N		# Last operation: this is (1,1) upper right corner
		RocP[j,2] <- TP/P
		
		#AUC and plot preparation
		RocP <- RocP[(RocP[,1]!=0 | RocP[,2]!=0),]		#delete non informative rows (0,0)
		RocP <- rbind(c(0,0),RocP)			#restore first point (0,0)
		RocP.len <- nrow(RocP)
		AUC.c <- 0
		
		#plot roc
		if(i==1){
			plot(0,0, xlim=c(0,1), ylim=c(0,1), type="n", xlab="false positive rate (1-specificity)", ylab="true positive rate (sensitivity)", axes=T)
			segments(x0=0,y0=0,x1=1,y1=1, col="blue", lty=2, lwd=1.2)		#diagonal line
		}
		
		for(j in 1:(RocP.len-1)){			#len-1 !
			x0 <- RocP[j,1]
			y0 <- RocP[j,2]
			x1 <- RocP[j+1,1]
			y1 <- RocP[j+1,2]
			
			#AUC calculation
			AUC.c <- AUC.c + ((y0+y1)/2)*(x1-x0)		#trapezoid formula
			
			if(line==F){
				if(j==1){
					points(x=c(x0,x1), y=c(y0,y1), type="p", col=lp.col[i], pch=1, cex=0.6)		#two points
				}else{
					points(x=x1, y=y1, type="p", col=lp.col[i], pch=1, cex=0.6)				#one point
				}
			}else{
				segments(x0=x0,y0=y0,x1=x1,y1=y1, col=lp.col[i], lwd=1)				#line segments
			}
		}
		
		AUC.v[i] <- round(AUC.c,2)
		cat("\n AUC = ",AUC.v[i])		#show on cmd line
		erg[[i]] <- list(RocP=RocP,AUC=round(AUC.c,2))
	}
	
	#manage title and legend options
	if(add==F){
		if(title==""){ mtext(text=paste(date(),sep=""), side=3, line=1.5, cex=1) }else{ mtext(text=title, side=3, line=1.5, cex=1) }
	}
	
	if(!is.null(names(x))){		#produce legend text with data names and AUC values
		leg.str <- paste(AUC.v,names(x),sep=", ")
	}else{						#else go only with AUC values
		leg.str <- as.character(AUC.v)
	}
	
	if(line){
		legend(x=0.7, y=0.4, legend=leg.str,
				col=lp.col, border="black", lty=1, lwd=1,
				bty="o", title="AUC", cex=0.5)
	}else{
		legend(x=0.7, y=0.4, legend=leg.str,
				col=lp.col, border="black", pch=1,
				bty="o", title="AUC", cex=0.5)
	}
	text(x=0.7, y=0.05, labels=paste("Curves below blue line denote invers correlation of the classifier"), adj=0, cex=0.5, col=NULL, font=NULL)
	
	#
	return(erg)
}


#ab <- roc.ek(x=list(aa), pos.level="+", line=T, title="test roc curve (+/- vs. IH score)", c.groups=1, add=F)