Difference between revisions of "Object-oriented programming in Opasnet"

library(OpasnetBaseUtils)
n <- 5

#################### Defines the S4 class "ovariable" which is the basic building block in open assessments.
setClass(
	"ovariable", 
	representation(
		sample = "data.frame", 
		data = "data.frame", 
		marginal = "vector", 
		formula = "function", 
		dependencies = "list"
	)
)

#################### Defines the S4 class "oassessment" which is the object type for open assessments.
setClass(
	"oassessment", 
	representation(
		names         = "data.frame", 
		decisions     = "data.frame",
		probabilities = "data.frame",
		stakeholders  = "data.frame",
		vars          = "list"
	)
)

########## Arithmetic operations of ovariables: first they are merged by index columns,
### then the operation is performed for the Result.x and Result.y columns.
### If one of the expressions is numeric, it is first transformed to ovariable.
setMethod(
	f = "Ops", 
	signature = signature(e1 = "ovariable", e2 = "ovariable"), 
	definition = function(e1, e2) {
		out <- merge(e1, e2)@sample
		colnames(out) <- gsub(".x", "", colnames(out))
		out$Result <- callGeneric(out$Result, out$Result.y)
		if(!is.null(out$Unit.y)) {out$Unit <- paste(out$Unit, "|(", out$Unit.y, ")", sep= "")}
		e1@sample <- out[, !colnames(out) %in% c("Result.y", "Unit.y")]
		return(e1)
	}
)

setMethod(
	f = "Ops", 
	signature = signature(e1 = "ovariable", e2 = "numeric"), 
	definition = function(e1, e2) {
		e2 <- make.ovariable(e2)
		e1 <- callGeneric(e1, e2)
		return(e1)
	}
)

setMethod(
	f = "Ops", 
	signature = signature(e1 = "numeric", e2 = "ovariable"), 
	definition = function(e1, e2) {
		e1 <- make.ovariable(e1)
		e1 <- callGeneric(e2, e1)
		return(e1)
	}
)

#################### Math defines basic mathematical operations (log, exp, abs, ...) for ovariables
setMethod(
	f = "Math", 
	signature = signature(x = "ovariable"), 
	definition = function(x) {
		x@sample$Result <- callGeneric(x@sample$Result)
		return(x)
	}
)

############ merge of ovariables merges the 'sample' slot by index columns except 'Unit'.
setMethod(f = "merge", 
	signature = signature(x = "ovariable", y = "ovariable"),
	definition = function(x, y) {
		test <- intersect(colnames(x@sample[x@marginal]), colnames(y@sample[y@marginal]))
		x@sample <- merge(x@sample, y@sample, by = test, all = TRUE)
		return(x)
	}
)

setMethod(f = "merge", 
	signature = signature(x = "ovariable", y = "numeric"),
	definition = function(x, y) {
		y <- make.ovariable(y)
		return(merge(x, y))
	}
)

setMethod(f = "merge", 
	signature = signature(x = "numeric", y = "ovariable"),
	definition = function(x, y) {
		x <- make.ovariable(x)
		return(merge(y, x))
	}
)

############ tapply of ovariables applies a function to each cell of a ragged array, that is to each (non-empty) group of
############ values given by a unique combination of the levels of certain factors.
### parameters (other parameters are as in generic tapply):
### X an ovariable

setMethod(f = "tapply", 
	signature = signature(X = "ovariable"),
	definition = function(X, INDEX, FUN = NULL, ..., simplify = TRUE) {
		out <- as.data.frame(as.table(tapply(X@sample$Result, INDEX, FUN, ..., simplify = TRUE)))
		colnames(out)[colnames(out) == "Freq"] <- "Result"
		X@sample <- out
		return(X)
	}
)

############################ orbind combines two ovariables using clever rbind

orbind <- function(x, y) {
		if(class(x) == "ovariable") {xsample <- x@sample} else {xsample <- x}
		if(class(y) == "ovariable") {ysample <- y@sample} else {ysample <- y}
		cols <- setdiff(colnames(ysample), colnames(xsample)) # Take all columns that do not exist in x and add them.
		col <- as.data.frame(array("NA", dim = c(1, length(cols))))
		colnames(col) <- cols
		if("Unit" %in% cols) {col[, "Unit"] <- "?"}
		temp <- cbind(xsample, col)

		cols <- setdiff(colnames(xsample), colnames(ysample)) # Take all columns that do not exist in y and add them.
		col <- as.data.frame(array("NA", dim = c(1, length(cols))))
		colnames(col) <- cols
		if("Unit" %in% cols) {col[, "Unit"] <- "?"}
		xsample <- rbind(temp, cbind(ysample, col)) # Combine x and y with rbind.
		return(xsample)
#Should this be made S4 function for ovariables? Then it could be named simply rbind.
	}


################# plot diagrams about ovariable data

setMethod(
	f = "plot",
	signature = signature(x = "ovariable"),
	definition = function(x) {
		x <- x@sample
		plot(x$Source, x$Result, ylab = x[x$Source == "Data", "Unit"][1])
	}
)

################# plot diagrams about oassessment data

setMethod(
	f = "plot",
	signature = signature(x = "oassessment"),
	definition = function(x) {
		for(i in names(x@vars)) {
			y<- x@vars[[i]]@sample
			plot(y$Source, y$Result, ylab = y[y$Source == "Data", "Unit"][1], main = i)
		}
	}
)

######################### print ovariable contents

setMethod(
	f = "print",
	signature = signature(x = "oassessment"),
	definition = function(x) {
		cat("Names\n")
		print(xtable(x@names), type = 'html')
		cat("Decisions\n")
		print(xtable(x@decisions), type = 'html')
		cat("Stakeholders\n")
		print(xtable(x@stakeholders), type = 'html')
		cat("Probabilities\n")
		print(xtable(x@probabilities), type = 'html')

		cat("\n\nThe list of variables in this assessment.\n")
		for(k in 1:length(x@vars)) {
			cat("<b>Variable", names(x@vars)[k], "</b> <p>\n\n")
			print(x@vars[[k]])
		}
	}
)


setMethod(
	f = "print",
	signature = signature(x = "ovariable"),
	definition = function(x) {
		cat("Sample\n")
		print(xtable(x@sample[x@sample$Iter %in% 1, ]), type = 'html')
		cat("Formula\n")
		print(x@formula)
		cat("Dependencies\n")
		print(x@dependencies)
	}
)

###################### fetch downloads a variable.

fetch <- function(x) {
	out <- tidy(op_baseGetData("opasnet_base", x), direction = "wide")
	return(out)
}


########### init.assessment creates S4 assessment from names data.frame, including decisions, stakeholders, probabilities, and variables. 
########### NOTE! You must include the formula code from each variable page, otherwise formulas and dependencies are not updated. 
########### Parameters:
## names: a data.frame that has the structure of oassessment@name (Columns: Name, Identifier, Direction, Result)
init.assessment <- function(names) 
{
	names         <- fetch(names)
	decisions     <- fetch(names[names$Result == "decisions", "Identifier"])
	stakeholders  <- fetch(names[names$Result == "stakeholders", "Identifier"])
	probabilities <- fetch(names[names$Result == "probabilities", "Identifier"])
	names         <- names[!names$Result %in% c("decisions", "stakeholders", "probabilities"), ]
	vars          <- list()
	for(x in 1:nrow(names)) { # Objects with names as aliases are created and filled with data from Opasnet Base.
		cat("Initialising variable ", as.character(names$Name[x]), ".\n", sep = "")
		ident <- as.character(names$Identifier[x])
		temp <- tidy(op_baseGetData("opasnet_base", ident), direction = as.character(names$Direction[x]))
		temp <- make.ovariable(temp)
		if(exists(paste("formula.", ident, sep = ""))) 
			{temp@formula <- get(paste("formula.", ident, sep = ""))}
		if(exists(paste("dependencies.", ident, sep = "")))
			{temp@dependencies <- get(paste("dependencies.", ident, sep = ""))}
		vars[as.character(names$Result[x])] <- temp
	}
	assessment <- new("oassessment", 
		names         = names, 
		decisions     = decisions,
		stakeholders  = stakeholders,
		probabilities = probabilities,
		vars          = vars
	)

	return(assessment)
# Does it cause problems to use the wide direction? Anyway, there should be exactly one Result column.
## No it doesn't if the Observation columns are such as Unit, Result, Description.
}

########### update updates the sample of an ovariable based on data and function.
setMethod(
	f = "update", 
	signature = "ovariable", 
	definition = function(object) {
		dat  <- data.frame(Source = "Data", interpret(object@data))
		form <- object@formula(object@dependencies)
		if(is.vector(form)) {form <- data.frame(Result = form)}
		if(class(form) == "ovariable") {form <- form@sample}
		form <- data.frame(Source = "Formula", form)
		object@sample <- orbind(dat, form)
		object@marginal <- c(TRUE, object@marginal) # This alone is not enough; orbind must operate with marginals as well.
		return(object)
	}
)


#### This code was taken out of update. It is still needed somewhere but not here.
#		dep <- object@dependencies
#		for(i in 1:length(dep)) {
#			if(class(dep[[i]]) == "ovariable") {
#				dep[[i]] <- dep[[i]]@sample
#			} else {
#				if(length(grep("Op_(en|fi)", dep[[i]])) > 0) {
#					dep[[i]] <- op_baseGetData("opasnet_base", dep[[i]])}
#				else {
#					if(class(dep[[i]]) != "data.frame" & !is.numeric(dep[[i]])) {
#						dep[[i]] <- get(dep[[i]])
#					}
#				}
#			}
#		}


#################### interpret takes a vector and makes a data.frame out of it (to be used in e.g. make.ovariable).
### It also changes abbreviations into probability samples.
interpret <- function(data) {
	sample <- NULL
	if(is.vector(data)) {data <- data.frame(Result = data)}
	if("Iter" %in% colnames(data)) {
		out <- data}
	else {
		test <- !is.na(as.numeric(as.character(data$Result)))
	
		for(i in 1:nrow(data)) {
			if(test[i]) {
				sample <- c(sample, rep(as.numeric(as.character(data[i, "Result"])), n))
			} else {
				samplingguide <- as.numeric(strsplit(gsub(" ", "", data[i, "Result"]), "-")[[1]])
				if(is.na(samplingguide[1]) | is.na(samplingguide[2])) {
					sample <- c(sample, rep(data[i, "Result"], n))
				} else {
					sample <- c(sample, runif(n, samplingguide[1], samplingguide[2]))
				}
			}
		}
	
		out <- as.data.frame(array(1:(n*nrow(data)*(ncol(data)+1)), dim = c(n*nrow(data), ncol(data) + 1)))
		colnames(out) <- c("Iter", colnames(data))
	
		for(i in colnames(data)) {
			out[i] <- rep(data[, i], each = n)
		}
		out$Iter <- 1:n
		out$Result <- sample
		
	}
	return(out)
}

#################################################################################
########### make.ovariable takes a vector or data.frame and makes an ovariable out of it.
make.ovariable <- function(
	data, 
	formula = function(dependencies){return(0)}, 
	dependencies = list(x = 0)
) {
	return(movariable(data, formula, dependencies))
}

setGeneric("make.ovariable") # Makes make.ovariable a generic S4 function.

setMethod(
	f = "make.ovariable",
	signature = signature(data = "data.frame"),
	definition = function(
		data, 
		formula = function(dependencies){return(0)}, 
		dependencies = list(x = 0)
	) {
cat("Data frame\n")
		return(movariable(data, formula, dependencies))
	}
)

setMethod(
	f = "make.ovariable",
	signature = signature(data = "vector"),
	definition = function(
		data, 
		formula = function(dependencies){return(0)}, 
		dependencies = list(x = 0)
	) {
cat("Vector\n")
		data <- data.frame(Result = data)
		return(make.ovariable(data, formula, dependencies))
	}
)

setMethod(
	f = "make.ovariable",
	signature = signature(data = "list"),
	definition = function(
		data, 
		formula = function(dependencies){return(0)}, 
		dependencies = list(x = 0)
	) {
		for(i in 1:length(data)) {
cat("List", i, "\n")
			data[[i]] <- make.ovariable(data[[i]], formula, dependencies)
		}
		return(data)
	}
)

setMethod(
	f = "make.ovariable",
	signature = signature(data = "ovariable"),
	definition = function(
		data, 
		formula = NULL, 
		dependencies = NULL
	) {
cat("ovariable\n")
		if(is.null(formula)) {formula <- data@formula}
		if(is.null(dependencies)) {dependencies <- data@dependencies}
		return(movariable(data@data, formula, dependencies))
	}
)


########### movariable takes a data.frame, a function, and a list and makes an ovariable out of them. It is a 
#####subfunction of make.ovariable and prevents infinite recursion of S4 methods.
movariable <- function(
	data, 
	formula, 
	dependencies
) {
	sample <- interpret(data)
	out <- new("ovariable", 
		sample = sample, 
		data = data,
		marginal = ifelse(colnames(sample) %in% c("Result", "Unit"), FALSE, TRUE), 
		formula = formula, 
		dependencies = dependencies)
	out <- update(out)
	return(out)
}

a <- interpret("1-2")
print(a)
a <- a * 3
print(a)

library(xtable)	
a <- data.frame(a = c("a", "b", "c", "d"), Obs = 1, Unit = "ug /m^3", Result = c(1, "3 - 4", 5, "7 - 8") )
print(xtable(a), type = 'html')
print(xtable(tidy(op_baseGetData("opasnet_base", "Op_en5483"))), type = 'html')

###########   Defining variable d

d <- make.ovariable(
	data = a, 
	formula = function(dependencies) {
		x <- tidy(dependencies$x)
		x <- x[x$Observation == "Per.person", ]
		x <- make.ovariable(x)
		x*10 + dependencies$b
	}, 
	dependencies = list(x = "Op_en5483", b = 1)
)

############## Calculating with d

d <- update(d)
plot(d)
print(xtable(d@sample[d@sample$Iter < 4, ]), type = 'html')

############################################################

cat("Initiation successful. Now starting the model.\n")
library(xtable)

make.ovariable("1-2") * 5

make.ovariable(data.frame(Result = "1-2"))


dependencies.Op_en5675 <- list(
	exposure = "Op_en5674", # formula.Op_en5674(dependencies.Op_en5674), # Training exposure
	erf = data.frame(Unit = "RR per ug/m3", Result = 1.5),
	population = population,
	background = 100 / 100000 # cases per 100000 person-years
)

exposure <- op_baseGetData("opasnet_base", dependencies.Op_en5675$exposure)
colnames(exposure)[colnames(exposure) == "obs.1"] <- "Obs"
exposure <- tidy(exposure, direction = "wide")
colnames(exposure)[colnames(exposure) == "result"] <- "Result"
exposure <- make.ovariable(exposure)
print(exposure)

formula.Op_en5675 <- function(x) {
	population <- make.ovariable(x$population)
	background <- make.ovariable(x$background)
	exposure <- op_baseGetData("opasnet_base", x$exposure)
	colnames(exposure)[colnames(exposure) == "obs.1"] <- "Obs"
	exposure <- tidy(exposure, direction = "wide")
	colnames(exposure)[colnames(exposure) == "result"] <- "Result"
	exposure <- make.ovariable(exposure)
	erf <- make.ovariable(x$erf)
	cases <- population * background * exp(exposure * log(erf))
	return(cases)
}

formula.Op_en5675(dependencies.Op_en5675)

make.ovariable(dependencies.Op_en5675)

out <- make.ovariable(
	data = "0 - 100000",
	formula = formula.Op_en5675,
	dependencies = dependencies.Op_en5675)

print(out)
out <- update(out)
print(out)
plot(out)