calc_dist_mat <- function(df) {
  n = nrow(df)
  dist.mat = matrix(0,nrow = n,ncol = n)
  for(i in 1:n){
    for(j in 1:n){
      dist.mat[i,j] = (sqrt((df[i, "x"] - df[j, "x"])^2 + 
                            (df[i, "y"] - df[j, "y"])^2 ) %>% as.numeric())
    }
  }
  return(round(dist.mat,4))
}

read_solomon <- function(instance, custno = 25, write.data = T, st = 1){
  
  x = readLines(paste0("C:/Users/", computer, "/Desktop/solomon_25/C",
                       instance, ".txt"))
  x = x[-(1:7)]
  x = x[-c(2,29)]
  x = trimws(x)
  x = gsub("\\s+", " ", x)
  df = data.frame(matrix(0,26,7))
  for (i in 2:length(x) ) {
    df[i-1,1:7] <- read.table(textConnection(x[[i]]))
  }
  colnames(df) = c("n", "x", "y", "demand", "et", "lt", "st")
  df[27,] = df[1,]
  if(custno != 25){
    n = sample(2:26, custno)
    n = sort(n)
    df2 = df[c(1,n,27),]
    df2[,1] = c(0,1:(custno+1))
  }else{
    df2 = df
    df2[27,1] = 26
  }
  df2[-c(1,custno+2), "st"] = st
  if(write.data == TRUE){
    xlsx::write.xlsx(df2, paste0("C:/Users/", computer,
                                 "/Documents/gamsdir/projdir/dataR.xlsx"),
                     row.names = F, showNA = F) 
  }
  return(df2)
}

create_dataR <- function(instance = 105, st = 1){
  x = readLines(paste0("C:/Users/", computer, "/Desktop/solomon_25/C",
                       instance, ".txt"))
  x = x[-(1:7)]
  x = x[-c(2,29)]
  x = trimws(x)
  x = gsub("\\s+", " ", x)
  df = data.frame(matrix(0,26,7))
  for (i in 2:length(x) ) {
    df[i-1,1:7] <- read.table(textConnection(x[[i]]))
  }
  colnames(df) = c("n", "x", "y", "demand", "et", "lt", "st")
  df[27,] = df[1,]
  df[27,1] = 26
  df$st[-c(1,27)] = st
  xlsx::write.xlsx(df, paste0("C:/Users/", computer,
                              "/Documents/gamsdir/projdir/dataR.xlsx"),
                   row.names = F, showNA = F) 
  find_drone_eligible(df)
  plot(df$x, df$y)
}

find_drone_eligible <- function(df, path = NULL){
  if(is.null(path)){
    path = paste0("C:/Users/", computer, "/Documents/gamsdir/projdir")
  }
  depot = dim(df)[1]
  all.nodes = 1:depot
  dist.mat = calc_dist_mat(df)
  
  rn = 0
  m = matrix("yes", depot^2, depot+1)
  m[,1] = paste(rep(all.nodes-1, each = depot), all.nodes-1, sep = ".")
  for(i in all.nodes){
    for(j in all.nodes){
      rn = rn + 1
      for(k in all.nodes){
        if((i %in% c(j,k,depot)) || (j %in% c(1,k,depot)) || (k == 1)){
          m[rn, k+1] = "no"
        }else if(drone.time[i,j] + drone.time[j,k] > drone.endurance * 60){
          m[rn, k+1] = "no"
        }
      }
    }
  }
  
  colnames(m) = c("", all.nodes-1)
  write.table(m, file = paste(path, "droneset.txt", sep = "/"),
              row.names = F, quote = F, sep = "\t")
  
  print(sum(m=="yes")/(depot^3))
}

plot_dronetruck_kvehicle_compare <- function(trs1, drs1, trs2, drs2){
  n = nrow(df)
  plot(x = df$x[-n], y = df$y[-n], type = "p", col = "black",
       xlab = "", ylab = "", pch = 21, bg = "black")
  #      xlim = c(min(df$x)-2, max(df$x)+2), ylim = c(min(df$y)-2, max(df$y)+2))
  points(df$x[1], df$y[1], pch = 15, col = "orange", lwd = 2, cex = 2)
  text(x = df$x[-n], y = df$y[-n], labels = 1:(n-1), pos = 3)
  n = length(trs1)
  for(v in 1:n){
    tr = trs1[[v]]
    tr2 = trs2[[v]]
    dr = drs1[[v]]
    dr2 = drs2[[v]]
    
    x1 = paste(tr[-length(tr)], tr[-1], sep="-")
    x2 = paste(tr2[-length(tr2)], tr2[-1], sep="-")
    
    
    
    plotdf = df[tr[-length(tr)], 2:3]
    plotdf[,3:4] = df[tr[-1], 2:3]
    arrows(plotdf$x, plotdf$y, plotdf$x.1, plotdf$y.1, lwd = 2,
           lty = 1, length = 0.1, col = v)
    if(sum(dr) != 0){
      arrows(df$x[dr[,1]], df$y[dr[,1]], df$x[dr[,2]], df$y[dr[,2]],
             col = v, lwd = 2, lty = 2, length = 0.15)
      arrows(df$x[dr[,2]], df$y[dr[,2]], df$x[dr[,3]], df$y[dr[,3]],
             col = v, lwd = 2, lty = 2, length = 0.15)
      both = dr[dr[,3] %in% dr[,1],3]
      points(df$x[both], df$y[both], pch = 23, bg = "blue", cex = 2)
      points(df$x[dr[!(dr[,3]%in%both),3]], df$y[dr[!(dr[,3]%in%both),3]],
             pch = 25, bg = "green", cex = 1.5)
      points(df$x[dr[!(dr[,1]%in%both),1]], df$y[dr[!(dr[,1]%in%both),1]],
             pch = 24, bg = "red", cex = 1.5)
    }
  }
}

plot_dronetruck_kvehicle <- function(truck.routes, drone.routes = NULL){
  n = nrow(df)
  plot(x = df$x[-n], y = df$y[-n], type = "p", col = "black",
       xlab = "", ylab = "", pch = 21, bg = "black")
  #      xlim = c(min(df$x)-2, max(df$x)+2), ylim = c(min(df$y)-2, max(df$y)+2))
  points(df$x[1], df$y[1], pch = 15, col = "orange", lwd = 2, cex = 2)
  text(x = df$x[-n], y = df$y[-n], labels = 1:(n-1), pos = 3)
  n = length(truck.routes)
  for(v in 1:n){
    tr = truck.routes[[v]]
    if(is.null(drone.routes[[v]])){
      dr = matrix(0, nrow = 1, ncol = 3)
    }else{
      dr = drone.routes[[v]]
    }
    
    plotdf = df[tr[-length(tr)], 2:3]
    plotdf[,3:4] = df[tr[-1], 2:3]
    arrows(plotdf$x, plotdf$y, plotdf$x.1, plotdf$y.1, lwd = 2,
           lty = 1, length = 0.1, col = v)
    if(sum(dr) != 0){
      arrows(df$x[dr[,1]], df$y[dr[,1]], df$x[dr[,2]], df$y[dr[,2]],
             col = v, lwd = 2, lty = 2, length = 0.15)
      arrows(df$x[dr[,2]], df$y[dr[,2]], df$x[dr[,3]], df$y[dr[,3]],
             col = v, lwd = 2, lty = 2, length = 0.15)
      both = dr[dr[,3] %in% dr[,1],3]
      points(df$x[both], df$y[both], pch = 23, bg = "blue", cex = 2)
      points(df$x[dr[!(dr[,3]%in%both),3]], df$y[dr[!(dr[,3]%in%both),3]],
             pch = 25, bg = "green", cex = 1.5)
      points(df$x[dr[!(dr[,1]%in%both),1]], df$y[dr[!(dr[,1]%in%both),1]],
             pch = 24, bg = "red", cex = 1.5)
    }
  }
}

plot_dronetruck_kvehicle2 <- function(truck.routes, drone.routes = NULL){
  n = nrow(df)
  plot(x = df$x[-n], y = df$y[-n], type = "p", col = "black",
       xlab = "", ylab = "", pch = 21, bg = "black", axes = FALSE, frame.plot = TRUE)
  #      xlim = c(min(df$x)-2, max(df$x)+2), ylim = c(min(df$y)-2, max(df$y)+2))
  points(df$x[1], df$y[1], pch = 15, col = "orange", lwd = 2, cex = 2)
  n = length(truck.routes)
  for(v in 1:n){
    tr = truck.routes[[v]]
    if(is.null(drone.routes[[v]])){
      dr = matrix(0, nrow = 1, ncol = 3)
    }else{
      dr = drone.routes[[v]]
    }
    
    plotdf = df[tr[-length(tr)], 2:3]
    plotdf[,3:4] = df[tr[-1], 2:3]
    arrows(plotdf$x, plotdf$y, plotdf$x.1, plotdf$y.1, lwd = 2,
           lty = 1, length = 0.1, col = v)
    if(sum(dr) != 0){
      arrows(df$x[dr[,1]], df$y[dr[,1]], df$x[dr[,2]], df$y[dr[,2]],
             col = v, lwd = 2, lty = 2, length = 0.15)
      arrows(df$x[dr[,2]], df$y[dr[,2]], df$x[dr[,3]], df$y[dr[,3]],
             col = v, lwd = 2, lty = 2, length = 0.15)
    }
  }
}


plot_dronetruck <- function(tr, dr = matrix(0,1,3)){
  n = nrow(df)
  custno = n - 2
  plotdf = df[tr[-length(tr)], 2:3]
  plotdf[,3:4] = df[tr[-1], 2:3]
  plot(x = df$x[-n], y = df$y[-n], type = "p", col = "black",
       xlab = "", ylab = "", pch = 21, bg = "black", cex=2, ylim=c(0,42))
  #      xlim = c(min(df$x)-2, max(df$x)+2), ylim = c(min(df$y)-2, max(df$y)+2))
  points(df$x[1], df$y[1], pch = 15, col = "brown", lwd = 2, cex = 3)
  arrows(plotdf$x, plotdf$y, plotdf$x.1, plotdf$y.1, lwd = 2,
         lty = 1, length = 0.2)
  text(x = df$x[-n], y = df$y[-n], labels = 1:(custno+1), pos = 3)
  title(sub = paste(tr, collapse = "-"))
  if(sum(dr) != 0){
    arrows(df$x[dr[,1]], df$y[dr[,1]], df$x[dr[,2]], df$y[dr[,2]],
           col = "blue", lwd = 2, lty = 2, length = 0.2)
    arrows(df$x[dr[,2]], df$y[dr[,2]], df$x[dr[,3]], df$y[dr[,3]],
           col = "blue", lwd = 2, lty = 2, length = 0.2)
    both = dr[dr[,3] %in% dr[,1],3]
    points(df$x[both], df$y[both], pch = 23, bg = "orange", cex = 2)
    points(df$x[dr[!(dr[,3]%in%both),3]], df$y[dr[!(dr[,3]%in%both),3]],
           pch = 21, bg = "green", cex = 2)
    points(df$x[dr[!(dr[,1]%in%both),1]], df$y[dr[!(dr[,1]%in%both),1]],
           pch = 21, bg = "blue", cex = 2)
  }
}

plot_costvector <- function(){
  plot(cost.vector, type = "l", ylim = c(best.cost-100, initial.cost*1.2))
  points(x=which.min(cost.vector), y=best.cost, col = 2, pch = 18, cex = 1.5)
  abline(v=iline, col = "green", lw =2)
  abline(v=line1, col = "yellow", lw =2)
  abline(v=line2, col = "orange", lw =2)
  abline(v=line3, col = "red", lw =2)
  if(!include.drone.initial){
    abline(v=dline, col = "blue", lw =2)
  }
  title(main = paste("incumbent cost: ", round(best.cost,2)))
}

plot_sortie <- function(i,j,k){
  i = i + 1
  j = j + 1
  k = k + 1
  arrows(df$x[i], df$y[i], df$x[j], df$y[j],
         col = "blue", lwd = 2, lty = 1, length = 0.15)
  arrows(df$x[j], df$y[j], df$x[k], df$y[k],
         col = "blue", lwd = 2, lty = 1, length = 0.15)
  points(df$x[k], df$y[k], pch = 21, bg = "green")
  points(df$x[i], df$y[i], pch = 21, bg = "red")
}

generate_instance <- function(ncust = 20, st = 5, grid.size = 41, tw = TRUE,
                              mean.cust.time = 30, plot.df = TRUE, demand = NULL,
                              tw.randomized = TRUE, loose = NULL, dataR = FALSE,
                              depot.region = 0){
  x = rep(seq(0,40, length.out = grid.size), grid.size)
  y = rep(seq(0,40, length.out = grid.size), each = grid.size)
  n = length(x)
  et = sample(0:10, n, TRUE)
  midx = mean(x)
  midy = mean(y)
  if(tw == FALSE){
    lt = rep(10,n)
  }else{
    if(tw.randomized == TRUE){
      et = round(runif(n, 0, ncust*mean.cust.time))
      lt = et + round(runif(n, mean.cust.time*4, mean.cust.time*10))
    }else{
      a = ncust * mean.cust.time / 4
      for(i in 1:n){
        if(x[i]<midx & y[i]<midy){
          et[i] = sample(0:1, 1)
        }else if(x[i]<midx & y[i]>midy){
          et[i] = sample(2:3, 1)
        }else if(x[i]>midx & y[i]>midy){
          et[i] = sample(4:5, 1)
        }else{
          et[i] = sample(6:7, 1) 
        }
      }
      et = et*60
      lt = et + sample(2:6,n,T,2:6)*60 + sample(0:59,n,T)
    }
  }


  df = data.frame(n = 1:n, x, y, demand = 0, et = et, lt = lt)
  if(depot.region==1){
    depot = sample(df$n[df$x<midx & df$y<midy],1)
  }else{
    depot = sample(df$n, 1)
  }
  
  slist = 1:n
  slist = slist[slist!=depot]
  df2 = df[sample(slist,ncust+2),]
  df2[1,] = c(0, df$x[depot], df$y[depot], 0, 0, max(df$lt) + 2*mean.cust.time)
  df2[ncust+2,] = c(0, df$x[depot], df$y[depot], 0, 0, max(df$lt) + 2*mean.cust.time)
  df2$n = 0:(ncust+1)
  df2$st = c(0,rep(st,ncust),0)
  if(is.null(demand)){
    df2$demand = c(0, rep(1, ncust), 0)
  }else{
    df2$demand = c(0, sample(demand, ncust, TRUE), 0)
  }

  if(!is.null(loose)){
    df2$et = df2$et / loose
    df2$lt = df2$lt * loose
  }
  if(plot.df == TRUE) { 
    plot(df2$x[-1], df2$y[-1])
    points(df2$x[1], df2$y[1], pch = 18, cex = 2)
    title(xlab = "x", ylab = "y")
  }
  if(dataR == TRUE){
    xlsx::write.xlsx(df2, paste0("C:/Users/", computer,
                                "/Documents/gamsdir/projdir/dataR.xlsx"),
                     row.names = F, showNA = F)
    find_non_drone_eligible(df2)
  }
  return(df2)
}

generate_instance_centered <- function(ncust = 20, st = 5, grid.size = 41, max.d = 10,
                                       tw = TRUE, mean.cust.time = 30, plot.df = TRUE,
                                       demand = NULL, tw.randomized = TRUE,
                                       loose = NULL, dataR = FALSE, depot.region = 0){
  x = rep(seq(0,40, length.out = grid.size), grid.size)
  y = rep(seq(0,40, length.out = grid.size), each = grid.size)
  n = length(x)
  et = sample(0:10, n, TRUE)
  midx = mean(x)
  midy = mean(y)
  if(tw == FALSE){
    lt = rep(10,n)
  }else{
    if(tw.randomized == TRUE){
      et = round(runif(n, 0, ncust*mean.cust.time))
      lt = et + round(runif(n, mean.cust.time*4, mean.cust.time*10))
    }else{
      a = ncust * mean.cust.time / 4
      for(i in 1:n){
        if(x[i]<midx & y[i]<midy){
          et[i] = sample(0:1, 1)
        }else if(x[i]<midx & y[i]>midy){
          et[i] = sample(2:3, 1)
        }else if(x[i]>midx & y[i]>midy){
          et[i] = sample(4:5, 1)
        }else{
          et[i] = sample(6:7, 1) 
        }
      }
      et = et * 60
      lt = et + sample(4:6,n,T,4:6)*60  + sample(0:59,n,T)
    }
  }
  
  
  df = data.frame(n = 1:n, x, y, demand = 0, et = et, lt = lt)
  if(depot.region==1){
    depot = sample(df$n[df$x<midx & df$y<midy],1)
  }else{
    depot = sample(df$n, 1)
  }
  dx = df$x[depot]
  dy = df$y[depot]
  slist = 1:n
  slist = slist[slist!=depot]
  
  center = sample(slist, 1)
  cx = df$x[center]
  cy = df$y[center]
  df2 = df %>% filter((x>cx-max.d & x<cx+max.d), (y>cy-max.d & y<cy+max.d))
  
  if(nrow(df2)<ncust+2){
    while(nrow(df2)<ncust+2){
      if(depot.region==1){
        depot = sample(df$n[df$x<midx & df$y<midy],1)
      }else{
        depot = sample(df$n, 1)
      }
      dx = df$x[depot]
      dy = df$y[depot]
      slist = 1:n
      slist = slist[slist!=depot]
      center = sample(slist, 1)
      cx = df$x[center]
      cy = df$y[center]
      df2 = df %>% filter((x>cx-max.d & x<cx+max.d), (y>cy-max.d & y<cy+max.d))

    }
    print(nrow(df2))
  }
  df3 = df2[sample(1:nrow(df2),ncust+2),]
  df4 = df %>% filter((x>cx-max.d & x<cx+max.d),
                      (y>cy-max.d & y<cy+max.d),
                      !(n%in%df3$n))
  df3[1,] = c(0, dx, dy, 0, 0, max(df2$lt) + 2*mean.cust.time)
  df3[ncust+2,] = c(0, dx, dy, 0, 0, max(df2$lt) + 2*mean.cust.time)
  df3$n = 0:(ncust+1)
  df3$st = c(0,rep(st,ncust),0)
  if(is.null(demand)){
    df3$demand = c(0, rep(1, ncust), 0)
  }else{
    df3$demand = c(0, sample(demand, ncust, TRUE), 0)
  }
  
  if(!is.null(loose)){
    df3$et = df3$et / loose
    df3$lt = df3$lt * loose
  }
  if(plot.df == TRUE) { 
    plot(df3$x[-1], df3$y[-1])
    points(df3$x[1], df3$y[1], pch = 18, cex = 2)
    title(xlab = "x", ylab = "y")
  }
  if(dataR == TRUE){
    xlsx::write.xlsx(df3, paste0("C:/Users/", computer,
                                 "/Documents/gamsdir/projdir/dataR.xlsx"),
                     row.names = F, showNA = F)
    find_non_drone_eligible(df3)
  }
  return(df3)
}

find_dist <- function(i,j,k){
  return(dist.mat[i,j]+dist.mat[j,k])
}

read_gams_sol_kvehicle <- function(path = "~/gamsdir/projdir", g = 1, drone = T,
                                   read.df = T, plot.route = T, vrptw = F){
  if(read.df){
    df <<- readxl::read_excel(paste0(path, "/dataR.xlsx"))
    dist.mat <<- calc_dist_mat(df)
    truck.time <<- round(dist.mat/truck.velocity,4)
    drone.time <<- round(dist.mat/drone.velocity,4)
    nodes <<- nrow(df)
    ncust <<- nodes - 2
  }
  if(vrptw){
    drone = F
    rtr = "vrptwroute"
  }else{
    rtr = "resultstr"
  }

  rdr = "resultsdr"
  gams.res = readxl::read_excel(paste0(path, "/", rtr, ".xls"))
  v = ncol(gams.res)-2
  gams.tr = vector("list", v)
  gams.dr = vector("list", v)

  colnames(gams.res) = c("from", "to", paste0("v",1:v))
  for(k in 1:v){
    x = !is.na(gams.res[,k+2])
    y = gams.res[x,]
    tr = character(sum(x)+1)
    tr[1] = "0"
    tr[2] = unlist(y[1,2])
    for(i in 2:sum(x)){
      j = which(y[,1] == tr[i])
      tr[i+1] = unlist(y[j,2])
    }
    tr[tr==ncust+1] = 0
    gams.tr[[k]] = as.numeric(tr)+1
    gams.dr[[k]] = matrix(0, 1, 3)
  }



  if(drone){
    gams.dr = vector("list", v)
    gams.res.dr = readxl::read_excel(paste0(path, "/", rdr, ".xls"))
    drone.used.routes = ncol(gams.res.dr)-3
    if(v != drone.used.routes){
      inds = as.numeric(colnames(gams.res.dr)[-c(1:3)])
    }else{
      inds = 1:v
    }
    j = 3
    for(k in 1:v){
      if(k %in% inds){
        j = j + 1
        x = !is.na(gams.res.dr[,j])
        y = gams.res.dr[x,]
        dr = apply(y[1:3], 2, as.numeric) %>% unname()
        dr[dr==ncust+1]=0
      }else{
        dr = matrix(-1,1,3)
      }
      if(length(dr)==3){dr = matrix(dr,1,3)}
      gams.dr[[k]] = dr+1
    }
  }
  
  f = is_solution_feasible2(gams.tr, gams.dr)
  print(f)
  if(plot.route){
    plot_dronetruck_kvehicle(gams.tr, gams.dr)
    title(main = paste("optimal cost:",
                       find_solution_cost(gams.tr, gams.dr, "sum")))
  }
  
  gams.tr <<- gams.tr
  gams.dr <<- gams.dr
  return(f)
}

write_vrp_sol <- function(path = "~/gamsdir/projdir", read.df = T){
  if(read.df){
    df <<- readxl::read_excel(paste0(path, "/dataR.xlsx"))
    dist.mat <<- calc_dist_mat(df)
    nodes <<- nrow(df)
    ncust <<- nodes - 2
  }
  read_gams_sol_kvehicle(path, plot.route = F, vrptw = T)
  k = length(gams.tr)
  a = character()
  x = 1
  for(v in 1:k){
    j = length(gams.tr[[v]])
    for(i in 1:(j-1)){
      to = gams.tr[[v]][i+1]
      a[x] = paste('x.l(', gams.tr[[v]][i]-1, ',',
                           ifelse(to==1,ncust+1,to-1), ',', v, ')', sep = "'")
      x = x+1
    }
  }

  write.table(sapply(1:length(a), function(x) paste0(a[x], " = 1;")),
              file = paste0("C:/Users/", computer,
                            "/Documents/gamsdir/projdir/vrproute.txt"),
              row.names = F, col.names = F, quote = F)
}

write_heuristic_sol <- function(trs, drs = matrix(0,1,3), type = "fx", path = NULL, soltype = "unique"){
  if(is.null(path)){
    path = paste0("C:/Users/", computer, "/Documents/gamsdir/projdir")
  }
  
  a = character()
  b = character()
  if(soltype=="unique"){
    tr = trs
    dr = drs
    tr = tr - 1
    tr[length(tr)] = nrow(df)-1
    for(i in 1:(length(tr)-1)){
      a[i] = paste(paste0('x.', type, '('), tr[i], ',', tr[i+1], ',', 1, ')', sep = "'")
    }
    if(sum(dr)!=0){
      dr = dr - 1
      for(i in 1:nrow(dr)){
        if(dr[i,3]!=0){
          cc = dr[i,3]
        }else{
          cc = nrow(df)-1
        }
        b[i] = paste(paste0('y.', type, '('), dr[i,1], ',', dr[i,2], ',',
                     cc, ',', 1, ')', sep = "'")
      }
    }
  }else{
    j = 1
    k = 1
    for(v in 1:length(trs)){
      tr = trs[[v]]
      dr = drs[[v]]
      tr = tr - 1
      tr[length(tr)] = nrow(df)-1
      for(i in 1:(length(tr)-1)){
        a[j] = paste(paste0('x.', type, '('), tr[i], ',', tr[i+1], ',', v, ')', sep = "'")
        j = j + 1
      }
      if(sum(dr)!=0){
        dr = dr - 1
        for(i in 1:nrow(dr)){
          if(dr[i,3]!=0){
            cc = dr[i,3]
          }else{
            cc = nrow(df)-1
          }
          b[k] = paste(paste0('y.', type, '('), dr[i,1], ',', dr[i,2], ',',
                       cc, ',', v, ')', sep = "'")
          k = k + 1
        }
      }
    }
  }
  
  write.table(sapply(1:(length(tr)-1), function(x) paste0(a[x], " = 1;")),
              file = paste(path, "heuristic_tr.txt", sep = "/"),
              row.names = F, col.names = F, quote = F)
  write.table(sapply(1:nrow(dr), function(x) paste0(b[x], " = 1;")),
              file = paste(path, "heuristic_dr.txt", sep = "/"),
              row.names = F, col.names = F, quote = F)
  
}

write_gams_files<- function(v, ncust, path = NULL){
  if(is.null(path)){
    path = paste0("C:/Users/", computer, "/Documents/gamsdir/projdir")
  }
  
  x = character(5)
  x[1] = "Set"
  x[2] = paste0("a 'all nodes' / 0*", ncust+1, " /")
  x[3] = paste0("i(a) 'start nodes' / 0*", ncust," /")
  x[4] = paste0("j(a) 'start nodes' / 1*", ncust," /")
  x[5] = paste0("k(a) 'start nodes' / 1*", ncust+1," /;")
  x = as.data.frame(x)
  write.table(x, paste(path, "sets.txt", sep = "/"),
              quote = F, row.names = F, col.names = F)
  
  vehicle = ifelse(v==1, "Set v /1/;", "Set v /1*2/;")
  write.table(vehicle, paste(path, "vehicle.txt",sep="/"), quote = F, row.names = F, col.names = F)
  
  depot = character(2)
  depot[1] = paste('truck_return_depot(v).. sum(j, x(j, ',ncust+1,', v)) =e= 1;',sep = '"')
  depot[2] = paste('truck_depot_no_leave(v).. sum(a, x(',ncust+1,', a, v)) =e= 0;',sep = '"')
  write.table(depot, paste(path, "depot_equation.txt", sep="/"),
              quote = F, row.names = F, col.names = F)
  
  coord = character(2)
  coord[1] = paste('recovery_time_coord1(k,v)$(not sameas(k,', ncust+1,
                   ')).. atd(k,v) =g= att(k,v) - bigM * (1 - sum((i,j)$droneset(i,j,k), y(i,j,k,v)));',
                   sep = '"')
  coord[2] = paste('recovery_time_coord2(k,v)$(not sameas(k,', ncust+1,
                   ')).. atd(k,v) =l= att(k,v) + bigM * (1 - sum((i,j)$droneset(i,j,k), y(i,j,k,v)));',
                   sep = '"')
  write.table(coord, paste(path, "recovery_equations.txt", sep = "/"),
              quote = F, row.names = F, col.names = F)
  
}

read_optimal <- function(dtset, where="gams", sol="best", file.name="vrpdk"){
  sol = ifelse(sol=="best", "Best", "Final")
  path = paste(wd.path, where, paste0("gams-tw-", dtset), paste0(file.name,".log"), sep = "/")
  con = file(path, "r")
  aa = readLines(con)
  close(con)
  for(j in length(aa):1){
    a = unlist(strsplit(aa[j], " "))[1]
    if(!is.na(a) & a[1] == sol){
      if(sol == "Best"){
        opt = unlist(strsplit(aa[j], " ")) %>% last() %>% as.numeric()
      }else{
        b = suppressWarnings(as.numeric(unlist(strsplit(aa[j], " "))))
        opt = b[!is.na(b)]
      }
      break()
    }
  }
  return(opt)
}

calc_stats <- function(df.formula, perc = TRUE, r = 2){
  mult = ifelse(perc, 100 ,1)
  res.df = data.frame(dtset = dtsets,
                      min = round(apply(eval(df.formula), 1, min), r) * mult,
                      mean = round(apply(eval(df.formula), 1, mean), r) * mult,
                      median = round(apply(eval(df.formula), 1, median), r) * mult,
                      max = round(apply(eval(df.formula), 1, max), r) * mult,
                      stringsAsFactors = F)
  meanrow = c("avg", round(apply(res.df[,-1], 2, mean), r))
  return(rbind(res.df, meanrow))
}

save_stats <- function(type = NULL, moves = TRUE){
  if(!is.null(type)){
    save.name = file.choose()
    resultsdf = xlsx::read.xlsx(save.name, 1)
  }
  # resultsdf = xlsx::read.xlsx("resultsdf-1mins-1.xlsx", 1)
  resultsdf2 = as.data.frame(resultsdf)
  colnames(resultsdf2) = c("dtset", paste0(rep("res"),1:nstart),
                           paste0(rep("initial"),1:nstart),
                           paste0(rep("times"),1:nstart),
                           paste0(rep("seed"),1:nstart),
                           paste0(rep("vrp"),1:nstart))
  resultsdf2$dtset = dtsets[dtlen]
  resultsdf2$opt = opt
  
  res <<- resultsdf2[,2:(nstart+1)]
  init <<- resultsdf2[,(nstart+2):(2*nstart+1)]
  times <<- resultsdf2[,(2*nstart+2):(3*nstart+1)]
  vrp <<- resultsdf2[,(4*nstart+2):(5*nstart+1)]
  
  res.df = calc_stats(quote(res/opt-1), T)
  
  print(res.df)
  xlsx::write.xlsx(res.df, save.name, row.names = F, showNA = F,
                   append = T, sheetName = "final")
  
  init.imp = calc_stats(quote(1-res/init), T)
  xlsx::write.xlsx(init.imp, save.name, row.names = F, showNA = F,
                   append = T, sheetName = "initial")
  
  times = calc_stats(quote(times), F)
  xlsx::write.xlsx(times, save.name, row.names = F, showNA = F,
                   append = T, sheetName = "times")
  
  vrp = calc_stats(quote(1-res/vrp), F)
  xlsx::write.xlsx(vrp, save.name, row.names = F, showNA = F,
                   append = T, sheetName = "vrp")
  
  if(moves==TRUE){
    movesdf2 = as.data.frame(movesdf)
    colnames(movesdf2) = c("dtset", paste0(rep("used"),1:max.move),
                           paste0(rep("feasible"),1:max.move),
                           paste0(rep("better"),1:max.move),
                           paste0(rep("incumbent"),1:max.move))
    movesdf2$dtset = dtsets
    xlsx::write.xlsx(movesdf, paste0(save.name2, collapse = "-"),
                     row.names = F, showNA = F, sheetName = "all") 
  }
}

save_solution <- function(trs, drs, path = NULL, data.folder = NULL, worst = T, where = "gams", opt = FALSE){
  if(is.null(path)) path = wd.path
  #folder.name = ifelse(smethod, paste(method, start.method, sep="-"), method)
  folder.name = paste0("heuristic/", start.method)
  initial.path = paste(path, folder.name, sep = "/")
  if(!dir.exists(initial.path)) dir.create(initial.path)
  new.path = paste(initial.path, paste("alns", dtset, sep = "-"), sep = "/")
  if(!dir.exists(new.path)){
    dir.create(new.path)
    file.copy(from = paste(data.folder, "dataR.xlsx", sep = "/"), to = new.path)
    write_gams_files(g, ncust, new.path)
    find_drone_eligible(df, new.path)
  }
  
  if(ns==1){
    # save the first solution as both best and worst
    save_plots(trs, drs, new.path, "best")
    save_plots(trs, drs, new.path, "worst")
    save_heuristic_route(trs, drs, new.path)
    write_heuristic_sol(trs, drs, "fx", new.path, "list")
    # save the plot of optimal solution from gams run
    if(opt){
      read_gams_sol_kvehicle(paste0(wd.path, "/", where ,"/gams-", dtset),
                             read.df = F, plot.route = F)
      png(filename = paste0(new.path, "/optimal_solution.png"), width = 600, height = 600)
      plot_dronetruck_kvehicle(gams.tr, gams.dr)
      title(main = paste("optimal cost:", find_solution_cost(gams.tr, gams.dr, "sum")))
      dev.off()
    }
  }else if(round(best.cost,2)<min(nresults[1:ns])){
    save_plots(trs, drs, new.path, "best")
    save_heuristic_route(trs, drs, new.path)
    write_heuristic_sol(trs, drs, "fx", new.path, "list")
  }else if(round(best.cost,2)>max(nresults[1:ns])){
    save_plots(trs, drs, new.path, "worst")
  }
}

save_plots <- function(trs, drs, path, type = "best"){
  # plotting and saving of solution
  sol.plot.name = paste0(path, "/", type, "_solution", ".png")
  if(file.exists(sol.plot.name)) file.remove(sol.plot.name)
  png(filename = sol.plot.name, width = 600, height = 600)
  plot_dronetruck_kvehicle(trs, drs)
  if(exists("opt")){
    title(main = paste(type, "cost:", round(best.cost,2),
                       ", perc dev:", round(best.cost/opt[data.ind]-1,4)*100))
  }else{
    title(main = paste(type, "cost:", round(best.cost,2)))
  }
  dev.off()
  
  # plotting and saving of accepted costs of simulated annealing
  cost.plot.name = paste0(path, "/", type, "_cost_vector", ".png")
  if(file.exists(cost.plot.name)) file.remove(cost.plot.name)
  png(filename = cost.plot.name, width = 600, height = 480)
  plot_costvector()
  dev.off()
}

save_heuristic_route <- function(trs, drs, path){
  #write_heuristic_sol(tr, dr, path)
  saveRDS(trs, paste(path, "xtr.RDS", sep = "/"))
  saveRDS(drs, paste(path, "xdr.RDS", sep = "/"))
}

update_matrices <- function(d){
  # for stat savings after completion of an instance run
  resultsdf[d,1] <<- d
  resultsdf[d,2:(1+nstart)] <<- nresults
  resultsdf[d,(nstart+2):(1+2*nstart)] <<- ninitials
  resultsdf[d,(2*nstart+2):(1+3*nstart)] <<- ntimes
  resultsdf[d,(3*nstart+2):(1+4*nstart)] <<- seeds
  resultsdf[d,(4*nstart+2):(1+5*nstart)] <<- nvrp
  
  movesdf[d,1] <<- d
  movesdf[d,2:(1+max.move)] <<- movesdf[d,2:(1+max.move)] + moves.used
  movesdf[d,(max.move+2):(1+2*max.move)] <<-
    movesdf[d,(max.move+2):(1+2*max.move)] + moves.feasible
  movesdf[d,(2*max.move+2):(1+3*max.move)] <<-
    movesdf[d,(2*max.move+2):(1+3*max.move)] + moves.successful
  movesdf[d,(3*max.move+2):(1+4*max.move)] <<-
    movesdf[d,(3*max.move+2):(1+4*max.move)] + moves.incumbent
}