##################################################################################################################
# Aluno: Tiago Oliveira Beraldi #
# Numero USP: 4896352 #
# Exerci苞io Programa 6: Programa誽o Gen彋ica #
# MAP2212 Laboratorio de Computacao e Simulacao BMAC 2009 #
##################################################################################################################
PG <- function(arquivo) {
# leitura e processamento do arquivo de input
input <- scan(arquivo)
alfa <- input[1]
beta <- input[2]
gama <- input[3]
M <- input[4]
N <- input[5]
b <- input[6]
mb <- floor(M/b) # numero ideal de linhas de cada cor
Smax <- 2*(M+N) # quantidade de iteracoes do simulated annealing
A <- array(0, dim = c(M, N)) # matriz a do problema
k <- floor(runif(1, 1, b + 2)) # cor aleatoria
pop <- min(b, 500) # tamanho da populacao
P <- array(k[1], dim = c(M, pop)) # matriz com populacao
W <- array(0, dim = c(b + 1, N, pop)) # numero de atividades da coluna j na cor k=1..b+1 para a populacao
# inicializa e atualizacao das matrizes A, W
i <- 7
while (i < length(input)) {
A[input[i], input[i+1]] <- 1
W[k[1], input[i+1],] <- W[k[1], input[i+1],] + 1
i <- i + 2
}
# apaga arquivo saida.txt
unlink("saida.txt")
# guarda melhor e pior solucao
sol <- array(0, dim=c(2, 2))
sol[1,] <- c(Inf, 0) # melhor solucao
sol[2,] <- c(0, 0) # pior solucao
nBrk <- 0 # Verifica se algoritmo PG esta congelado
# inicio do processamento
while (pop > 1) {
# simulated annealing
for(i in 1:pop) {
ret <- HSA(alfa, beta, gama, A, P[,i], W[,,i], M, N, b, mb, Smax)
nBrk <- nBrk + 1
# atualiza individuo
P[,i] <- ret$X[,1]
W[,,i] <- ret$W
# atualiza melhor solucao
if (sum(ret$fx[1:3]) < sol[1,1]) {
sol[1,] <- c(sum(ret$fx[1:3]), i)
nBrk <- 0
}
# atualiza pior solucao
if (sum(ret$fx[1:3]) > sol[2,1])
sol[2,] <- c(sum(ret$fx[1:3]), i)
}
# gera arquivo de saidas
saida <- array()
saida[1] <- sol[1,1]
saida[2] <- length(P[,sol[1,2]])
saida[3] <- P[1, sol[1,2]]
for (i in 2:length(P[,sol[1,2]])) {
saida[3] <- paste(saida[3], P[i,sol[1,2]], sep = " ") }
# grava melhor solucao parcial
write(saida, file="saida.txt", sep = "\n", append=TRUE)
# finaliza algoritmo se PG congelou
if (nBrk > b) break
# programacao genetica
# calculo da probabilidade de metropolis para aceitacao do nascimento/mutacao ou morte do individuo
metrop <- exp(-(sol[2,1] - sol[1,1]))
# nascimento/mutacao
if (runif(1) > metrop) {
# crossing-over
for (i in 1:M) {
# coloca parte da melhor solucao na pior
if ((i %% 2) == 0)
P[i, sol[2,2]] <- P[i, sol[1,2]]
}
# reinicializa pior solucao
sol[2,] <- c(0, 0)
}
# mata pior individuo
else {
pop <- length(P[1,])-1
novoX <- array(0, dim=c(M, pop))
novoX[, 1:sol[2,2]] <- P[, 1:sol[2,2]]
novoX[, sol[2,2]:pop] <- P[, (solucao[2,2]+1):(pop + 1)]
X <- novoX
}
}
}
##################################################################################################################
# Simulated Annealing #
##################################################################################################################
HSA <- function(alfa, beta, gama, A, X1, W1, M, N, b, mb, Smax) {
W <- array(0, dim = c(b + 1, N)) # numero de atividades da coluna j na cor k=1..b+1
X <- array(0, dim = c(M, 1)) # numero da cor k=1..b+1 da linha i=1..M
W <- W1
X[,1] <- X1
# inicializacao de fx
fx <- vector()
aux <- as.vector(table(X, exclude = b + 1))
fx[1] <- alfa * (sum((aux - mb)^2) + ((b - length(aux)) * mb^2))
fx[2] <- 0
fx[3] <- gama * sum(X[, 1] == b + 1)
fx[4] <- 0
# parametros do simulated annealing
theta <- b
mi <- 1 / sqrt(theta)
e1 <- 1 / (b + 2)
e2 <- 1 / (b + 1)
stop <- M * N
step <- 1
# iteracoes
for(i in 1:Smax) {
k <- floor(runif(1, 1, b + 2)) # cor aleatoria
linha <- floor(runif(1, 1, M + 1)) # linha aleatoria
corold <- X[linha, 1]
cornew <- k
if (corold != cornew) {
fx1 <- CalculaFX(A, X, W, alfa, beta, gama, b, mb, corold, cornew, linha, theta, mi)
# calculo da probabilidade de metropolis
delta <- (sum(fx1) - sum(fx))
if (delta <= 0) {
metrop <- 1
}
else {
metrop <- exp(-(theta * delta))
}
# aceita novo fx com probabilidade de metropolis e atualiza fx, W, X
if (metrop > runif(1)) {
W[corold, ] <- W[corold, ] - A[linha, ]
W[cornew, ] <- W[cornew, ] + A[linha, ]
X[linha, 1] <- cornew
step <- i
fx <- fx1
}
}
# criterio de resfriamento
if ((i %% stop) == 0) {
theta <- theta * (1 + e1)
mi <- mi * (1 + e2)
}
# criterio de parada
if ((i - step) > (M * stop)) break
}
list(X=X, W=W, fx=fx)
}
##################################################################################################################
# Calculo da funcao de custo #
##################################################################################################################
CalculaFX <- function(A, X, W, alfa, beta, gama, b, mb, corold, cornew, linha, theta, mi) {
# atualiza valores de X, W
W[corold, ] <- W[corold, ] - A[linha, ]
W[cornew, ] <- W[cornew, ] + A[linha, ]
X[linha, 1] <- cornew
fx <- vector()
# atualiza hx2
aux <- as.vector(table(X, exclude = b + 1))
fx[1] <- alfa * (sum((aux - mb)^2) + ((b - length(aux)) * mb^2))
# atualiza cx
I <- array(1, dim=c(1,b)) %*% (W[1:b,] > 0)
fx[2] <- beta * sum(I > 1)
# atualiza rx
fx[3] <- gama * sum(X[, 1] == b + 1)
# atualiza ux
fx[4] <- sum(ifelse(I > 1, I, 0))/mi
return (fx)
}
PG("entrada.txt")