In derek-corcoran-barrios/RtoAmpl:

knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)

multisppBiomass <- function(Stacklist, maxalpha = 5, maxbiomass = 2, maxcapacidad = 10, name = "Stack", Dist = 500000, Threshold, costlayer, bf = c(1,2)){

  Masklayer <- costlayer
  values(Masklayer) <- ifelse(is.na(values(Masklayer)), NA, 1)
  TempRaster <- list()
  for (i in 1:length(Stacklist)){
    Stacklist[[i]] <- Stacklist[[i]] * Masklayer
    TempStack <- Stacklist[[i]]

    values(TempStack)[values(TempStack) < Threshold] = 0
    values(TempStack)[values(TempStack) >= Threshold] = 1

    TempRaster[i] <- sum(TempStack)
  }
  TempRaster <- do.call("sum", TempRaster)
  TempRaster[values(TempRaster) > 0] = 1
  TempRaster[values(TempRaster) == 0] = NA
  for (i in 1:length(Stacklist)){
    Stacklist[[i]] <- Stacklist[[i]]*TempRaster
  }


  Alphas <- list()
  for(j in 1:length(Stacklist)){
    Alpha <- list()
    for (i in 1:nlayers(Stacklist[[j]])){
      temp <- data.frame(Alpha = values(Stacklist[[j]][[i]]), ID = 1:length(values(Stacklist[[j]][[i]])), Time = i-1)
      Alpha[[i]] <- temp[complete.cases(temp),]
    }
    Alphas[[j]]<- do.call("rbind", Alpha)
    Alphas[[j]]$Spp <- names(Stacklist)[j]
  }

  Alpha <-  do.call(rbind, Alphas)
  Alpha$Alpha <- Alpha$Alpha/max(Alpha$Alpha)*maxalpha
  s <- Alpha %>% group_by(ID) %>% summarise(SUMA = sum(Alpha)) %>% filter(SUMA > 0)
  Alpha <- Alpha[Alpha$ID %in% s$ID,]
  Spps <- unique(Alpha$Spp)
  Nodos <- unique(Alpha$ID)

  Alphas <- list()
  for (i in Spps){
    Alphas[[i]] <- dplyr::filter(Alpha, Spp == i)

    temp <-  split(Alphas[[i]], Alphas[[i]]$Time)
    Alphas[[i]] <- do.call(cbind, lapply(1:length(temp), function(i){
      if (i == 1){
        setNames(data.frame(paste("[",temp[[i]][["Spp"]],",",temp[[i]][["ID"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["Alpha"]),
                 c("V", paste("T", i, sep = ""), i-1))
      } else if (i == length(temp)){
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Alpha"], rep("\n", NROW(temp[[i]]))),
                 c(paste("T", i, sep = ""), i-1, "line"))
      } else {
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Alpha"]),
                 c(paste("T", i, sep = ""), i-1))
      }
    }))
  }

  Alpha <-do.call("rbind", Alphas)

  Biomasa <- list()

  for (i in 1:length(Stacklist)){
    Biomasa[[i]] <- data.frame(Sp = rep(Spps[i], times = length(values(Stacklist[[i]][[1]]))), ID = 1:length(values(Stacklist[[i]][[1]])), Biomasa = (values(Stacklist[[i]][[1]])/max(values(Stacklist[[i]][[1]]), na.rm = T))*maxbiomass)
    Biomasa[[i]] <-  Biomasa[[i]][complete.cases(Biomasa[[i]]),]
    Biomasa[[i]]$Biomasa <- round(Biomasa[[i]]$Biomasa, 4)
    Biomasa[[i]]$Sp <- Spps[i]
  }

  Biomasa <- do.call(rbind, Biomasa)
  Biomasa$line <- "\n"
  Biomasa$ID  <-  paste0("[",Biomasa$Sp, ",", Biomasa$ID, "]")
  Biomasa <- Biomasa[,-1]

  Capacidades <- list()
  for(j in 1:length(Stacklist)){
    Capacidad <- list()
    for (i in 1:nlayers(Stacklist[[j]])){
      temp <- data.frame(Capacidad = values(Stacklist[[j]][[i]]), ID = 1:length(values(Stacklist[[j]][[i]])), Time = i-1)
      Capacidad[[i]] <- temp[complete.cases(temp),]
    }
    Capacidades[[j]]<- do.call("rbind", Capacidad)
    Capacidades[[j]]$Spp <- names(Stacklist)[j]
  }

  Capacidad <-  do.call(rbind, Capacidades)
  Capacidad$Capacidad <- Capacidad$Capacidad/max(Capacidad$Capacidad)*maxcapacidad
  Capacidad <- Capacidad[Capacidad$ID %in% s$ID,]

  Capacidades <- list()
  for (i in Spps){
    Capacidades[[i]] <- dplyr::filter(Capacidad, Spp == i)

    temp <-  split(Capacidades[[i]], Capacidades[[i]]$Time)
    Capacidades[[i]] <- do.call(cbind, lapply(1:length(temp), function(i){
      if (i == 1){
        setNames(data.frame(paste("[",temp[[i]][["Spp"]],",",temp[[i]][["ID"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["Capacidad"]),
                 c("V", paste("T", i, sep = ""), i-1))
      } else if (i == length(temp)){
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Capacidad"], rep("\n", NROW(temp[[i]]))),
                 c(paste("T", i, sep = ""), i-1, "line"))
      } else {
        setNames(data.frame(temp[[i]]["Time"], temp[[i]]["Capacidad"]),
                 c(paste("T", i, sep = ""), i-1))
      }
    }))
  }

  Capacidad <-do.call("rbind", Capacidades)


  cost <-data.frame(ID = Nodos,cost = values(costlayer)[Nodos])
  cost$ID <- paste0("[", cost$ID, "]")
  cost$line <- "\n"

  BF <- data.frame(Spp = Spps, BF = bf, Space = "\n")


  Beta <- list()
  for(i in 1:length(Stacklist)){
    Beta[[i]] <- DistConect(Stacklist[[i]][[1]], Distance = Dist, Time = nlayers(Stacklist[[i]]))
    Beta[[i]] <- dplyr::filter(Beta[[i]], from %in% Nodos, to %in% Nodos)
    Beta[[i]]$Sp <- Spps[i]
  }

  Beta <- do.call(rbind, Beta)

  temp <-  split(Beta, Beta$Time)
  Betas <- do.call(cbind, lapply(1:length(temp), function(i){
    if (i == 1){
      setNames(data.frame(paste("[",temp[[i]][["Sp"]], ",",temp[[i]][["from"]], ",", temp[[i]][["to"]], ",*","]", sep = ""), temp[[i]]["Time"], temp[[i]]["beta"]),
               c("V", paste("T", i, sep = ""), i-1))
    } else if (i == length(temp)){
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["beta"], rep("\n", NROW(temp[[i]]))),
               c(paste("T", i, sep = ""), i-1, "line"))
    } else {
      setNames(data.frame(temp[[i]]["Time"], temp[[i]]["beta"]),
               c(paste("T", i, sep = ""), i-1))
    }
  }))


  sink(paste0(name, ".dat"))
  cat(c("set V :=", Nodos, ";"))
  cat("\n")
  cat(c("set Sp :=", names(Stacklist), ";"))
  cat("\n")
  cat(c("set E :=", paste0("(",unique(unite_(Beta, col = "V", sep = ",", from = c("from", "to"))$V), ")"), ";"))
  cat("\n")
  cat("param T:= 7;")
  cat("\n")
  cat("param alpha :=")
  cat("\n")
  cat(do.call(paste, Alpha))
  cat(";")
  cat("\n")
  cat("param beta :=")
  cat("\n")
  cat(do.call(paste, Betas))
  cat(";")
  cat("\n")
  cat("param b0 :=")
  cat("\n")
  cat(do.call(paste, Biomasa))
  cat(";")
  cat("\n")
  cat("param u :=")
  cat("\n")
  cat(do.call(paste, Capacidad))
  cat(";")
  cat("\n")
  cat("param bf := ")
  cat("\n")
  cat(do.call(paste, BF))
  cat(";")
  cat("\n")
  cat("param c :=")
  cat("\n")
  cat(do.call(paste, cost))
  cat(";")
  cat("\n")
  sink()
  return(list(Nodos = Nodos, Biomasa = Biomasa, Alphas = Alphas, Alpha = Alpha))
}
library(RtoAmpl)

library(RtoAmpl)
data(Cost)
data("univariate")
data("bivariate")

Stacklist <- list(univariate, bivariate)
names(Stacklist) <- c("SPA", "SPB")

setwd("/home/derek/Documents/PostdocPablo/AMPL")

for (i in 6:9){
  multisppBiomass(Stacklist = Stacklist, costlayer = Cost, Threshold = 0.2, name = paste0("bla","_", i), bf = c((2*i),i), maxalpha = 4)
}

cd
cd Documents/PostdocPablo/AMPL

./ampl

include RunLoop.run

setwd("/home/derek/Documents/PostdocPablo/AMPL")

SOLS1 <- list.files(pattern = "bla_[0-9].txt")
SOLS2 <-  list.files(pattern = "bla_[0-9][0-9].txt")
SOLS <- c(SOLS1, SOLS2)
library(readr)

bla <- list()
for (i in 1:length(SOLS)){
  DF <- read_delim(SOLS[i], " ", escape_double = FALSE, col_names = FALSE,  trim_ws = TRUE)

  colnames(DF) <- c("ID", "z")

  temp <- univariate[[1]]

  values(temp) <- NA
  values(temp)[DF$ID] <- DF$z
  bla[[i]] <- temp
}

Costs1 <-  list.files(pattern = "bla_[0-9]cost.txt")
Costs2 <-  list.files(pattern = "bla_[0-9][0-9]cost.txt")
Costs <- c(Costs1, Costs2)

DF <- list()
for(i in 1:length(Costs)){
  DF[[i]] <- read_delim(Costs[i],  " ", escape_double = FALSE, col_names = FALSE, trim_ws = TRUE)
  colnames(DF[[i]]) <- c("bf", "cost")
  DF[[i]]$SP <- c("A", "B")
}

DF <- do.call(rbind, DF)

library(ggplot2)
library(dplyr)
DF <- filter(DF, cost > 0)

ggplot(DF, aes(x = cost, y = bf)) + geom_line(aes(color = SP)) + geom_point(aes(color = SP)) + theme_classic()

setwd("/home/derek/Documents/RtoAmpl")

animation::saveGIF(for(i in 1:length(bla)){plot(bla[[i]])}, movie.name = "Yupi.gif")

derek-corcoran-barrios/RtoAmpl documentation built on May 14, 2019, 10:33 a.m.