GSA.R

initialization <- function(dims, N, low, up) {
    data <- runif(N * dims, low, up)
	matrix(data, nrow=N)
}


zeros <- function(N, dims) {
	matrix(rep(0, dims * N), nrow=N)
}


ones <- function(N, dims) {
	matrix(rep(1, dims * N), nrow=N)
}


space.bound <- function(X, low, up) {
    size = dim(X)
    N = size[1]
    dims = size[2]
    for(i in 1:N) {
        if (!all(X[i,] >= low & X[i,] <= up)) X[i,] <- runif(dims, low, up)
    }
    return(X)
}


move <- function(X, a, V) {
    size = dim(X)
    N = size[1]
    dims = size[2]
    V = matrix(data=runif(dims*N), nrow=N, ncol=dims) * (V + a)
    X = X + V
    return(list(X=X, V=V))
}


mass.calculation <- function(fit, minflag) {

    N = length(fit)

    Fmax = max(fit)
    Fmin = min(fit)

    if (Fmax == Fmin) {
        M = ones(N, 1)
    } else {
                                        # minimization
        if (minflag) {
            best = Fmin
            worst = Fmax
        } else {
                                        # maximization
            best = Fmax
            worst = Fmin
        }

        M = (fit - worst) / (best - worst)
    }

    M / sum(M)
}


Gconstant <- function(iteration, maxit) {
    alpha = 20
    G0 = 100
    G0 * exp(-alpha * iteration/maxit)
}


Gfield <- function(M, X, G, Rnorm, Rpower, ElitistCheck, iteration, maxit) {
    size = dim(X)
    N = size[1]
    dims = size[2]

    finalper = 2 # In Xthe last iteration, only 2 percent of agents apply force to the others.

                                        # Total force calculation
    if (ElitistCheck) {
        kbest = finalper + (1 - iteration/maxit) * (100 - finalper)
        kbest = ceiling(N * kbest/100)
    } else {
        kbest = N
    }

    Ms = sort(M, decreasing=TRUE, index.return=TRUE)

    E = matrix(nrow=N, ncol=dims)

    for (i in 1:N) {
        E[i,] = zeros(1,dims)
        for (ii in 1:kbest) {
            j = Ms$ix[ii]
            if (j != i) {
                R = dist(rbind(X[i,],X[j,]))
                for (k in 1:dims) {
                    E[i,k] = E[i,k] + runif(1) * M[j] * (X[j,][k] - X[i,][k]) / (R^Rpower + 0.001)
                }
            }
        }
    }

                                        # Acceleration
    a = E * G
    return(a)
}


GSA <- function(f, N, dims, maxit, low, up, ElitistCheck, minflag=TRUE, Rpower) {

    Rnorm = 2

    X = initialization(dims, N, low, up)
    V = zeros(N, dims)

    for (iteration in 1:maxit) {
        X <- space.bound(X, low, up)

        fitness = apply(X, 1, f)

        if (minflag) {
                                        # minimization
            best = min(fitness)

        } else {
                                        # maximization
            best = max(fitness)
        }

        bestX = which(fitness == best)

        if (iteration == 1) {
            Fbest = best
            Lbest = X[bestX,]
        }
        if (minflag == 1) {
                                        # minimization
            if (best < Fbest) {
                Fbest = best;
                Lbest = X[bestX,]
            }
                                        # maximization
            if (best > Fbest) {
                Fbest = best;
                Lbest = X[bestX,]
            }
        }

                                        # Calculation of masses
        M = mass.calculation(fitness, minflag)

                                        # Calculation of Gravitational constant
        G = Gconstant(iteration, maxit)

                                        # Calculation of acceleration in gravitational field
        a = Gfield(M, X, G, Rnorm, Rpower, ElitistCheck, iteration, maxit)

        res = move(X, a, V)

        X = res$X
        V = res$V

    }

    return(list(X=Lbest, Fbest=Fbest))

}