jb_pubdistance.r

jb_pubdistance <- function(
  v_pubs=c("The Cambridge Brew House","The Cambridge Blue",
           "King Street Run P.H.","The Elm Tree","The Regal"),
  v_location="Cambridge, UK",
  v_zoom=13,
  listpubs=FALSE,
  cam_pubs=TRUE,
  crow_distances=FALSE,
  units="minutes"){
  
  ## Read in pub data
  temporaryFile <- tempfile()
  download.file("https://gist.githubusercontent.com/epijim/47da5037b6bb0d750b9b/raw/1f5608b826d039384775f2de87fcab9f5c8311d9/cam_pubs.csv",
                destfile=temporaryFile, method="curl")
  pubs <- read.csv(temporaryFile)
  
  
  if(listpubs==TRUE){
    print(pubs$name)
    stop("listpubs TRUE: function stopped and pubs listed")
  }
  
  # required package
  library(ggmap)
  library(hash)
  # load map
  basemap <- get_map(location=v_location, 
                     zoom=v_zoom, source='google', 
                     maptype="terrain", color="bw") 
  
  
  if(cam_pubs==TRUE){
    if(crow_distances==FALSE){
      ## Read in distance data
      temporaryFile <- tempfile()
      download.file("https://gist.githubusercontent.com/epijim/1b049708cef27d9ddc34/raw/02a4f8559541b7dd5589ef1f530966759d81a711/campubdistances.csv",
                    destfile=temporaryFile, method="curl")
      googled_distances <- read.csv(temporaryFile)
      distance_hash <- hash(keys=paste(googled_distances$latlon_1,googled_distances$latlon_2),
                            values=googled_distances$minutes)
      
      pubs$latlon <- paste0(pubs$lat,",",pubs$lon)
      
      ## Choose pubs you want #########################################
      selected.pubs <- subset(pubs, name %in% v_pubs)
      
      ## Variables to set #############################################
      n <- as.numeric(nrow(selected.pubs)) # number of pubs
      M <- 5000 # number of iterations
      temp <- 100 # initial temperature
      finaltemp <- 0.1
      tempfactor <- (temp/finaltemp)^(-1/M)
      fixed.itr <- 20 # number of iteration for a fixed temperature 
      
      ## FUNCTIONS ####################################################
      
      ## Check if integer 
      check.integer <- function(x) {
        x == round(x)
      }
      
      ## Difference between 2 points
      jb_difference <- function(testv1,testv2){
        d <- distance_hash[[paste(testv1,testv2)]]
        if(is.null(d))
        {
          d <- distance_hash[[paste(testv2,testv1)]]
        }
        d <- as.numeric(d)
        return(d)
      }
      
      distance <- function(ord){
        p <- length(ord)
        d <- rep(0,p)
        for(i in 1:(p-1)){
          d[i] <- jb_difference(selected.pubs$latlon[ord[i]],
                                selected.pubs$latlon[ord[i+1]])
        }
        d[p] <- jb_difference(selected.pubs$latlon[ord[1]],
                              selected.pubs$latlon[ord[p]] )
        distance=sum(d)
        return(distance) 
      }
      
      ## Initial values ################################################
      
      cost1 <- rep(0,M)
      numaccept <- 0 # to keep track of accepting
      templist <- rep(0,M) # to keep track of temperature
      ord <- seq(1:n) # initial order for the path
      max.distance <- 15
      
      ## Run the loop   ################################################
      for(i in 1:M){
        cost1[i]=distance(ord)
        if (i %% 10 == 0)
        {
          cat(sprintf("Iteration %d of %d\n",i,M))
        }
        for(t in 1:(fixed.itr)){
          # propose to move
          a=sample(1:n,2) # randomly choosing two cities
          result=ord
          result[a[2]]=ord[a[1]] # swapping
          result[a[1]]=ord[a[2]] # swapping
          U=runif(1)
          if( 
            U < exp( (distance(ord)-distance(result))/temp )){
            ord<-result # accept/reject
            numaccept <- numaccept + 1 # to compute the acceptance rate +}
          }
        }
        templist[i] <- temp
        temp <- temp*0.999
        temp <- tempfactor*temp
      }
      
      temperature <- data.frame(1:M,templist)
      
      length <- 1:length(cost1)
      fitting <- data.frame(length,cost1)
      
      pubs_inorder <- selected.pubs[ord,]
      
      map <- ggmap(basemap) +
        geom_segment(aes(xend=c(tail(lon, n=-1), NA), 
                         yend=c(tail(lat, n=-1), NA)),
                     data=pubs_inorder)
      
      result <- list(distance=cost1[M],
                     pubs_inorder=pubs_inorder,
                     fittingprocess=fitting,
                     temperature=temperature)
    }
    if(crow_distances==TRUE){
      ## Choose pubs you want #########################################
      selected.pubs <- subset(pubs, name %in% v_pubs)
      
      ## Variables to set #############################################
      n <- as.numeric(nrow(selected.pubs)) # number of pubs
      M <- 4000 # number of iterations
      temp <- 100 # initial temperature
      finaltemp <- 0.1
      tempfactor <- (temp/finaltemp)^(-1/M)
      fixed.itr <- 20 # number of iteration for a fixed temperature 
      
      ## FUNCTIONS ####################################################
      
      ## Check if integer 
      check.integer <- function(x) {
        x == round(x)
      }
      
      ## Difference between 2 points
      jb_difference <- function(lat_1,lon_1,lat_2,lon_2){
        rad <- pi/180
        a1 <- lat_1*rad
        a2 <- lon_1*rad
        b1 <- lat_2*rad
        b2 <- lon_2*rad
        dlon <- b2 - a2
        dlat <- b1 - a1
        a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
        c <- 2 * atan2(sqrt(a), sqrt(1 - a))
        R <- 6378.14500001
        d <- R * c
        return(d)
      }
      
      distance <- function(ord){
        p <- length(ord)
        d <- rep(0,p)
        for(i in 1:(p-1)){
          d[i] <- jb_difference(selected.pubs$lat[ord[i]],
                                selected.pubs$lon[ord[i]],
                                selected.pubs$lat[ord[i+1]],
                                selected.pubs$lon[ord[i+1]])
        }
        d[p] <- jb_difference(selected.pubs$lat[ord[1]],
                              selected.pubs$lon[ord[1]],
                              selected.pubs$lat[ord[p]],
                              selected.pubs$lon[ord[p]] )
        distance=sum(d)
        return(distance) 
      }
      
      ## Initial values ################################################
      
      cost1 <- rep(0,M)
      numaccept <- 0 # to keep track of accepting
      templist <- rep(0,M) # to keep track of temperature
      ord <- seq(1:n) # initial order for the path
      max.distance <- 15
      
      ## Run the loop   ################################################
      for(i in 1:M){
        cost1[i]=distance(ord)
        for(t in 1:(fixed.itr)){
          # propose to move
          a=sample(1:n,2) # randomly choosing two cities
          result=ord
          result[a[2]]=ord[a[1]] # swapping
          result[a[1]]=ord[a[2]] # swapping
          U=runif(1)
          if( 
            U < exp( (distance(ord)-distance(result))/temp )){
            ord<-result # accept/reject
            numaccept <- numaccept + 1 # to compute the acceptance rate +}
          }
        }
        templist[i] <- temp
        temp <- temp*0.999
        temp <- tempfactor*temp
      }
      
      temperature <- data.frame(1:M,templist)
      
      length <- 1:length(cost1)
      fitting <- data.frame(length,cost1)
      
      pubs_inorder <- selected.pubs[ord,]
      
      map <- ggmap(basemap) +
        geom_segment(aes(xend=c(tail(lon, n=-1), NA), 
                         yend=c(tail(lat, n=-1), NA)),
                     data=pubs_inorder)
      
      result <- list(distance=cost1[M],
                     pubs_inorder=pubs_inorder,
                     fittingprocess=fitting,
                     temperature=temperature)
    }
  }
  
  if(cam_pubs==FALSE){
    
    ## Variables to set #############################################
    n <- as.numeric(nrow(v_pubs)) # number of pubs
    M <- 4000 # number of iterations
    temp <- 100 # initial temperature
    finaltemp <- 0.1
    tempfactor <- (temp/finaltemp)^(-1/M)
    fixed.itr <- 20 # number of iteration for a fixed temperature 
    
    ## FUNCTIONS ####################################################
    
    ## Check if integer 
    check.integer <- function(x) {
      x == round(x)
    }
    
    ## Difference between 2 points
    jb_difference <- function(lat_1,lon_1,lat_2,lon_2){
      rad <- pi/180
      a1 <- lat_1*rad
      a2 <- lon_1*rad
      b1 <- lat_2*rad
      b2 <- lon_2*rad
      dlon <- b2 - a2
      dlat <- b1 - a1
      a <- (sin(dlat/2))^2 + cos(a1) * cos(b1) * (sin(dlon/2))^2
      c <- 2 * atan2(sqrt(a), sqrt(1 - a))
      R <- 6378.14500001
      d <- R * c
      return(d)
    }
    
    distance <- function(ord){
      p <- length(ord)
      d <- rep(0,p)
      for(i in 1:(p-1)){
        d[i] <- jb_difference(v_pubs[ord[i],1],
                              v_pubs[ord[i],2],
                              v_pubs[ord[i+1],1],
                              v_pubs[ord[i+1],2])
      }
      d[p] <- jb_difference(v_pubs[ord[1],1],
                            v_pubs[ord[1],2],
                            v_pubs[ord[p],1],
                            v_pubs[ord[p],2] )
      distance=sum(d)
      return(distance) 
    }
    
    ## Initial values ################################################
    
    cost1 <- rep(0,M)
    numaccept <- 0 # to keep track of accepting
    templist <- rep(0,M) # to keep track of temperature
    ord <- seq(1:n) # initial order for the path
    max.distance <- 15
    
    ## Run the loop   ################################################
    for(i in 1:M){
      cost1[i]=distance(ord)
      for(t in 1:(fixed.itr)){
        # propose to move
        a=sample(1:n,2) # randomly choosing two cities
        result=ord
        result[a[2]]=ord[a[1]] # swapping
        result[a[1]]=ord[a[2]] # swapping
        U=runif(1)
        if( 
          U < exp( (distance(ord)-distance(result))/temp )){
          ord<-result # accept/reject
          numaccept <- numaccept + 1 # to compute the acceptance rate +}
        }
      }
      templist[i] <- temp
      temp <- temp*0.999
      temp <- tempfactor*temp
    }
    
    temperature <- data.frame(1:M,templist)
    
    length <- 1:length(cost1)
    fitting <- data.frame(length,cost1)
    
    pubs_inorder <- v_pubs[ord,]
    
    map <- ggmap(basemap) +
      geom_segment(aes(xend=c(tail(lon, n=-1), NA), 
                       yend=c(tail(lat, n=-1), NA)),
                   data=pubs_inorder)+
      xlab("Latitude")+
      ylab("Longitude")
    
    result <- list(distance=cost1[M],
                   pubs_inorder=pubs_inorder,
                   fittingprocess=fitting,
                   temperature=temperature)
  }
  
  print(map)
  return(result)
  
}