example

example.r

# install.packages("shape") # in case shape is not installed
require(shape)
# install.packages("plotrix") # in case plotrix is not installed
require(plotrix)

##################################################################################################################################
##################################################################################################################################
############################################# Helper functions   #################################################################
##################################################################################################################################
##################################################################################################################################
plottitle <- function(t){
  ### Plot "Folie ###"
  plot(NA,NA,frame.plot=FALSE, xlab="", ylab="", xlim= c(-1,1), ylim=c(-1,1), axes=FALSE)
  text(0,0,t, cex=2, )
}


getPolar <- function(n0,n1){
  # Polarwinkel berechnung zum vekürzen der Pfeile     
  x1 <- n1$x
  x0 <- n0$x
  y1 <- n1$y
  y0 <- n0$y
  value<-NULL
    
  x<-x1-x0
  y<-y1-y0
  if(x==0){
    if(y>0){
      value<-pi/2
    }else{
      if(y<0){
        value<- 3*pi/2
      }else{
        # warning("Angle can't be calculated because both coordinates are the same. Returning NA.")
        return(NA)
      }
    }
  }else{
    value<-atan(y/x)
    if(x<0){
      value<-value+pi
    } 
  }  
  return(value)
} 

getNewPoint <- function(p1, p2, short){
  # verkürzen der Pfeile damit es "hübscher aussieht"
  m <- getPolar(p1,p2)
  if(is.na(m)){
    return(NA)
  }else{
    r <- sqrt((p1$x - p2$x )^2 + (p1$y - p2$y)^2)
    
    x <- p1$x + (r - short) * cos(m)  
    y <- p1$y + (r - short) * sin(m) 
    
    return (data.frame(x=x, y=y))
  }
}
getNumberOfStops <- function(tour){
  # wieviele Stops führen zum Depot zurück?
  n <- which(tour==1) # 1=depot 
  return (length(n))
}
calcTourplan <- function(cij, tour){
  
  F <- 0
  
  # Speichere die Orte, die keinen Vorgänger haben --> Also vom Depot aus angesteuert werden
  hasPreDecessor <- rep(FALSE, nrow(locations)) 
  for(i in 1:ncol(tour)){
    
    F<- F + cij[i+1, tour[i]]
    hasPreDecessor[tour[i]] <- TRUE 
    
    
  }
  
  I <- which(hasPreDecessor==FALSE) 
  # Welche knoten sind mit dem Depot verbunden?
  
  if(length(I)>0){ 
    for(i in I){ 
      F<- F + cij[1, i] # i = 1 is depot
    }
  }
  return (F)
}

printTour <- function(locations, tour, t, short= 3.5){
  # Erzeuge einen Plot mit dem Depot  
  plot(locations$x[1], locations$y[1], pch=15, cex=4, col=4, xlim=range(locations$x)*1.1, ylim=range(locations$y)*1.1, xlab="", ylab="", asp=1)
  title(t) #füge den Titel hinzu
  
  # finde die Orte, die keinen Vorgänger haben --> Also vom Depot aus angesteuert werden
  hasPreDecessor <- rep(FALSE, nrow(locations)) 
  for(i in 1:ncol(tour)){
  
    p1 <- locations[i+1,] # i = 1 is depot
    p2 <- locations[tour[i], ] 
    p<-getNewPoint(p1,p2,short) 
    
    if(length(p)>1){
      if(tour[i]>1){
        hasPreDecessor[tour[i]] <- TRUE # wird nicht mit dem Depot verbunden
      }
      #verbinde mit Nachfolger
      arrows( p1$x, p1$y, p$x, p$y, length =.2, lwd= 1 , col=2 , angle = 20 )    
    }
  }
   
  I <- which(hasPreDecessor==FALSE)
  # verbinde Startpunkte von Touren mit dem Depot (werden blau eingezeichnet)
  if(length(I)>0){
    p1<-locations[1,] # i = 1 is depot
    for(i in I){ 
      p2 <- locations[i,] # i = 1 is depot
      p<-getNewPoint(p1,p2,short)
      if(length(p)>1){       
        arrows( p1$x, p1$y, p$x, p$y, length =.2, lwd= 1, col=4 , angle = 20 ,lty=1 )    
      }
    }
  }
  #Zeichne die Knoten ein
  points(locations$x,locations$y, pch=21, cex=3, bg="white")
  #und schreibe den Namen in die Knoten
  text(locations$x, locations$y, locations$id, cex=0.8)
  # ALternative labeling methode
  # spread.labels(locations$x, locations$y, locations$id, 0.2 , cex=0.8)
}

##################################################################################################################################
##################################################################################################################################
############################################# Genetischer Algorithmus  ###########################################################
##################################################################################################################################
##################################################################################################################################


recombine<- function(tours, f, m, recomb){
  # wie auf Folie 174-177 beschrieben
  # tours = Genpool
  # f = father index
  # m = mother index
  # recomb = Predefined Recombination vector
  
  t.f <- tours[f,] #get tourplan father
  t.m <- tours[m,] #get tourplan mother
  
  t.c <- t.f # use father as template
  
  # Recombing
  for(i in 1:ncol(recomb)){
    if(recomb[i] == 1) {
      # wenn mutter chrom verwendet werden soll,
      # ersetze dies
      t.c[i] <- t.m[i]
    }
  }
  # prepare Return
  t.c <- matrix(t.c, nrow=1)
  colnames(t.c)  <- colnames(tours)
  return (t.c)
}


repair1 <- function(tour){
  # wie auf Folie 178 beschrieben
  #cond 1/2 once arrive/leave
  
  #substitute multiple indices by 0(1) --> ALL!
  #getunique stops
  y<-unique(tour) 
  y<- y[y!=1]  # 1 --> 0 # elminate Depot-Stops
  
  #iterate over potentail multiple Stop indices
  for(i in y){
    # which stops have this index?
    x <- which(tour == i)
    
    if(length(x)>1){ 
      # if i is multiple stop       
      tour[x] <- 1 # return to depot
    }
  }
  # Return 
  return (tour)
}


repair2 <- function(tour){
  
  #repair short cycles (DFJc)
  
  # iteriere über jeden Stop im Tourenplan 
  for(i in 1:length(tour)){
    # wenn der Nachfolger nicht das Depot ist (hier index = 1)
    # dann kann es Kurzzyklen geben
    if(tour[i]>1){ # not returning to depot
      
      # finde sukzezzive jeden nachfolger in der Tour --> Rekursiv
      r<-searchFor(tour, tour[i], i+1) #recursive call
      
      # Ein kurzzyklus liegt vor, wenn der Startknoten wieder erreicht wird.  
      # in diesem Fall ist r = FALSE
      if(r == FALSE){ #if shortcycle found, return current node to depot
        tour[i] <- 1
      }
    }
  }
  #elimiate cycles
  return (tour)
}

searchFor <- function(tour, index, start){
  #helper function für recursiven Aufruf
  
  if(tour[index -1] == 1) # 1 == depot
  {
    return(TRUE)
  }else{
    if(tour[index -1] == start) # Kurzzyklus gefunden
    {
      return(FALSE)
    }else{
      return(searchFor(tour, tour[index -1], start)) # erneuter Rekursiver aufruf
    }
  }
}

# Capacity Check
repair3 <- function(tour, locations=NULL, capacity = NULL){
  
  #not implemented yet
  return (tour)
}

# mutation
mutate <- function(tour){
  # vgl. Folie 179
  # suche zufällig einen Start knoten aus, der mit einem Endknoten verbunden wird.
  
  # combine which stop o[1] with wich stop o[2]
  o <- sample(1:length(tour), 2)
  tour[o[1]] <- o[2] + 1 # 1 is depot 
  
  # mutated tour 
  #immediate Repair
  tour <- repair1(tour)
  tour <- repair2(tour)
  tour <- repair3(tour)
  
  return (tour) # mutated + repaired tour
}

##################################################################################################################################
##################################################################################################################################
############################################# Skript starts here #################################################################
##################################################################################################################################
##################################################################################################################################

set.seed(1) # make sure to get same random results
n<-9
p <- data.frame(id= paste("P",1:n, sep=""),  x = round( runif(n,-100,100)), y= round(runif(n,-100,100)))
p$x[1] <- 0
p$y[1] <- 0

p$x[2:n] <- c(-40, -10, 10, 30, 50,  40,  10, -40)
p$y[2:n] <- c( 40,  50, 40, 50, 20, -20, -30, -20)

p$id <- c("L0", paste("K",1:(n-1), sep=""))
rownames(p) <- p$id

cij <- matrix(rep(0, n^2), nrow=n)
rownames(cij) <- p$id
colnames(cij) <- p$id


for(i in 1:(n-1)){  
  for(j in (i+1):n){ 
    cij[i,j] <- round(sqrt((p$x[i]-p$x[j])^2 + (p$y[i]-p$y[j])^2))
    cij[j,i] <- cij[i,j]      
  }  
}
recomb <- matrix(c(0,0,0,1,1,0,1,1), nrow=1)
colnames(recomb)  <- rownames(p)[-1]

T<-4

tours<- matrix(rep(NA, T*(n-1)), ncol=n-1)
colnames(recomb)  <- rownames(p)[-1]
tours[1,] <- c(2,3,4,0,0,7,0,0)+1
tours[2,] <- c(0,3,5,6,0,0,8,0)+1
tours[3,] <- c(2,0,0,3,4,0,8,0)+1
tours[4,] <- c(0,0,2,3,4,0,8,0)+1  
 
l<- layout(matrix(c(1,1,2:5), byrow=TRUE, ncol=2),heights=c(1,10,10)) # adjust plot
par(mar = rep(0.1,4)) # c(bottom, left, top, right) # adjust margin
plottitle("Folie 175 - initial set of solutions (gene-pool)")
par(mar = c(2.5,2.5,2.5,0.5)) # c(bottom, left, top, right) # adjust margin
## Folie 175
for(i in 1:T){
  t <- paste("Parent Tour", i) 
  printTour(p, t(as.matrix(tours[i,], nrow=1)), t)     
}
cat ("Press [enter] to continue")
line <- readline()
##################### Stop here for first image ##########################

# combine parent 1 + 4 --> vgl. Folie 177
p1<-1
p2<-4
t.c <- recombine(tours, p1,p2, recomb) 
# add to genpool
tours <- rbind(tours, t.c) 

#choose Tourenplans for next image
T <- c(p1, p2, nrow(tours))
# set Layout
l<- layout(matrix(c(1,1,2,3,4,4), byrow=TRUE, ncol=2),heights=c(1,10,10)) # adjust plot
par(mar = rep(0,4)) # c(bottom, left, top, right) # adjust margin

##### new image ###
plottitle("Folie 177 -- re-combine (parent x parent >child)")
par(mar = c(2.5,2.5,2.5,0.5)) # c(bottom, left, top, right) # adjust margin

for(i in T){ 
  # adjust title
  if(i ==  nrow(tours)){ 
    t <- paste("Child Tour: [", paste(tours[i,]-1, sep="",  collapse=","), "]\n","Recombined with: [", paste(recomb, sep="",  collapse=","), "]", sep="")
  }else{
    t <- paste("Parent Tour ", i,": [", paste(tours[i,]-1, sep="",  collapse=","), "]", sep="")
  } 
  printTour(p, t(as.matrix(tours[i,], nrow=1)), t)  
} 
cat ("Press [enter] to continue")
line <- readline()
##################### Stop here for second image ##########################

l<- layout(matrix(c(1,1,2:5), byrow=TRUE, ncol=2),heights=c(1,10,10)) # adjust plot
par(mar = rep(0,4)) # c(bottom, left, top, right) # adjust margin
plottitle("Folie 178 -- repair chromosomes")
par(mar = c(2.5,2.5,2.5,0.5)) # c(bottom, left, top, right) # adjust margin

# plot child tour
t <- paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "]\n","Recombined with: [", paste(recomb, sep="",  collapse=","), "]", sep="")
printTour(p, t.c, t)

#Repair Step 1
t.c <- repair1(t.c) 
t <-paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "] \nRepair 1: Substitue by 0", sep="")
printTour(p, t.c, t)

#Repair Step 2
t.c <- repair2(t.c) 
t <-paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "] \nRepair 2: elim. short cycles", sep="")
printTour(p, t.c, t)

#Repair Step 3
t.c <- repair3(t.c) 
t <-paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "] \nRepair 3: Capacity Check (ommited)", sep="")
printTour(p, t.c, t)
cat ("Press [enter] to continue")
line <- readline()
##################### Stop here for Thrid image ##########################

l<- layout(matrix(c(1,1,2:3), byrow=TRUE, ncol=2),heights=c(1,5)) # adjust plot
par(mar = rep(0,4)) # c(bottom, left, top, right) # adjust margin
plottitle("Folie 179 -- mutation")
par(mar = c(2.5,2.5,2.5,0.5)) # c(bottom, left, top, right) # adjust margin

t <-paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "] \n3-Step-Repair processed", sep="")
printTour(p, t.c, t)
 
t.c <- mutate(t.c) 
t <-paste("Child Tour: [", paste(t.c-1, sep="",  collapse=","), "] \nmutated", sep="")
printTour(p, t.c, t)

tours[nrow(tours),] <- t.c # update tourplan in genpool
cat ("Press [enter] to continue")
line <- readline()
##################### Stop here for fourth image ##########################


l<- layout(matrix(c(1:6), byrow=TRUE, ncol=2)) # adjust plot
par(mar = rep(0,4)) # c(bottom, left, top, right) # adjust margin
plottitle("Folie 180 \nCalculate Fitness\nof genpool")
par(mar = c(2.5,2.5,2.5,0.5)) # c(bottom, left, top, right) # adjust margin

for( i in 1: nrow(tours)){
  tp <- t(as.matrix(tours[i,], nrow=1))
  costs <- calcTourplan(cij, tp)
  n <- getNumberOfStops(tp)
  t <-paste("Tour ",i," [", paste(tp-1, sep="",  collapse=","), "] \nlength: ",round(costs)," -- # Tours: ", n, sep="")
  printTour(p, tp, t)
}

##################### Stop here for fifth image ##########################