
###################################################
### recommended RStudio to run the following code
### NOTICE: sometimes you need to set working directory, e.g.:
### setwd("/home/yang/share/empire")
###################################################

source("plot.R")

# library(car)
# library(perturb)
# library(plyr, quietly=TRUE)
# library(usdm, quietly=TRUE)
# library(effsize, quietly=TRUE)

RemoveConfoundingEffect <- function(data, ynames, confeff.name="LOC")
{
	#### confeff.name: name of the variable in data that each of ynames have confounding effect with
	for (yname in ynames) {
		if (yname==confeff.name) { next }

		formula <- as.formula(paste(yname, confeff.name, sep=" ~ "))
		model <- glm(formula=formula, data=data, family=gaussian)
		data[, yname] <- data[, yname]-model$fitted.values
	}
	data
}

### stepwise logistic regression variable selection
Stepwise <- function(data, yname, xnames, confounding.effect.name="LOC", rule="aic", direction="forward", remove.flag=FALSE, trace=FALSE)
{
	### stepwise variable selection
	Stepwise.One <- function(data, yname, xnames, rule="aic", direction="forward", trace=FALSE)
	{
		cat("step", rule, direction, "\n")

		rule <- tolower(rule)
		direction <- tolower(direction)

		null.formula <- as.formula(sprintf("%s ~ 1", yname))
		full.formula <- as.formula(sprintf("%s ~ %s", yname, paste(xnames, collapse=" + ")))

		model <- null.model <- full.model <- NULL
			null.model <- glm(null.formula, family=binomial, data=data)
			full.model <- glm(full.formula, family=binomial, data=data)
		
		if (rule=="aic") { k <- 2 } else { k <- log(nrow(data)) }
		scope <- list(upper=full.model, lower=null.model)
		if (direction=="forward") {
			model <- step(null.model, scope=scope, direction="forward",  k=k, trace=trace)
		} else if (direction=="backward") {
			model <- step(full.model, scope=scope, direction="backward", k=k, trace=trace)
		} else if (direction=="forwardboth") {
			model <- step(null.model, scope=scope, direction="both",     k=k, trace=trace)
		} else if (direction=="backwardboth") {
			model <- step(full.model, scope=scope, direction="both",     k=k, trace=trace)
		}

		return(model)
	}

	if (length(xnames)==0) { stop("Length of independent variable is zero. Please check.\n") }

	data.new <- data
	if (remove.flag==TRUE) {
		if (trace) { cat("Removing confounding effect for independent variables before stepwise variable selection.\n") }
		data.new <- RemoveConfoundingEffect(data.new, ynames=xnames, confeff.name=confounding.effect.name)
	}

	model <- Stepwise.One(data=data.new, yname=yname, xnames=xnames, rule=rule, direction=direction, trace=trace)

	return(model)
}

### this function sort the data frame base on the predicted value and second by other policies
SortData <- function(data, effortaware=FALSE, sorted=FALSE, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE)
{
	if (!all(c("NUM", "REL", "LOC", "PRE") %in% colnames(data))) { stop("ERROR: NUM, REL, LOC, or PRE is not colnames of data frame") }

	key.2nd <- "REL"
	if (effortaware) {
        key.2nd <- "density"
		if (!(key.2nd %in% colnames(data))) {
			data[[key.2nd]] <- data$NUM/(data$LOC+1)
		}
		sorted <- FALSE
	}

	if (!sorted) {
		if (worstcase) {
			data <- data[order(-data$PRE, +data[[key.2nd]], -data$LOC), ]
		} else if (bestcase) {
			data <- data[order(-data$PRE, -data[[key.2nd]], +data$LOC), ]
		} else if (LOCUP) {
			data <- data[order(-data$PRE, +data$LOC), ]
		} else {
			data <- data[order(-data$PRE), ]
		}
		sorted <- TRUE
	}

	return(data)
}

#####################################################################
### ER : effort reduction compared to random model (model evaluation criteria for cost-effectiveness)
#####################################################################
ComputeERs <- function(data, cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), effortaware=TRUE, sorted=FALSE, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE, allpercentcutoff=TRUE)
{
	ComputeER <- function(data, cutoff=0.2, effortaware=TRUE, sorted=FALSE, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE)
	{
		if (!sorted) {
			data <- SortData(data=data, effortaware=effortaware, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)
			sorted <- TRUE
		}

		data$PREBIN <- 0
		len <- nrow(data) * cutoff
		data$PREBIN[1:len] <- 1

		totalLOC   <- sum(data$LOC)
		totalLOC.M <- sum(data$PREBIN*data$LOC)

		totalNUM   <- sum(data$NUM)
		totalNUM.R <- sum(data$PREBIN*data$REL*data$NUM)

		if (totalNUM==0) { stop("of this testing data set, total number of NUM is zero") }
		if (totalLOC==0) { stop("of this testing data set, total number of LOC is zero") }

		Effort.M <- totalLOC.M / totalLOC ### effort of model m
		Effort.R <- totalNUM.R / totalNUM ### effort of random model

		ER <- Effort.R - Effort.M ### 绝对的ER 相对于整个系统代码审查的百分比
		return(ER)
	}

	ERs <- NULL
	for (cutoff in cutoffs) {
		aER <- ComputeER(data=data, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP, cutoff=cutoff)
		# names(aER) <- paste(names(aER), "P", cutoff, sep=".")
		# ERs <- c(ERs, unlist(aER))
		ERs <- c(ERs, aER)
	}

	names(ERs) <- paste("ER.P", cutoffs, sep=".")

	return(ERs)
}

#####################################################################
### CE : cost-effectiveness based on Alberg diagram (model evaluation criteria for cost-effectiveness)
#####################################################################
ComputeCEs <- function(data, cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), effort.percent.flag=TRUE, effortaware=TRUE, sorted=FALSE, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE)
{
	### compute area under curve for CE
	ComputeCEAreawithCutoffs <- function(data, cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), effort.percent.flag=TRUE, effortaware=TRUE, sorted=FALSE, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE)
	{
		if (!sorted) {
			data <- SortData(data=data, effortaware=effortaware, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)
			sorted <- TRUE
		}

		len <- nrow(data)
		### cumulative sums of LOC or NUM
		cumXs <- cumsum(data$LOC) # x: LOC%
		cumYs <- cumsum(data$NUM) # y: Bug%

		Xs <- cumXs/cumXs[len]
		Ys <- cumYs/cumYs[len]

		poses <- vector(length=length(cutoffs))
		cutoffs.new <- NULL
		if (effort.percent.flag) {
			for (i in seq(length.out=length(cutoffs))) {
				### find the lowest index which is larger or equal than cutoff
				### that's to say cutoff is between pos-1 and pos
				cutoff   <- cutoffs[i]
				poses[i] <- min(which(Xs >= cutoff))
			}
			cutoffs.new <- cutoffs ### x-axis is the corresponding cutoff (effort)
		} else {
			poses <- len*cutoffs ### which
			cutoffs.new <- Xs[poses] ### x-axis
		}

		Yposes <- vector(length=length(cutoffs))
		if (effort.percent.flag) {
			for (i in seq(length.out=length(cutoffs))) {
				cutoff <- cutoffs[i]
				pos <- poses[i]

				if (pos==1) {
					Yposes[i] <- cutoff * (Ys[1]/Xs[1])
				} else {
					Yposes[i] <- (cutoff-Xs[pos-1]) * ((Ys[pos]-Ys[pos-1])/(Xs[pos]-Xs[pos-1])) + Ys[pos-1]
				}
			}
		} else {
			Yposes <- Ys[poses]
		}

		fix_subareas <- vector(length=len)
		fix_subareas[1] <- 0.5 * Ys[1] * Xs[1]
		fix_subareas[2:len] <- 0.5 * (Ys[1:(len-1)] + Ys[2:len]) * abs(Xs[1:(len-1)] - Xs[2:len])

		areas      <- vector(length=length(cutoffs))
		RECALLs    <- vector(length=length(cutoffs))
		for (i in seq(length.out=length(poses))) {
			pos <- poses[i]
			Ypos <- Yposes[i]
			cutoff <- cutoffs.new[i]

			subareas <- vector(length=pos)
			if (pos==1) {
				subareas[1] <- 0.5 * cutoff * Ypos
			} else if (pos==2) {
				subareas[1] <- fix_subareas[1]
				subareas[2] <- (Ypos+Ys[1])*(abs(cutoff-Xs[1]))*0.5
			} else {
				subareas[1:(pos-1)] <- fix_subareas[1:(pos-1)]
				subareas[pos]  <- (Ypos+Ys[pos-1])*(abs(cutoff-Xs[pos-1]))*0.5
			}

			areas[i]      <- sum(subareas)
		}

		return(list(AREA=areas))
	}

	if (!sorted) {
		data <- SortData(data=data, effortaware=effortaware, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)
		sorted <- TRUE
	}

	data.mdl <- data
	data.opt <- data[order(-data$density, +data$LOC), ]
	data.min <- data[order(+data$density, -data$LOC), ]

	areas.opt <- ComputeCEAreawithCutoffs(data=data.opt, cutoffs=cutoffs, effort.percent.flag=effort.percent.flag, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)
	areas.mdl <- ComputeCEAreawithCutoffs(data=data.mdl, cutoffs=cutoffs, effort.percent.flag=effort.percent.flag, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)

	areas.random <- 0.5 * cutoffs * cutoffs

	CEs <- (areas.mdl$AREA - areas.random)/(areas.opt$AREA - areas.random)

	CEs <- pmin(1, CEs)
	CEs[is.na(CEs)] <- 0

	# if (any(is.na(ARs)) | any(is.na(CEs))) { cat("one element of ARs or CEs is NA") }

	if (effort.percent.flag) {
		names(CEs) <- paste("CE.E", cutoffs, sep=".")
	} else {
		names(CEs) <- paste("CE.P", cutoffs, sep=".")
	}

	return(CEs)
}

Compute <- function(data.train, data.valid, yname, xnames=NULL, xnames1=NULL, xnames2=NULL, run=0, fold=0, model="lr", formula=NULL, family="binomial", effort.name="LOC", number.name="BUG", cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), trace=FALSE)
{
	if (trace) { cat("cross validation at (run: ", run, ", fold: ", fold, ")", "\n") }

	data.new.valid <- data.valid
	data.new.train <- data.train

	### Note: do not need to select variable agian(It was determined by previous selection
	### based on previous model formula and train dataset to obtain model, and valid it by valid dataset
	valid.dt.divide <- train.dt.divide <- NULL
	valid.pre <- train.pre <- NULL

	if (is.null(formula)) {
        if (!is.null(xnames)) { formula <- as.formula(sprintf("%s ~ %s", yname, paste(xnames, collapse="+"))) }
        
		if (!is.null(xnames1)) { if (length(xnames1)==0) { formula1 <- as.formula(sprintf("%s ~ 1", yname)) } else { formula1 <- as.formula(sprintf("%s ~ %s", yname, paste(xnames1, collapse="+"))) } }
        if (!is.null(xnames2)) { if (length(xnames2)==0) { formula2 <- as.formula(sprintf("%s ~ 1", yname)) } else { formula2 <- as.formula(sprintf("%s ~ %s", yname, paste(xnames2, collapse="+"))) } }
	}

	model <- tolower(model)
	if (model=="lr" || model=="lr.2phase") {
		if (is.null(xnames) && model=="lr.2phase") {
			mdl <- NULL

			data.new.train.1 <- data.new.train[data.new.train$C.in.ALLD==1, ]
			data.new.train.2 <- data.new.train[data.new.train$C.in.ALLD==0, ]
			data.new.train <- rbind(data.new.train.1, data.new.train.2)

			mdl.1 <- glm(formula=formula1, family=family, data=data.new.train.1)
			mdl.2 <- glm(formula=formula2, family=family, data=data.new.train.2)

			data.new.valid.1 <- data.new.valid[data.new.valid$C.in.ALLD==1, ]
			data.new.valid.2 <- data.new.valid[data.new.valid$C.in.ALLD==0, ]
			data.new.valid <- rbind(data.new.valid.1, data.new.valid.2)

			valid.pre.1 <- predict(mdl.1, data.new.valid.1, type="response")
			valid.pre.2 <- predict(mdl.2, data.new.valid.2, type="response")
			valid.pre <- c(valid.pre.1, valid.pre.2)

			train.pre.1 <- mdl.1$fitted.values
			train.pre.2 <- mdl.2$fitted.values
			train.pre <- c(train.pre.1, train.pre.2)
		} else {
			mdl <- glm(formula=formula, family=family, data=data.new.train)

			valid.pre <- predict(mdl, data.new.valid, type="response")
			train.pre <- mdl$fitted.values
		}
	}

	### divide (effort-aware)
	train.dt.divide.effort <- data.frame(REL=data.new.train[, yname], PRE=train.pre/(data.new.train[, effort.name] + 1), LOC=data.new.train[, effort.name], NUM=data.new.train[, number.name])
	valid.dt.divide.effort <- data.frame(REL=data.new.valid[, yname], PRE=valid.pre/(data.new.valid[, effort.name] + 1), LOC=data.new.valid[, effort.name], NUM=data.new.valid[, number.name])

	if (!is.numeric(train.dt.divide.effort$REL)) {
		train.dt.divide.effort$REL <- as.numeric(as.character(train.dt.divide.effort$REL))
	}

	if (!is.numeric(valid.dt.divide.effort$REL)) {
		valid.dt.divide.effort$REL <- as.numeric(as.character(valid.dt.divide.effort$REL))
	}

	if (any(is.na(valid.pre))) { stop("valid.pre is NA") }
	if (any(is.na(train.pre))) { stop("train.pre is NA") }

	subcompute <- function(train.dt, valid.dt, worstcase=FALSE, bestcase=FALSE, LOCUP=FALSE)
	{
		train.dt.ea <- SortData(data=train.dt, effortaware=TRUE, sorted=FALSE) ### effort-aware
		valid.dt.ea <- SortData(data=valid.dt, effortaware=TRUE, sorted=FALSE) ### effort-aware

		sorted <- TRUE
		allpercentcutoff <- TRUE

		effortaware <- TRUE
		CE <- ComputeCEs(data=valid.dt.ea, cutoffs=c(0.2), effortaware=effortaware, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)
		ER <- ComputeERs(data=valid.dt.ea, cutoffs=c(0.2), effortaware=effortaware, sorted=sorted, worstcase=worstcase, bestcase=bestcase, LOCUP=LOCUP)

		return(list(ER=ER, CE=CE))
	}

	divide.effort <- subcompute(train.dt=train.dt.divide.effort, valid.dt=valid.dt.divide.effort) ### column LOC with different values

	return(list(divide.effort=divide.effort))
}

### one fold of cross-validation
OneFold <- function(data, yname, xnames=NULL, xnames1=NULL, xnames2=NULL, run=0, fold=0, model="lr", formula=NULL, family="binomial", effort.name="LOC", number.name="BUG", cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), remove.at.fold.flag=FALSE, remove.flag=FALSE, trace=FALSE)
{
	data.valid <- data[data$fold==fold, ]
	data.train <- data[data$fold!=fold, ]

	if (remove.at.fold.flag) { ### default is FALSE
		if (remove.flag=="separate") {
			if (length(xnames)!=0) {
				data.valid <- RemoveConfoundingEffect(data.valid, ynames=xnames, confeff.name=effort.name)
				data.train <- RemoveConfoundingEffect(data.train, ynames=xnames, confeff.name=effort.name)
			}
		} else if (remove.flag=="together") {
			if (length(xnames)!=0) {
				data.new <- RemoveConfoundingEffect(data, ynames=xnames, confeff.name=effort.name)
				data.valid <- data.new[data.new$fold==fold, ]
				data.train <- data.new[data.new$fold!=fold, ]
			}
		}
	}

	Compute(data.train=data.train, data.valid=data.valid, fold=fold, model=model, formula=formula, family=family, yname=yname, xnames=xnames, xnames1=xnames1, xnames2=xnames2, effort.name=effort.name, number.name=number.name, run=run, cutoffs=cutoffs, trace=trace)
}

### one time of XXX-fold cross-validation
OneRun <- function(run, data, yname, xnames=NULL, xnames1=NULL, xnames2=NULL, totalRuns=10, totalFolds=10, model="lr", formula=NULL, family="binomial", effort.name="LOC", number.name="BUG", cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), remove.at.fold.flag=FALSE, remove.flag=FALSE, parallel=FALSE, trace=FALSE)
{
	### split each class in data into number folds
	SplitData <- function(data, yname, folds, seedval)
	{
		classes <- unique(data[, yname]) ### multiple classes

		set.seed(seedval)
		seedval <- sample(1:1000000, 1) ### global seed

		data$fold <- 0
		for (aclass in classes) { ### 对数据集中每一类都分类fols等份
			index <- which(data[, yname]==aclass)

			set.seed(seedval)
			seq <- sample(1:length(index))

			data$fold[index] <- rep(1:folds, length.out=length(index))[seq]
		}
		return(data)
	}

	seed <- totalRuns * 100 + totalFolds * 10 + run
	set.seed(seed)
	seedval  <- sample(1:1000000, 1)

	data.new <- SplitData(data=data, yname=yname, folds=totalFolds, seedval=seedval)

	if (trace) { cat("cross validation at run: ", run, "(total runs:", totalRuns, ") for each of fold(total folds:", totalFolds, ").\n") }

	resList.FOLDS <- NULL
	if (parallel && (.Platform$OS.type=="unix")) {
		require(parallel) ### this is for multi-threading
		numWorkers <- 4
		resList.FOLDS <- mclapply(c(1:totalFolds), FUN=OneFold, data=data.new, model=model, formula=formula, family=family, yname=yname, xnames=xnames, xnames1=xnames1, xnames2=xnames2, effort.name=effort.name, number.name=number.name, run=run, remove.at.fold.flag=remove.at.fold.flag, remove.flag=remove.flag, cutoffs=cutoffs, , trace=trace, mc.cores=numWorkers)
	} else {
		resList.FOLDS <-   lapply(c(1:totalFolds), FUN=OneFold, data=data.new, model=model, formula=formula, family=family, yname=yname, xnames=xnames, xnames1=xnames1, xnames2=xnames2, effort.name=effort.name, number.name=number.name, run=run, remove.at.fold.flag=remove.at.fold.flag, remove.flag=remove.flag, cutoffs=cutoffs, trace=trace)
	}

	names(resList.FOLDS) <- paste("FOLD_", 1:totalFolds, sep="")
	return(resList.FOLDS)
}

###########################################################################
### cross-validation for machine learning models or logistic regression model
### remove flag can be "separate", "together", FALSE which is implemented in the 'Compute' function
######## parameters description ########
### data: data frame {column names with yname(dependent variable name) and xnames(indpendent variable names)}
### yname: dependent variable name, e.g. yname="buggy" indicate whether a module is buggy or not
### xnames: independent variables names, e.g. xnames=c("x1", "x2", "x3")
### effort.name: the proxy for computing cost-effectiveness, e.g. effort.name="LOC" indicate the number lines of code (effort.name should be a column in the 'data' data.frame)
### number.name: the number of bugs, e.g. number.name="bugs" indicate the number of bugs in a class (number.name should in the 'data' data.frame)
### totalFolds: number of folds, e.g. if m*n folds, the totalFolds indicate n
### totalRuns: number of runs, e.g. if m*n folds, the totalRuns is m
### others could be set as default
###########################################################################
CrossValidation <- function(data, yname, xnames=NULL, xnames1=NULL, xnames2=NULL, totalFolds=10, totalRuns=10, model="lr", formula=NULL, family="binomial", effort.name="LOC", number.name="BUG", cutoffs=c(0.01, 0.05, 0.1, 0.15, seq(from=0.2, to=1.0, by=0.1)), remove.at.fold.flag=FALSE, remove.flag=FALSE, parallel=FALSE, trace=FALSE)
{
	### runs n times m-fold cross-validation
	if (remove.at.fold.flag==FALSE) {
		### if not remove at fold thus remove confounding effect for the entire data set
		if ((remove.flag!=FALSE)&&(length(xnames)!=0)) {
			data <- RemoveConfoundingEffect(data, ynames=xnames, confeff.name=effort.name)
		}
	}

	if (trace) { cat("Cross Validation for each run.\n") }

	resList.RUNS <- lapply(c(1:totalRuns), FUN=OneRun, data=data, model=model, formula=formula, family=family, yname=yname, xnames=xnames, xnames1=xnames1, xnames2=xnames2, effort.name=effort.name, number.name=number.name, totalRuns=totalRuns, totalFolds=totalFolds, remove.at.fold.flag=remove.at.fold.flag, remove.flag=remove.flag, cutoffs=cutoffs, parallel=parallel, trace=trace)

	names(resList.RUNS) <- paste("RUN_", 1:totalRuns, sep="")

	subunlist <- function(theList, type.name="normal") {
		all.theList <- sapply(sapply(theList, '['), '[', type.name)

		CEs <- ERs <- NULL
		for (i in seq(all.theList)) {
			CEs <- rbind(CEs, all.theList[[i]][["CE"]])
			ERs <- rbind(ERs, all.theList[[i]][["ER"]])
		}

		return(list(ER=ERs, CE=CEs))
	}

	divide.effort <- subunlist(theList=resList.RUNS, type.name="divide.effort")

	return(list(divide.effort=divide.effort))
}

main <- function(filename, xnames, feature.names, class.name="BUG", number.name="BUG", model="lr.2phase", rule="aic", direction="forward", remove.flag=FALSE, cv.remove.flag=NULL, totalFolds = 10, totalRuns = 10, cutoffs=seq(from=0.1, to=1.0, by=0.1), trace=FALSE, parallel=FALSE)
{
	cat("RunMain for", filename, "\n")
	data <- read.csv(file=filename)
	fbasename <- basename(filename)

	for(feature.group.name in names(feature.names)) {
		cat("Analysis ", feature.group.name, "for ", fbasename, "\n")

		cat("RunMain...\n")
		if (is.null(cv.remove.flag)) { cv.remove.flag <- remove.flag }

		prefix <- sprintf("file(%s)_model(%s)_rule(%s)_direction(%s)_remove(%s)_%s", basename(filename), model, rule, direction, remove.flag, feature.group.name)
		fout <- sprintf("output/%s.out", prefix)

		### class name as y-axis name ###
		yname <- class.name

		### filter all valid xname
		all.xnames <- feature.names[[feature.group.name]]

		temp.xnames <- NULL
		for (xname in all.xnames) {
			if (length(unique(as.vector(data[, xname]))) > 1) {
				temp.xnames <- c(temp.xnames, xname)
			} else {
				cat(xname, " is removed since all instances with the same value\n")
			}
		}

		all.xnames <- temp.xnames
		if (length(all.xnames)==0) { return(FALSE) }
		if (!file.exists("output")) { dir.create("output") }

		subsets <- subsets.i <- subsets.o <- NULL
		
		outPath <- file.path("output", "model")
		if (!file.exists(outPath)) { dir.create(file.path(outPath)) }
		if (model=="lr.2phase") {
			if (feature.group.name=="B+C") {
				cat("Cross validation for ", feature.group.name, model, "\n")

				fmodel.o <- sprintf("output/model/%s-o.mdl", prefix)
				if (file.exists(fmodel.o)) {
					subsets.o <- dget(file=fmodel.o)
				} else {
					data.o <- data[data$C.in.ALLD==0, ]
					cat("stepwise for instances in the second class\n")
					mdl.o <- Stepwise(data=data.o, yname=yname, xnames=all.xnames, rule=rule, direction=direction, remove.flag=remove.flag, trace=trace)
					subsets.o <- attr(mdl.o$terms, "term.labels")
					dput(subsets.o, file=fmodel.o)
				}
				
				fmodel.i <- sprintf("output/model/%s-i.mdl", prefix)
				if (file.exists(fmodel.i)) {
					subsets.i <- dget(file=fmodel.i)
				} else {
					data.i <- data[data$C.in.ALLD==1, ]
					cat("stepwise for instances in the first class\n")
					mdl.i <- Stepwise(data=data.i, yname=yname, xnames=all.xnames, rule=rule, direction=direction, remove.flag=remove.flag, trace=trace)
					subsets.i <- attr(mdl.i$terms, "term.labels")
					dput(subsets.i, file=fmodel.i)
				}
			} else {
				fmodel <- sprintf("output/model/%s.mdl", prefix)
				if (file.exists(fmodel)) {
					subsets <- dget(file=fmodel)
				} else {
					mdl <- Stepwise(data=data, yname=yname, xnames=all.xnames, rule=rule, direction=direction, remove.flag=remove.flag, trace=trace)
					subsets <- attr(mdl$terms, "term.labels")

					dput(subsets, file=fmodel)
				}
			}
		}

		if (length(subsets)==0 && length(subsets.i)==0 && length(subsets.o)==0) { cat("subsets is empty\n"); return(FALSE) }
		outs <- CrossValidation(model=model, data=data, yname=yname, xnames=subsets, xnames1=subsets.i, xnames2=subsets.o, effort.name="LOC", number.name=number.name, totalFolds=totalFolds, totalRuns=totalRuns, cutoffs=cutoffs, remove.flag=cv.remove.flag, trace=trace, parallel=parallel)
		dput(outs, file=fout)
	}
}

Des <- function(filenames, trace=FALSE, parallel=FALSE)
{
	stats <- NULL
	for (filename in filenames) {
		cat("RQ1 for ", filename, "\n")

		data <- read.csv(file=filename)

		TF <- nrow(data) ### 函数总数
		TL <- sum(data$SLOC) ### 代码行总数

		tb <- sum(data$BUG) ### 缺陷函数总数
		tf <- TF - tb ### 无缺陷函数总数

		ti <- sum(data$C.in.ALLD)   ### 簇内函数总数
		to <- TF - ti ### 簇外函数总数

		tib <- sum(data$BUG*data$C.in.ALLD) ### 在簇内并且有缺陷
		tob <- sum(data$BUG*(1-data$C.in.ALLD)) ### 在簇外并且有缺陷

		tif <- sum((1-data$BUG)*data$C.in.ALLD) ### 在簇内并且没有缺陷
		tof <- sum((1-data$BUG)*(1-data$C.in.ALLD)) ### 在簇外并且没有缺陷

		ft <- fisher.test(x=rbind(c(tib, tif), c(tob, tof)), alternative="two.sided") ### fisher test 比较缺陷密度
		pv <- ft$p.value

		p <- tib/ti ### 簇中缺陷比例
		q <- tob/to ### 簇外缺陷比例
		OR <- (p/(1-p))/(q/(1-q))

		tbl <- table(data$membership.ALLD)
		cnum <- length(which(tbl>=2)) ### 簇个数
		lsize <- tbl[which.max(tbl)] ### 最大簇的size

		stat <- unlist(c("SLOC"=TL, "#functions"=TF, "#F is buggy"=tb, "%F is buggy"=tb/TF*100, "#clusters"=cnum, "Max cluster size"=lsize, "#F inside clusters"=ti, "%F inside clusters"=ti/TF*100, "#F outside clusters"=to, "%F outside clusters"=to/TF*100, "#F inside clusters is buggy"=tib, "#F outside clusters is buggy"=tob, "%F inside clusters is buggy"=tib/ti*100, "%F outside clusters is buggy"=tob/to*100))
		stats <- rbind(stats, stat)
	}

	rownames(stats) <- basename(filenames)
	write.csv(stats, file="stat.csv")
}

RQ1 <- function(filenames, trace=FALSE, parallel=FALSE)
{
	df.RQ1 <- NULL
	for (filename in filenames) {
		cat("RQ1 for ", filename, "\n")

		data <- read.csv(file=filename)

		cnames <- c("C.size.ALLD", "C.ties.ALLD")
		memb <- data$membership.ALLD ### numeric vector giving the cluster id to which each vertex belongs
		tblf <- table(memb)
		memb.cs <- as.numeric(names(which(tblf>=2))) ### which cluster id has more than 1 nodes
		
		df.cluster <- NULL
		for (memb.c in memb.cs) {
			index <- which(memb==memb.c)
			bug_density <- sum(data[index, "BUG"])/length(index)
			df.cluster <- rbind(df.cluster, unlist(c("N"=length(memb.cs), data[index[1], cnames], "bugrate"=bug_density)))
		}

		test.res <- NULL
		for (aname in cnames) {
			ct <- cor.test(df.cluster[, aname], df.cluster[, "bugrate"], method="spearman")
			test.res <- c(test.res, unlist(c(ct$estimate, "p"=ct$p.value)))
		}
		test.res <- round(test.res, 3)
		df.RQ1 <- rbind(df.RQ1, c(length(memb.cs), test.res))
	}

	rownames(df.RQ1) <- basename(filenames)
	write.csv(df.RQ1, file="RQ1.csv")
}

RQ2 <- function(filenames, trace=FALSE, parallel=FALSE)
{
	df.RQ2 <- NULL
	for (filename in filenames) {
		cat("RQ2 for ", filename, "\n")

		data <- read.csv(file=filename)

		TF <- nrow(data) ### 函数总数
		TL <- sum(data$SLOC) ### 代码行总数

		tb <- sum(data$BUG) ### 缺陷函数总数
		tf <- TF - tb ### 无缺陷函数总数

		ti <- sum(data$C.in.ALLD)   ### 簇内函数总数
		to <- TF - ti ### 簇外函数总数

		tib <- sum(data$BUG*data$C.in.ALLD) ### 在簇内并且有缺陷
		tob <- sum(data$BUG*(1-data$C.in.ALLD)) ### 在簇外并且有缺陷

		tif <- sum((1-data$BUG)*data$C.in.ALLD) ### 在簇内并且没有缺陷
		tof <- sum((1-data$BUG)*(1-data$C.in.ALLD)) ### 在簇外并且没有缺陷

		ft <- fisher.test(x=rbind(c(tib, tif), c(tob, tof)), alternative="two.sided") ### fisher test 比较缺陷密度
		pv <- ft$p.value

		p <- tib/ti ### 簇中缺陷比例
		q <- tob/to ### 簇外缺陷比例
		OR <- (p/(1-p))/(q/(1-q))

		tbl <- table(data$membership.ALLD)
		cnum <- length(which(tbl>=2)) ### 簇个数
		lsize <- tbl[which.max(tbl)] ### 最大簇的size

		stat <- unlist(c("SLOC"=TL, "#functions"=TF, "#F is buggy"=tb, "%F is buggy"=tb/TF*100, "#clusters"=cnum, "Max cluster size"=lsize, "#F inside clusters"=ti, "%F inside clusters"=ti/TF*100, "#F outside clusters"=to, "%F outside clusters"=to/TF*100, "#F inside clusters is buggy"=tib, "#F outside clusters is buggy"=tob, "%F inside clusters is buggy"=tib/ti*100, "%F outside clusters is buggy"=tob/to*100, "Fisher's Exact Test(p value)"=pv, "OR"=OR))
		
		df.RQ2 <- rbind(df.RQ2, c("inside"=p*100, "outside"=q*100, "Fisher's Exact Test(p value)"=pv, "OR"=OR))
	}

	df.RQ2[, 3] <- p.adjust(df.RQ2[, 3], method="bonferroni")
	rownames(df.RQ2) <- basename(filenames); df.RQ2 <- round(df.RQ2, 3)
	write.csv(df.RQ2, file="RQ2.csv")
}

RQ3 <- function(filenames, trace=FALSE)
{
	### univariate logistic regression analysis
	UnivariateLogisticAna <- function(model="lr", family=binomial, data, yname, xnames, method="normal", effort.name="LOC", number.name=yname, remove.flag=FALSE, trace=FALSE)
	{
		ulr.results <- data.frame(matrix(NA, nrow=length(xnames), ncol=10))
		rownames(ulr.results) <- xnames
		colnames(ulr.results) <- c("Constant", "Coefficient", "p", "OR", "Std.err")

		if (remove.flag) {
			data <- RemoveConfoundingEffect(data, ynames=xnames, confeff.name=effort.name)
		}

		for (i in 1:length(xnames)) {
			xname <- xnames[i]
			if (trace) { cat("Independent name:", xname, "\n") } ### if (mode(data[, xname])=="complex") { next }

			tdata <- data
			fm <- paste(yname, xname, sep="~") ### ("%s ~ %s", yname, xname)
			ul <- glm(as.formula(fm), data=tdata, family=binomial)

			ul.sum <- summary(ul)
			pre <- predict(ul, data, type="response")
			if (tolower(method)=="divide") {
				valid.dt <- data.frame(REL=data[, yname], LOC=data[, effort.name], NUM=data[, number.name], PRE=pre/data[, effort.name])
			} else {
				valid.dt <- data.frame(REL=data[, yname], LOC=data[, effort.name], NUM=data[, number.name], PRE=pre)
			}

			if (is.na(ul$coefficients[xname])) {
				ulr.results[i, "Constant"] <- ul$coefficients["(Intercept)"]
			} else {
				ulr.results[i, "Constant"]    <- ul.sum$coefficients["(Intercept)", "Estimate"]
				ulr.results[i, "Coefficient"] <- ul.sum$coefficients[xname, "Estimate"]
				ulr.results[i, "Std.err"] <- ul.sum$coefficients[xname, "Std. Error"]
				ulr.results[i, "p"] <- ul.sum$coefficients[xname, "Pr(>|z|)"]

				coef <- ul$coefficients[xname] # coefficient
				asd  <- sd(tdata[, xname])
				delta_OR <- exp(coef * asd)
				ulr.results[i, "OR"] <- delta_OR # delta_OR
			}
		}

		ulr.results[, "p-adj"] <- p.adjust(p=ulr.results[, "p"], method="BH")
		return(ulr.results)
	}

	c.xnames <- c("C.betweenness.sub.ALLD","C.centr_betw.sub.ALLD","C.centr_clo.out.sub.ALLD","C.centr_degree.out.sub.ALLD","C.centr_eigen.sub.ALLD","C.closeness.out.sub.ALLD","C.constraint.sub.ALLD","C.degree.out.sub.ALLD","C.eccentricity.out.sub.ALLD","C.page_rank.sub.ALLD") ### all independent variable names

	df.RQ3 <- df.RQ3.p <- NULL
	for (filename in filenames) {
		cat("RQ3 for", filename, "\n")

		data <- read.csv(file=filename)

		ulr <- UnivariateLogisticAna(data=data[data$C.in.ALLD==1, ], yname="BUG", xnames=c.xnames, effort.name="SLOC", trace=FALSE); gc()
		ulr[is.na(ulr)] <- -1
		ulr <- round(ulr, 3)
		df.RQ3 <- cbind(df.RQ3, ulr$OR)
		df.RQ3.p <- cbind(df.RQ3.p, ulr$p)
		rownames(df.RQ3) <- rownames(ulr)
		rownames(df.RQ3.p) <- rownames(ulr)
	}

	for (i in nrow(df.RQ3.p)) {
		df.RQ3.p[i, ] <- p.adjust(df.RQ3.p[i, ], method = "bonferroni")
	}
	colnames(df.RQ3) <- paste(basename(filenames), "Delta OR")
	colnames(df.RQ3.p) <- paste(basename(filenames), "pvalue")
	
	df.RQ3 <- cbind(df.RQ3, df.RQ3.p)
	write.csv(df.RQ3, file="RQ3.csv")
}

RQ4 <- function(filenames, model="lr.2phase", rule="aic", direction="forward", remove.flag=FALSE, trace=FALSE)
{
	B.xnames <- c("SLOC","FANIN","FANOUT","NPATH","Cyclomatic","CyclomaticModified","CyclomaticStrict","Essential","Knots","MaxNesting","MaxEssentialKnots","MinEssentialKnots","n1","n2","N1","N2","added","changed","deleted","E.nWeakComp.out.ALLD","E.Broker.out.ALLD","E.EgoBetw.out.ALLD","E.constraint.out.ALLD","E.Size.out.ALLD","E.Ties.out.ALLD","E.nBroker.out.ALLD","E.nEgoBetw.out.ALLD","E.Pairs.out.ALLD","E.Density.out.ALLD","E.pWeakComp.out.ALLD","E.2StepReach.out.ALLD","E.ReachEffic.out.ALLD","E.effsize.out.ALLD","E.efficiency.out.ALLD","G.degree.out.ALLD","G.dwReach.out.ALLD","G.closeness.out.ALLD","G.eigen_centr.ALLD","G.betweenness.ALLD","G.power_centr.ALLD")
	C.xnames <- c("C.betweenness.sub.ALLD","C.centr_betw.sub.ALLD","C.centr_clo.out.sub.ALLD","C.centr_degree.out.sub.ALLD","C.centr_eigen.sub.ALLD","C.closeness.out.sub.ALLD","C.constraint.sub.ALLD","C.degree.out.sub.ALLD","C.eccentricity.out.sub.ALLD","C.page_rank.sub.ALLD")

	### all independent variable names
	xnames <- c(B.xnames, C.xnames)

	### all investigated models
	feature.names <- list("B"=B.xnames, "B+C"=c(B.xnames, C.xnames))

	totalRuns  <- 30
	totalFolds <- 3
	
	### this should be setted by user
	number.name <- "BUG"
	cutoffs <- 0.2

	for (filename in filenames) {
		main(filename=filename, xnames=xnames, number.name=number.name, feature.names=feature.names, model=model, rule=rule, direction=direction, remove.flag=remove.flag, cutoffs=cutoffs, totalFolds=totalFolds, totalRuns=totalRuns, trace=trace)	
	}	
}

projects <- c("bash-3.2", "gcc-4.0.0", "gimp-2.8.0", "glibc-2.1.1", "gstreamer-1.0.0")

filenames <- file.path("input", paste(projects, "-combine.csv", sep=""))

Des(filenames=filenames)
RQ1(filenames=filenames)
RQ2(filenames=filenames)
RQ3(filenames=filenames)
RQ4(filenames=filenames)

for (filename in filenames) {	
	plot.single.effort.ranking(filename)
	plot.single.effort.classification(filename)
}

table.effort.ranking(filenames)
table.effort.classification(filenames)
