# Multivariate random effects wraper file for demonstration purposes
# see readme.txt for overview
# load up all the associated functions
source("MVRA-functions.r")

# the data is from JoAnn Buscaglia of the FBI laboratory and consists of 310 measurements
# on 62 panes of glass, five replicates on each.
# the varibles are:
# Si - silicon
# K  - potassium
# Ca - calcium
# Fe - Iron
# Group - not relevant here
# log(CaK) - log10 of calcuim to potassium
# log(CaSi) - log10 of calcium to silicon
# log(CaFe) - log10 of calcium to iron
# window - which of the 62 panes the measurement was from
dat <-  read.table("glass-data.txt", header=T)

# define the population
population <- dat

# define the continuous variables to be used - in this case it is the log ratios
variables <- c(6,7,8)

# define a vector to put the output in
output <- rep(0, 4)

# define the variable which is to be used as the grouping factor
grouping.variable <- 9

# group is a vector 
group <- population[,grouping.variable]

# sort out the unique identifiers for the windows
windows <- unique(group)
# and how many windows we are dealing with
no.windows <- length(windows)

# we have five measurements on each pane so to generate a set of measurements for a
# control pane, and a set of measurements for a recovered pane we simply divide the
# five into two groups of three, with one shared measurement. This way we know we have
# two sets of three measurements which come from the same pane
control.items <- c(1,2,3)
recovered.items <- c(3,4,5)

# calculate the within and between pane co-variences and means for each pane
UC <- calculate.UC(population, variables, grouping.variable)
# between pane
C <- UC$C
# within pane
U <- UC$U
# means
group.means <- UC$group.means

# which pane acts as the control window
# please modify at your pleasure
control.window <- 10

# extract the data for the control.window
base <- population[which(population[,grouping.variable] == control.window),]
# control is a matrix of measurement data from the control window
control <- base[control.items,]

	# calculate a LR for each pane in the data
	for(recovered.window in 1:no.windows)
#	for(recovered.window in 1:3)

		{
		# extract a suitable subset of the data to be the recovered sample
		base <- population[which(dat[,grouping.variable] == recovered.window),]
		recovered <- base[recovered.items,]

		# calculate the means of variables for control and recovered	
		Z <- y1.y2.means(control, variables, recovered, variables)
		N.control <- Z$Nc
		N.recovered <- Z$Nr

		# calculate a suitable window width for the density estimation
		h.opt <- calculate.h.opt(no.windows, length(variables))

		# calculate the likelihood ratio using kernel density approach to between sources
		LR.MVRAD <- MVRA.density(Z$control.mean, Z$recovered.mean, N.control, N.recovered, variables, U, C, group.means, h.opt)

		# calculate the likelihood ratio using normal assumptions to between sources
		recovered.minus.suspect <- Z$recovered.mean - Z$control.mean
		LR.norm <- MVRA.normal(U, C, recovered.minus.suspect, N.control, N.recovered)

		# assign the answers
		output.row <- c(control.window, recovered.window, LR.MVRAD, LR.norm)
		output <- rbind(output, output.row)
		}

colnames(output) <- c("control", "recovered", "MVRAD", "LR-Normal")


# get rid of the leading row in the output matrix
output <- output[2:nrow(output),]
rownames(output) <- seq(1:nrow(output))
print(round(output, 4))