# We depend on GenSA, as we want to run one of the algorithms from the paper.

if (!require(GenSA)) install.packages("GenSA")
library(GenSA)

# This is a port of fitness.m / initialize.m, taken from http://staff.itee.uq.edu.au/marcusg/msg.html
# Note that both files are a part of this function:
# - initialize.m comes before the function to be returned,
# - fitness.m comes in function(x) { ... }.

# Please also note that fitness is not signature-compatible with the Matlab version,
# e.g. it accepts exactly one individual unlike Matlab which takes a vector of individuals,
# the lower-upper arguments come in this order, and ratio control is more refined.

# Of course we don't use global variables anymore.

GaussianLandscapeGenerator <- function(dimension, nGaussian, lower, upper, globalValue, ratioMin, ratioMax, optVariance) {
    if (dimension <= 1) {
        stop(paste("Illegal value for argument 'dimension': shall be > 1, found", dimension))
    }
    if (nGaussian <= 0) {
        stop(paste("Illegal value for argument 'nGaussian': shall be > 0, found", nGaussian))
    }
    if ((length(upper) != dimension) || (length(lower) != dimension)) {
        stop(paste("'dimension' and dimensions of 'lower' and 'upper' do not match:", dimension, length(lower), length(upper)))
    }
    if (sum(upper <= lower) > 0) {
        lower_txt <- paste(lower, collapse=" ")
        upper_txt <- paste(upper, collapse=" ")
        stop(paste("Illegal value for arguments 'upper' and 'lower': they shall be componentwise '>', found", upper_txt, "and", lower_txt))
    }
    if (globalValue <= 0) {
        stop(paste("Illegal value for argument 'globalValue': shall be > 0, found", globalValue))
    }
    if (ratioMin <= 0 || ratioMin > ratioMax || ratioMax >= 1) {
        stop(paste("Illegal value for arguments 'ratioMin', 'ratioMax': shall be 0 < ratioMin <= ratioMax < 1, found", ratioMin, ratioMax))
    }

    rotation <- array(dim=c(dimension,dimension,nGaussian))
    for (i in 1:nGaussian) {
        vi <- diag(nrow = dimension)
        for (j in 1:(dimension-1)) {
            for (k in (j+1):dimension) {
                r <- diag(nrow = dimension)
                alpha <- runif(1) * pi / 2 - pi / 4
                r[j, j] <- cos(alpha)
                r[j, k] <- sin(alpha)
                r[k, j] <- -sin(alpha)
                r[k, k] <- cos(alpha)
                vi <- vi %*% r
            }
        }
        rotation[,,i] <- vi
    }

    varianceRange <- (upper - lower) / 20
    variance <- matrix(runif(nGaussian * dimension), nrow = nGaussian, ncol = dimension) * varianceRange + 0.05
    meanVector <- matrix(runif(nGaussian * dimension), nrow = nGaussian, ncol = dimension) * (upper - lower) + lower

    optimumValue <- runif(nGaussian) * globalValue * (ratioMax - ratioMin) + ratioMin

    variance[1,] <- optVariance
    optimumValue[1] <- globalValue
    meanVector[1,] <- meanVector[2,] + 1e-3

    covMatrixInv <- array(dim=c(dimension,dimension,nGaussian))
    for (i in 1:nGaussian) {
        covMatrix <- diag(variance[i,])
        covMatrixInv[,,i] <- solve(t(rotation[,,i]) %*% covMatrix %*% rotation[,,i])
    }

    list(fun = function(x) {
        tmp <- array(0, dim = nGaussian)
        for (i in 1:nGaussian) {
            newx <- x - meanVector[i,]
            z <- (newx %*% covMatrixInv[,,i]) * newx
            tmp[i] <- t(sum(z))
        }
        -max(exp(-0.5 * tmp / dimension) * optimumValue)
    }, optimum = meanVector[1,])
}

dimension <- 2

lower <- sapply(1:dimension, function(i) {0})
upper <- sapply(1:dimension, function(i) {1})

runOne <- function(test, name) {
    function(i) {
        result <- GenSA(fn = test, lower = lower, upper = upper)
        message(paste("  ", name, "run", i, ":", result$value))
        result$value
    }
}

measurements <- 1:6

for (attempt in 1:10) {
    test <- GaussianLandscapeGenerator(dimension, 10, lower, upper, 1, 0.5, 0.5 + 1e-9, 1e-7)
    message(paste("Attempt", attempt))
    lhsI <- sapply(measurements, function(i) {1})
    rhsI <- sapply(measurements, function(i) {2})
    lhsV <- sapply(1:6, runOne(test$fun, "GenSA v1"))
    rhsV <- sapply(1:6, runOne(test$fun, "GenSA v2"))
    df <- data.frame(algo <- c(lhsI, rhsI), value <- c(lhsV, rhsV))
    print(summary(aov(value ~ algo, data = df)))
    message("")
}