# libVMR 3.01
#
# Copyright (C) 2014,  Yury V. Reshetov
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>. 
#
#
# R code ported from JAVA https://sourceforge.net/projects/libvmr/
# original JAVA code by Yury V. Reshetov


########################
#Parameters:
 
trainFile     <- "C:\\train.csv"
fronttestFile <- "C:\\test.csv"

iterPerColumn <- 1000000

kernelTrickMode <- 0 #default 0
#0 - Automatic (<=11 Big; <=45 Oairs; >45 Small)
#1 - Big
#2 - Pairs
#3 - Small

seed <- 7777

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



#сложное полиминальное ядерное преобразование
#@author Yury V. Reshetov
#@version 3.01
BigKernelMachine <- function(samples) {
	count <- 1;
	for (i in 1:ncol(samples)) {
		count <- count * 2;
	}
	result <- matrix(0, nrow = nrow(samples), ncol = count)
	for (i in i:nrow(samples)) {
		for (j in 1:count) {
			pp <- 1;
			z <- j-1; #original JAVA code was "z=j". Changed it to "j-1" to have same result, as R array indexing start with 1
			for (n in 1:ncol(samples)) {
				if (z %% 2 == 1) { 
					pp = pp * samples[i,n]
				}
				z <- floor(z / 2)
			}
			result[i,j] = pp;               
		}
	}
	return(result)
}


#Парное ядерное преобразование
#@author Yury V. Reshetov
#@version 3.01
PairsKernelMachine <- function(samples){
	tempresult <- cbind(matrix(1, ncol=1, nrow=nrow(samples)), samples)
	
	result <- matrix(0, nrow = nrow(samples), ncol = floor((ncol(samples) + 1) * (ncol(samples) + 2) / 2))
	for (k in 1:nrow(samples)) {
		n = 1;
		for (i in 1:(ncol(samples)+1)) {
			for (j in 1:(ncol(samples)+1)) {
				#Проверка наличия дубликатов
				if (j >= i) {
					if ((i > 1) && (i == j)) {
						result[k, n] = 2 * tempresult[k, i] * tempresult[k, j] - 1
					} else {
						result[k, n] = tempresult[k, i] * tempresult[k, j]
					}
					n <- n + 1
				}
			}
		}
	}
	return(result)
}




#простое ядерное преобразование
#@author Yury V. Reshetov
#@version 3.01
SimpleKernelMachine <- function(samples) {
	return (cbind(matrix(1, ncol=1, nrow=nrow(samples)), samples));
}




#Преобразование переменных из диапазона от минимального значения до
#максимального включительно к диапазону от -1 до 1 включительно
convertToRange <- function(values) {
	minVal <- min(values)
	maxVal <- max(values)
	if(minVal==maxVal){
		return(rep(0,length(values)))
	}else{
		return(2*(values-minVal)/(maxVal-minVal)-1)
	}
}






#Инвариантное преобразование выборки в результате чего два бинарных класса
#после преобразования становятся одним общим классом
getInvariant <- function(samples, y) {
	samples <- as.matrix(samples)
	result <- matrix(0, ncol=ncol(samples), nrow=nrow(samples))
	for (i in 1:nrow(samples)) {
		result[i,] <- samples[i,] / y[i]
	}
	return(result)
}





#Replace missing values
#Normalize
#Kernel trick
TransformSamples <- function(samples){
	samples[is.na(samples) || is.null(samples)] <- 0

	for(i in 1:ncol(samples)){
		samples[, i] <- convertToRange(samples[, i])
	}
	
	if(kernelTrickMode == 0){
		if (ncol(samples) <= 11) {
			samples <- BigKernelMachine(samples)
		}else if(ncol(samples) <= 45){
			samples <- PairsKernelMachine(samples)
		}else{
			samples <- SimpleKernelMachine(samples)
		}
	}else if(kernelTrickMode == 1){
		samples <- BigKernelMachine(samples)
	}else if(kernelTrickMode == 2){
		samples <- PairsKernelMachine(samples)
	}else if(kernelTrickMode == 3){
		samples <- SimpleKernelMachine(samples)
	}
	return(samples)
}




#Обучение искусственного нейрона по методу Брауна, Робинсон, Решетова
#
#author Yury V. Reshetov
#version 3.01
BrownRobinsonReshetovAlgorithmTrain <- function(x, y){
	x <- getInvariant(x, y)
	x <- cbind(x, -x)
	#Количество строк
	rows <- nrow(x)
	#Количество столбцов
	columns <- ncol(x)
	#Счётчик для игрока по строкам - оппонента
	counterrows <- rep(0, rows)
	#Счётчик для игрока по столбцам - генератор гипотез
	countercolumns <- rep(0, columns)
	#Накопитель для строк
	rowsadder <- rep(0, rows)
	#Накопитель для столбцов
	collsadder <- rep(0, columns)
	#Индекс выбранной строки
	selectedrow = sample(1:rows, 1)
	#Индекс выбранного столбца
	seletctedcolumn <- 1;
	#Количество холостых итераций
	delta <- floor(columns / 2);
	if (delta <= 1) {
		delta <- 2;
	}
	for (u in 1:columns) {
		if (u == delta) {
			#После завершения холостых итераций очищаем счётчик столбцов
			countercolumns[] <- 0
		}
		for (t in 1:iterPerColumn) {
			for (j in 1:columns) {
				collsadder[j] <- collsadder[j] + x[selectedrow, j]
			}
			seletctedcolumn <- 1;
			#Ищем лучший столбец
			for (j in 2:columns) {
				if ((collsadder[j] == collsadder[seletctedcolumn]) && (countercolumns[j] < countercolumns[seletctedcolumn])) {
					seletctedcolumn <- j;
				}
				if (collsadder[j] > collsadder[seletctedcolumn]) {
					seletctedcolumn <- j;
				}
			}
			#Улучшаем гипотезу
			countercolumns[seletctedcolumn] <- countercolumns[seletctedcolumn] + 1;
			#Ищем пример, наименее соответствующий гипотезе
			for (i in 1:rows) {
				rowsadder[i] <- rowsadder[i] + x[i, seletctedcolumn];
			}
			selectedrow <- 1;
			for (i in 2:rows) {
				if ((rowsadder[i] == rowsadder[selectedrow]) && (counterrows[i] < counterrows[selectedrow])) {
					selectedrow <- i;
				}
				if (rowsadder[i] < rowsadder[selectedrow]) {
					selectedrow <- i;
				}
			}
			#Увеличиваем счётчик для найденного примера
			counterrows[selectedrow] <- counterrows[selectedrow] + 1;
		}
		cat("Progress:", 100*u/columns,"%\n")
	}

	count <- floor(columns / 2);
	#Результаты подсчёта для столбцов
	result <- rep(0, count)
	#Весовые коэффициенты искусственного нейрона
	returnWeights <- rep(0, count)
	#Максимальное значение для нормировки весовых коэффициентов
	maxVal <- 0;
	for (i in 1:count) {
		#Вычисляем разницу счётчика i-го столбца между ЗА и ПРОТИВ
		result[i] <- countercolumns[i] - countercolumns[i + count];
		#Присваиваем разницу i-му весовому коэффициенту искусственного
		#нейрона
		returnWeights[i] <- result[i];
		#Если найдено максимальное значение для нормировки
		if (abs(returnWeights[i]) > maxVal) {
			#Запоминаем это значение
			maxVal <- abs(returnWeights[i]);
		}
	}

	#this.constantnonzero = Math.abs(result[0]) > 0;

	for(i in 1:count) {
		#Нормируем весовые коэффициенты искусственного нейрона
		returnWeights[i] <- returnWeights[i] / maxVal;
		#Записываем результаты в текстовом виде
		#Если весовой коэффициент не нулевой
		if (abs(result[i]) > 0) {
			#То указываем его
			#variables[i] = returnWeights[i] + variables[i]
		} else {
			#Обнулённые весовые коэффициенты нам не нужны
			#variables[i] = "";
			#accounting.clearPredictor(i)
		}
	}
	return (returnWeights)
}


predict <- function(VMRweights, newdata){
	result <- rep(0, nrow(newdata))
	if(length(VMRweights) == ncol(newdata)){
		for(i in 1:nrow(newdata)){
			result[i] <- sum(VMRweights * newdata[i,])
			if(result[i]>1){
				result[i] = 1
			}else if(result[i]<-1){
				result[i] <- -1
			}
		}
	}
	return(result)
}



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

startTime <- Sys.time()

set.seed(seed)


trainData <- as.matrix(read.csv(trainFile, header=TRUE, sep=";", na.strings=c(".", "NA", "", "?"), strip.white=TRUE, encoding="UTF-8", check.names=FALSE))
fronttestData <- as.matrix(read.csv(fronttestFile, header=TRUE, sep=";", na.strings=c(".", "NA", "", "?"), strip.white=TRUE, encoding="UTF-8", check.names=FALSE))


predictorNames <- colnames(trainData)[1:(ncol(trainData)-1)]
targetName <- tail(colnames(trainData),1)


trainX <- trainData[,predictorNames]
trainX <- TransformSamples(trainData[,predictorNames]) #Normalize predictors and do kernel trick
trainY <- convertToRange(trainData[,targetName])       #Scale target to [-1;1]

fronttestX <- fronttestData[,predictorNames]
fronttestX <- TransformSamples(fronttestData[,predictorNames]) #Normalize predictors and do kernel trick
fronttestY <- convertToRange(fronttestData[,targetName])       #Scale target to [-1;1]


VMR <- BrownRobinsonReshetovAlgorithmTrain(x = trainX, y = trainY)

trainYpredicted <- predict(VMR, trainX)
fronttestYpredicted <- predict(VMR, fronttestX)

trainWinrate <- (sum(trainYpredicted>=0 & trainY>=0) + sum(trainYpredicted<0 & trainY<0)) / length(trainY)
fronttestWinrate <- (sum(fronttestYpredicted>=0 & fronttestY>=0) + sum(fronttestYpredicted<0 & fronttestY<0)) / length(fronttestY)

trainError <- mean(abs(trainY-trainYpredicted))
fronttestError <- mean(abs(fronttestY-fronttestYpredicted))

cat("Win rate on train sample:", trainWinrate, "\n")
cat("Average error on train sample:", trainError, "\n")
cat("Win rate on fronttest sample:", fronttestWinrate, "\n")
cat("Average error on fronttest sample:", fronttestError, "\n")
cat("Error < 1 is ok; >= 1 is a failure of the model\n")

endTime <- Sys.time()

endTime - startTime