In this project we will simulate the evolution of a colony of pipe spiders. Pipe spiders live by eating the flies that they catch flying down their pipes. An individual spider’s survival fitness depends in a large part on the geometry of its web, as this determines how effective it is in catching flies. We will suppose that our spiders live in pipes of radius 1. Webs are two dimensional, stretched across the pipe. That is, we can consider a web as existing in a circle of radius 1. A web will be described by a sequence of angles, which represent the points on the edge of the circle where the web is attached.
This project consists in programming functions about these pipe spiders and their webs and ability to catch flies. We will split our project in some tasks, that will be shown in the next sessions.



Task 1

Write a function webdraw that takes a web vector and draws it.

library(graphics)

webdraw = function(web) {
  
  plot(c(-1,1),c(-1,1),asp=1,type = "n",xlab="",ylab="")
  c1 = 0
  c2 = 0
  r = 1
  symbols(c1,c2,circles = r,inches=FALSE,add = TRUE)
  
  
  if(length(web)>1){
    for (i in 1:(length(web)-1)){
      segments(1*cos(web[i]),1*sin(web[i]),1*cos(web[i+1]),1*sin(web[i+1]))
    }
    
  }
}
##Testing the function above:

web <- c(0, 2*pi/3, 3*pi/2, pi/4)
webdraw(web)


Task 2

Write a function ldist that calculates the distances of a number of points from a given line. The function should take three arguments, a matrix p0 with two rows, a vector p1 of length 2, and a vector p2 of length 2. Each column of p0 represents a point (x0, y0), while p1 and p2 give the points (x1, y1) and (x2, y2) that define a line. The function should return a vector of length equal to the number of columns of p0, with the k-th entry equal to the distance of the k-th point (x0, y0) to the line given by (x1, y1) and (x2, y2). If p1 and p2 are identical then ldist should return the (Euclidean) distance between p0 and p1.

ldist = function (p0,p1,p2) {
  
  x0 = p0[1,]
  y0 = p0[2,]
  
  x1 = p1[1]
  y1 = p1[2]
  
  x2 = p2[1]
  y2 = p2[2]
  
  
  if((x1==x2) && (y1==y2)) { ##in case p1==p2 we should calculate the euclidian distance to one of these points
    answer = c()
    
    for (i in 1:ncol(p0)){
      answer[i] = sqrt(sum((p0[,i] - p1)^2))
    }
    return(answer)
    
  }else {
    return(abs((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1))/sqrt(((x2-x1)^2)+((y2-y1)^2)))
  }
}
##Testing the function above with the following input to check if it is working properly:

set.seed(90045)
tha <- runif(5, 0, 2*pi)
R <- sqrt(runif(5, 0, (1 - 0.05)^2))
p0 <- rbind(R*cos(tha), R*sin(tha))
ldist(p0, c(1, 0), c(0, 1))
## [1] 1.3085946 0.8867301 0.1969623 0.1138217 0.3329909
ldist(p0, c(0, 0), c(0, 0))
## [1] 0.6659744 0.8523597 0.9153808 0.9109671 0.4329158


Task 3

Write a function caught that takes as input a web vector and a fly matrix, and returns a logical vector indicating which flies have been caught. We suppose that a fly is caught if it touches any of the strands that make up the web. That is, if the centre of the fly is within distance 0.05 of any one of the strands.

caught = function (web,p0) {
  
  x = 1*cos(web)
  y = 1*sin(web)
  
  ans = rep(FALSE,ncol(p0))
  ind = c()
    for (j in 1:(length(web)-1)){
      aux = ldist(p0,c(x[j],y[j]),c(x[j+1],y[j+1])) <= 0.05
      for (i in 1:length(aux)) {
        if(aux[i] == TRUE){
          ind = c(ind,i)
        }
      }
    }
  ans[ind] = TRUE
  
  return (ans)
  
}
## To check that this function is working, test it as follows.

web <- c(0, 2*pi/3, 3*pi/2, pi/4)
webdraw(web)
cw <- caught(web, p0)
for (i in 1:length(cw)) {
  if (cw[i]) {
    symbols(p0[1,i], p0[2,i], circles = 0.05, inches = FALSE,
            add = TRUE, bg = "red")
  } else {
    symbols(p0[1,i], p0[2,i], circles = 0.05, inches = FALSE,
            add = TRUE)
  }
}

Task 4

The efficiency of a web is a trade off between its ability to catch flies and the length of thread used. For a given web w, let p be the probability that it catches any given fly, and a the length of the web, then we define the efficiency of the web to be:

100p − a


Write a function efficiency that takes as input a web vector, and returns an estimate of its efficiency. This requires an estimate of p and the calculation of a. For the former, generate an iid sample of 1000 fly centres, and use the proportion caught in the web to estimate p. The length a can be calculated exactly.

efficiency = function(web){

  
  theta = runif(1000,0,2*pi)
  R = sqrt(runif(1000,0,(1-0.05)^2))
   
  flyMatrix = rbind(R*cos(theta),R*sin(theta))
  
  aux = caught(web,flyMatrix)
  
  
  
  
  lengthWeb = 0
  for(i in 1:(length(web)-1)){ 
    lengthWeb =lengthWeb + sqrt(sum((c(1*cos(web[i]),1*sin(web[i])) - c(1*cos(web[i+1]),1*sin(web[i+1])))^2))
  }
  
  
  p = sum(as.numeric(aux))/ncol(flyMatrix)         
  a = lengthWeb
  
  efficiency = 100*p - a
  
  if (efficiency < 0) {              ##making sure that it never assumes negative values
    efficiency = 0
  }
  
  return(efficiency)
}
##To check that the code is working, we will test it as follows

efficiency(web)
## [1] 12.38834
n <- 100
effvec <- rep(0, n)
for (i in 1:100) effvec[i] <- efficiency(web)
r <- 2*sd(effvec)/sqrt(n)
mean(effvec) + c(-r, 0, r)
## [1] 11.89658 12.16334 12.43010


We would like to optimise webs for efficiency. To do this we will use an evolutionary algorithm. Evolutionary algorithms are stochastic optimisation algorithms, meaning they employ randomness in their search for an optimal solution. They can also cope with random errors in the function they are optimising, which is the situation we have since our efficiency function is a Monte Carlo estimate of the true efficiency. Most local search techniques have a current solution which is updated iteratively. An evolutionary algorithm maintains a population of solutions, which are updated iteratively. The n-th iteration is often referred to as the n-th generation. Evolutionary algorithms have three main components, which mimic the process of natural selection. Each generation is subject to

mutation: each individual in the population has a chance of changing randomly
recombination: randomly chosen pairs of individuals are combined to produce new individuals, which we add to the population
selection: the population is reduced to its original size by selecting the most efficient individuals


Task 5

Write a function mutate that takes as input a web vector and returns it altered according to the following rules:

• With probability 0.25 remove an element of the vector at random. If it has length 2 then do nothing.
• With probability 0.25 insert an element into the vector. The new point should be uniformly distributed on the interval [0, 2π), and be equally likely to appear at any point in the new vector.
• With probability 0.5 each element of the vector is perturbed, independently of the others, by an amount normally distributed with mean 0 and standard deviation 0.1.

mutate = function(web) {
  
  
  p1 = runif(1,0,1)
  
  
  # removing
  if(p1 <= 0.25) {
    if(length(web)==2){
      web = web  #break
    }else{
      remove = sample(1:length(web),1)
      web = web[-remove]
    }
  }
  
  #inserting
  if(p1 > 0.25 && p1 <= 0.5) {
    
    newPoint = runif(1,0,2*pi)
    while(newPoint==(2*pi)) {         ##to make sure that newPoint assumes values in the interval [0,2pi)
      newPoint = runif(1,0,2*pi)
    }
    
    ind = sample((1:(length(web)+1)),1)
    
    web <- append(web, newPoint , after=(ind-1))
  }
  
  #pertubating
  if (p1 > 0.5) {
    for(i in 1:length(web)) {
      web[i] = web[i] + rnorm(1,0,0.1)
    }
  }
  
  return(web)
  
}


Task 6

Write a function recombine that takes as input two web vectors, a and b say, and returns a list of two web vectors, obtained by recombining the input vectors. The recombination should work as follows:

• Generate p ∈ (0, 1) uniformly.
• Split a into two parts, a1 and a2, with lengths roughly in the proportion to p : 1 − p, ensuring that both a1 and a2 have length at least 1. Do the same for b (using the same p).
• The two new vectors are obtained by combining a1 and b2, and combining b1 and a2

Note that the recombination step makes it possible to get repeated values in web vectors. It is for this reason that ldist needs to be able to deal with identical p1 and p2.
When you have the evolutionary algorithm working, use it to obtain a web with efficiency of at least 40.

recombine = function(webA,webB) {
  
  
  p = runif(1,0,1)
  

  
  
  splitA = round(length(webA)*p)
  a1 = sample(webA,splitA,replace=FALSE)
  a2 = setdiff(webA,a1)
  
  #making sure that none of these vectors have 0 elements:
  if(splitA == 0) {
    a1 = sample(webA,1,replace=FALSE)
    a2 = setdiff(webA,a1)
    
  }
  if(splitA == length(webA)){
    a1 = sample(webA,(length(webA) - 1),replace=FALSE)
    a2 = setdiff(webA,a1)
    
  }

 
  
  splitB = round(length(webB)*p)
  b1 = sample(webB,splitB,replace=FALSE)
  b2 = setdiff(webB,b1)
  
  #making sure that none of these vectors have 0 elements:
  if(splitB == 0) {
    b1 = sample(webB,1,replace=FALSE)
    b2 = setdiff(webB,b1)
    
  }
  if(splitB == length(webB)){
    b1 = sample(webB,(length(webB) - 1),replace=FALSE)
    b2 = setdiff(webB,b1)
    
  }
  
  
  
  newWebA = c(a1,b2)
  newWebB = c(b1,a2)
  
  
  newvectors = list(newWebA,newWebB)
  
  return(newvectors)
}
###Using the evolutionary algorithm below to obtain a web with efficiency of at least 40###


popsize <- 16 # size of population, must be even
r_mutate <- 0.25 # proportion of each generation mutated
set.seed(2064) #setting one seed that will create a web with the desired efficiency.

# set up intial population of webs, each consisting of a single strand:
population <- list()
for (i in 1:popsize) {
  population <- c(population, list(runif(2, 0, 2*pi)))
}

# main loop:
n_gen <- 100 # num of generations that will be generated
total_gen <- 0 # total generations so far


for (gen in 1:n_gen) {
  # mutate
  for (i in 1:popsize) {
    if (runif(1) < r_mutate) {
      population[[i]] <- mutate(population[[i]])
    }
  }
  # recombine
  x <- sample(1:popsize, popsize/2)
  y <- (1:popsize)[-x]
  for (i in 1:(popsize/2)) {
    population <- c(population, recombine(population[[x[i]]],
                                          population[[y[i]]]))
    
  }
  # select
  fitness <- rep(0, 2*popsize)
  for (i in 1:(2*popsize)) {
    fitness[i] <- efficiency(population[[i]])
  }
  fittest <- order(fitness, decreasing=TRUE)
  population <- population[fittest[1:popsize]]
}

# report:
total_gen <- total_gen + n_gen
cat("best efficiency after", total_gen,
    "generations is", fitness[fittest[1]], "\n")
## best efficiency after 100 generations is 40.75817
webdraw(population[[fittest[1]]]) #the web that had the best efficiency (and higher than 40) in this simulation process


Task 7

See if you can design a web with efficiency greater than 45.

After many trials and simulation using different algorithms configuration and different web sizes, the web with the best efficiency that we have found is shown below:

aux = seq(0,2,by=1/16)
aux = aux[-33]
webx = aux*pi



set.seed(12787)


popsize <- 28# size of population, must be even

# set up intial population of webs, each consisting of a single strand:
population <- list()
for (i in 1:popsize) {
  population <- c(population, list(sample(webx)))
}

# main loop:
n_gen <- 45 # num of generations that will be generated
total_gen <- 0 # total generations so far


for (gen in 1:n_gen) {
  
  # recombine
  x <- sample(1:popsize, popsize/2)
  y <- (1:popsize)[-x]
  for (i in 1:(popsize/2)) {
    population <- c(population, recombine(population[[x[i]]],
                                          population[[y[i]]]))
    
  }
  # select
  fitness <- rep(0, 2*popsize)
  for (i in 1:(2*popsize)) {
    fitness[i] <- efficiency(population[[i]])
  }
  fittest <- order(fitness, decreasing=TRUE)
  population <- population[fittest[1:popsize]]
}

# report:
total_gen <- total_gen + n_gen
cat("best efficiency after", total_gen,
    "generations is", fitness[fittest[1]], "\n")
## best efficiency after 45 generations is 43.2209
webdraw(population[[fittest[1]]]) #the web that had the best efficiency of this simulation process



Conclusion


So, with all these 7 tasks, this project is done. I’ve learned a lot about programming, logic, statistics and simulation tools during this project. And it felt very good to practice using all of this knowledge.