#!/usr/bin/env RScript

#require(ggplot2)
require(plotly)
require(reshape2)

data <- data.frame(
	y=c(13, 5, 22),
	x1=c(4, 13, 15),
	x2=c(35, 24, 67))
#general
a_ <- matrix(c(
3,6,8,9,3,0,
7,9,3,6,0,1,
7,9,6,5,4,7),
ncol=3)
#continuous
a <-matrix(c(
5,5,9,9,3,3,
9,3,5,3,5,9,
3,9,3,5,9,5),
ncol=3)
#convex
aq <- matrix(c(
2,2,6,6,7,7,
6,7,2,7,2,6,
7,6,7,2,6,2),
ncol=3)
#one
az <-matrix(c(
1,1,1,1,1,1,
1,1,1,1,1,1,
1,1,1,1,1,1),
ncol=3)

L <- function(beta, data, a) {
	betaperm <- logical(6)
	permutations <- (matrix(c(
	1, 1, 2, 2, 3, 3,
	2, 3, 1, 3, 1, 2,
	3, 2, 3, 1, 2, 1),
	ncol=3))
	for (pi in 1:6) {
# 		print((data$y[permutations[pi,1]] - (beta[1]*data$x1[permutations[pi,1]] + beta[2]*data$x2[permutations[pi,1]])))
# 		print((data$y[permutations[pi,2]] - (beta[1]*data$x1[permutations[pi,2]] + beta[2]*data$x2[permutations[pi,2]])))
# 		print((data$y[permutations[pi,3]] - (beta[1]*data$x1[permutations[pi,3]] + beta[2]*data$x2[permutations[pi,3]])))
# 		print("---")
# 		print(paste(permutations[pi,1], permutations[pi,2], permutations[pi,3]))
# 		print("===")
		#browser()
		betaperm[pi] <-
			((data$y[permutations[pi,1]] - (beta[1]*data$x1[permutations[pi,1]] + beta[2]*data$x2[permutations[pi,1]])) <=
			(data$y[permutations[pi,2]] - (beta[1]*data$x1[permutations[pi,2]] + beta[2]*data$x2[permutations[pi,2]]))) &&
			((data$y[permutations[pi,2]] - (beta[1]*data$x1[permutations[pi,2]] + beta[2]*data$x2[permutations[pi,2]])) <=
			(data$y[permutations[pi,3]] - (beta[1]*data$x1[permutations[pi,3]] + beta[2]*data$x2[permutations[pi,3]])))
	}
	aipi <- match(TRUE, betaperm)
# 	browser()
	Lres <- a[aipi, 1]*(data$y[1]-(beta[1]*data$x1[1]+beta[2]*data$x2[1])) +
			a[aipi, 2]*(data$y[2]-(beta[1]*data$x1[2]+beta[2]*data$x2[2])) +
			a[aipi, 3]*(data$y[3]-(beta[1]*data$x1[3]+beta[2]*data$x2[3]))
# 	print(betaperm)
# 	print(aipi)
# 	print(paste(a[aipi, 1], (data$y[1]-(beta[1]*data$x1[1]+beta[2]*data$x2[1]))))
# 	print(paste(a[aipi, 2], (data$y[2]-(beta[1]*data$x1[2]+beta[2]*data$x2[2]))))
# 	print(paste(a[aipi, 3], (data$y[3]-(beta[1]*data$x1[3]+beta[2]*data$x2[3]))))
# 	print((data$y[1]-(beta[1]*data$x1[1]+beta[2]*data$x2[1])) +
# 			(data$y[2]-(beta[1]*data$x1[2]+beta[2]*data$x2[2])) +
# 			(data$y[3]-(beta[1]*data$x1[3]+beta[2]*data$x2[3])))
	return(list(L=Lres, aipi=aipi))
}

plotL <- function(beta1, beta2) {
	Lperm <- sapply(beta1, function(x1){
		sapply(beta2, function(x2) {
			L(c(x1, x2), data, a)
		})
	})
	Lbeta <- Lperm[seq(1,nrow(Lperm)-1,2),]
	Cells <- Lperm[seq(2,nrow(Lperm),2),]
# 	print(Lbeta)
	colnames(Lbeta) <- as.character(beta1)
	rownames(Lbeta) <- as.character(beta2)
	colnames(Cells) <- as.character(beta1)
	rownames(Cells) <- as.character(beta2)
#	dimnames(Lbeta) <- list(as.character(beta1), as.character(beta2))
	dataWidePlot <- cbind(expand.grid(dimnames(Lbeta)), value = unlist(as.vector(Lbeta)))
	colnames(dataWidePlot) <- c("beta1", "beta2", "L")
	ggData <- ggplot(dataWidePlot, aes(x=beta1, y=beta2)) +
		geom_tile(aes(fill=L))

	cellWidePlot <- cbind(expand.grid(dimnames(Lbeta)), value = unlist(as.vector(Cells)))
	colnames(cellWidePlot) <- c("beta1", "beta2", "C")
	cellWidePlot$C <- as.factor(cellWidePlot$C)
	ggCell <- ggplot(cellWidePlot, aes(x=beta1, y=beta2)) +
		geom_tile(aes(fill=C))

	return(list(ggData, ggCell))
# 	return(plot_ly() %>% add_surface(x=~beta1, y=~beta2, z=~Lbeta))
}

# plot_ly() %>% add_surface(x=c(1,2), y=c(1,2), z=matrix(c(1,2,3,4), ncol=2))
# plot_ly(data, x= ~x1, y= ~y)
# export(plotL(seq(-1, 2, 0.1), seq(-0.5, 1.5, 0.1)), "l.png")
renPlot <- plotL(seq(-1, 1, 0.03), seq(-1, 1, 0.03))
L(c(4,6), data, a)
L(c(4,4.1), data, a)
L(c(4.1,4), data, a)