## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = requireNamespace("vegan", quietly = TRUE)
)

## ----vegan, eval = !requireNamespace("vegan", quietly = TRUE), echo = FALSE, comment = NA----
# message('This vignette requires the "vegan" package. Please, install it: install.packages("vegan").')

## ----install, eval=FALSE------------------------------------------------------
# install.packages("ecoregime")
# devtools::install_github(repo = "MSPinillos/ecoregime", dependencies = T, build_vignettes = T)

## ----setup--------------------------------------------------------------------
library(ecoregime)

## ----citation-----------------------------------------------------------------
citation("ecoregime")

## ----data---------------------------------------------------------------------
# Matrix including the state variables (sp1-sp12) of the EDR trajectories
edr <- EDR_data$EDR1$abundance

# The first target is composed of one state resulting from averaging the state
# variables of two states in the reference EDR
target1 <- data.frame(matrix(colMeans(edr[traj == 3 & state %in% 1:2, paste0('sp', 1:12)]),
                             ncol = 12, 
                             dimnames = list(1, paste0('sp', 1:12))))
target1$traj <- 'target1'
target1$state <- 1

# The second target is composed of three states resulting from averaging the 
# state variables of four states in the reference EDR
target2 <- data.frame(t(sapply(1:3, function(istate){
  matrix(colMeans(edr[traj == 6 & state %in% istate:(istate+1), 
                      paste0('sp', 1:12)]))
})))
names(target2) <- paste0('sp', 1:12)
target2$traj <- 'target2'
target2$state <- 1:3

# For the third target, we will consider a trajectory of a different EDR
target3 <- EDR_data$EDR2$abundance[1:5, 3:ncol(EDR_data$EDR2$abundance)]
target3$traj <- 'target3'


## ----state_var----------------------------------------------------------------
# Select the columns containing the state variables (sp1, ..., sp12) and 
# include the information of EDR and target states in the same data.frame 
state_var1 <- data.frame(rbind(edr[, paste0("sp", 1:12)], 
                               target1[, paste0("sp", 1:12)]))
head(state_var1)

## ----trajectories1------------------------------------------------------------
trajectories1 <- c(edr$traj, target1$traj)
head(trajectories1)

## ----states1------------------------------------------------------------------
states1 <- as.integer(c(edr$state, target1$state))
head(states1)

## ----dStates1-----------------------------------------------------------------
# Compute state dissimilarities from state_var
dStates1 <- vegan::vegdist(x = state_var1, method = "bray")

## ----petra k, fig.width=8, fig.height=4, warning=FALSE------------------------
# Compute petra_edr using a small k
petra_k1 <- petra_edr(state_var = state_var1, 
                      trajectories = trajectories1,
                      states = states1, 
                      targets = "target1",
                      d_function = "vegan::vegdist", 
                      d_args = list(x = state_var1, method = "bray"),
                      k = 2L,
                      minPts = 2L,
                      return_args = T)

# Compute petra_edr using a large k
petra_k2 <- petra_edr(state_var = state_var1, 
                      trajectories = trajectories1,
                      states = states1, 
                      targets = "target1",
                      d_function = "vegan::vegdist", 
                      d_args = list(x = state_var1, method = "bray"),
                      k = 20L,
                      minPts = 2L,
                      return_args = T)

# Use plot to see the results
par(mfrow = c(1, 2))
plot(x = petra_k1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "k = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_k2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "k = 20")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2, 
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

## ----petra k1-----------------------------------------------------------------
petra_k1$k_dist

## ----petra k2-----------------------------------------------------------------
tail(petra_k2$k_dist)

## ----petra eps, fig.width=8, fig.height=4, warning=FALSE----------------------
# Compute petra_edr using a small eps
petra_eps1 <- petra_edr(state_var = state_var1, 
                        trajectories = trajectories1, 
                        states = states1,
                        targets = "target1",
                        d_function = "vegan::vegdist", 
                        d_args = list(x = state_var1, method = "bray"),
                        k = 20L, 
                        minPts = 2L, 
                        eps = 0.03, 
                        return_args = T)

# Compute petra_edr using a large eps
petra_eps2 <- petra_edr(state_var = state_var1, 
                        trajectories = trajectories1, 
                        states = states1, 
                        targets = "target1",
                        d_function = "vegan::vegdist", 
                        d_args = list(x = state_var1, method = "bray"),
                        k = 20L, 
                        minPts = 2L, 
                        eps = 0.1, 
                        return_args = T)

# Plot PETRA outputs
par(mfrow = c(1, 2))
plot(x = petra_eps1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "eps = 0.03")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_eps2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "eps = 0.1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")
         

## ----petra eps1---------------------------------------------------------------
petra_eps1$k_dist

## ----petra minPts, fig.width=8, fig.height=4, warning=FALSE-------------------
# Compute petra_edr using a small minPts
petra_minPts1 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 6L, 
                           minPts = 2L, 
                           return_args = T)

# Compute petra_edr using a large minPts
petra_minPts2 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 6L, 
                           minPts = 6L, 
                           return_args = T)

# Plot PETRA outputs
par(mfrow = c(1, 2))
plot(x = petra_minPts1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "minPts = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_minPts2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "minPts = 10")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = c(NA, 2, 2), 
       pch = c(20, NA, NA),
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

## ----petra minPts1 minPts2----------------------------------------------------
petra_minPts1$predicted_dist[, c("target", "predicted_state", "N")]
petra_minPts2$predicted_dist[, c("target", "predicted_state", "N")]

## ----method-------------------------------------------------------------------
# Compute petra_edr using method = "mean"
petra_method1 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1,
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 20L, 
                           minPts = 2L, 
                           method = "mean",
                           return_args = T)

# Compute petra_edr using method = "medoid"
petra_method2 <- petra_edr(state_var = state_var1, 
                           trajectories = trajectories1, 
                           states = states1, 
                           targets = "target1",
                           d_function = "vegan::vegdist", 
                           d_args = list(x = state_var1, method = "bray"),
                           k = 20L, 
                           minPts = 2L, 
                           method = "medoid",
                           return_args = T)

## ----method1, fig.width=8, fig.height=4, warning=FALSE------------------------
par(mfrow = c(1, 2))
plot(petra_method1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlim = c(-0.1, 0.3), 
     ylim = c(0.1, 0.4),
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "method = mean")

plot(petra_method2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlim = c(-0.1, 0.3), 
     ylim = c(0.1, 0.4),
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "method = medoid")

## ----direction, fig.width=9, fig.height=3, warning=FALSE----------------------
# Compute petra_edr before the target
petra_direction1 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = -1, 
                              return_args = T)

  # Compute petra_edr after the target
petra_direction2 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = 1, 
                              return_args = T)

  # Compute petra_edr before and after the target
petra_direction3 <- petra_edr(state_var = state_var1, 
                              trajectories = trajectories1, 
                              states = states1, 
                              targets = "target1",
                              d_function = "vegan::vegdist", 
                              d_args = list(x = state_var1, method="bray"),
                              k = 20L, 
                              minPts = 2L, 
                              direction = 2,
                              return_args = T)

  # Plot PETRA outputs
par(mfrow = c(1, 3))
plot(x = petra_direction1, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = -1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_direction2, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = 1")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = c(NA, 2, 2), 
       pch = c(20, NA, NA),
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

plot(x = petra_direction3, 
     target.colors = "red", 
     petra.colors = "blue", 
     xlab = "MDS D1", 
     ylab = "MDS D2", 
     main = "direction = 2")
legend("bottomleft", 
       c("Target", "Predicted trajectory", "EDR trajectories"), 
       lwd = 2,
       col = c("red", "blue", "grey"), 
       cex = 0.8, 
       bty = "n")

## ----state_var2---------------------------------------------------------------
petra_k2$state_var

## ----trajectories2------------------------------------------------------------
petra_k2$trajectories

## ----states2------------------------------------------------------------------
petra_k2$states

## ----state_var trajectories states--------------------------------------------
state_var1_k2 <- petra_k2$state_var
state_var1_k2$traj <- petra_k2$trajectories
state_var1_k2$state <- petra_k2$states
state_var1_k2[, c("traj", "state", paste0("sp", 1:12))]

## ----k_dist-------------------------------------------------------------------
petra_eps1$k_dist

## ----predicted_dist-----------------------------------------------------------
petra_eps1$predicted_dist

## ----dStates------------------------------------------------------------------
# State variables and dissimilarity metric of the EDR and the first target
state_var1 <- rbind(edr[, paste0("sp", 1:12)], 
                    target1[, paste0("sp", 1:12)])
dStates1 <- vegan::vegdist(state_var1, method = "bray")

# State variables and dissimilarity metric of the EDR and the second target
state_var2 <- rbind(edr[, paste0("sp", 1:12)], 
                    target2[, paste0("sp", 1:12)])
dStates2 <- vegan::vegdist(state_var2, method = "bray")

# State variables and dissimilarity metric of the EDR and the third target
state_var3 <- rbind(edr[, paste0("sp", 1:12)], 
                    target3[, paste0("sp", 1:12)])
dStates3 <- vegan::vegdist(state_var3, method = "bray")


## ----plot_edr, fig.width=9, fig.height=3, warning=FALSE-----------------------
par(mfrow = c(1, 3))

# Number of trajectories in the EDR
Ntraj <- length(unique(edr$traj))

# Location of target1. 
# As target1 is composed of one state, we need to specify type = 'states'
plot_edr(x = dStates1, 
         trajectories = c(edr$traj, target1$traj),
         states = as.integer(c(edr$state, target1$state)),
         type = "states",
         state.colors = c(rep("grey", length(edr$traj)), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the first target")
legend("bottomleft", 
       c("EDR trajectories", "Target 1"), 
       lwd = c(2, NA), 
       pch = c(NA, 20), 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")


# Location of target2
plot_edr(x = dStates2, 
         trajectories = c(edr$traj, target2$traj),
         states = as.integer(c(edr$state, target2$state)),
         traj.colors = c(rep("grey", Ntraj), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the second target")
legend("bottomleft", 
       c("EDR trajectories", "Target 2"), 
       lwd = 2, 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")

# Location of target3
plot_edr(x = dStates3, 
         trajectories = c(edr$traj, target3$traj),
         states = as.integer(c(edr$state, target3$state)),
         traj.colors = c(rep("grey", Ntraj), "red"),
         xlab = "MDS D1", 
         ylab = "MDS D2", 
         main = "Location of the third target")
legend("bottomleft", 
       c("EDR trajectories", "Target 3"), 
       lwd = 2, 
       col = c("grey", "red"), 
       cex = 0.8, 
       bty = "n")

## ----dDis---------------------------------------------------------------------
# We cannnot compute dDis for the first target because it is composed of a single state. We will calculate an equivalent metric by calculating the average dissimilarity between the target and the EDR states.
dDis1 <- mean(as.matrix(dStates1)[-nrow(state_var1), nrow(state_var1)])
names(dDis1) <- "dDis (ref. target1)"

# For targets 2 and 3, we use the function dDis
dDis2 <- dDis(d = dStates2, 
              d.type = 'dStates', 
              trajectories = c(edr$traj, target2$traj),
              states = c(edr$state, target2$state), 
              reference = 'target2')
dDis3 <- dDis(d = dStates3, 
              d.type = 'dStates', 
              trajectories = c(edr$traj, target3$traj),
              states = c(edr$state, target3$state), 
              reference = 'target3')

dDis1; dDis2; dDis3


## ----state_var3---------------------------------------------------------------
# data.table including the state variables, trajectories, and states of the EDR and the targets
data <- rbind(edr[, -1], target1, target2, target3)

# state_var needs to be a data.frame with only the state variables
state_var <- data.frame(data[, paste0("sp", 1:12)])

## ----target_pars, echo=FALSE--------------------------------------------------
target_pars <- data.frame(Target = paste0("target", 1:3),
                          k = c(10L, 50L, 50L),
                          eps = c(NA, 0.05, 0.5), 
                          minPts = c(3L, 2L, 2L),
                          w_function = c(NA, "exponential", "linear"), 
                          alpha = c(NA, 3, NA))

knitr::kable(target_pars, row.names = F)

## ----petra--------------------------------------------------------------------
petra <- petra_edr(state_var = state_var,
                   trajectories = data$traj,
                   states = as.integer(data$state),
                   targets = c("target1", "target2", "target3"),
                   d_function = "vegan::vegdist", 
                   d_args = list(x = state_var, method = "bray"),
                   k = c(10L, 50L, 50L), 
                   eps = c(NA, 0.05, 0.5), 
                   minPts = c(3L, 2L, 2L),
                   w_function = c(NA, "exponential", "linear"), 
                   alpha = c(NA, 3, NA),
                   direction = 2, 
                   method = "mean", 
                   return_args = T)

## ----petra k_dist 1-----------------------------------------------------------
petra$k_dist[target == "target1"]

## ----petra k_dist 2-----------------------------------------------------------
petra$k_dist[target == "target2"]

## ----petra k_dist 3-----------------------------------------------------------
petra$k_dist[target == "target3"]

## ----predicted_dist 1---------------------------------------------------------
petra$predicted_dist[target == "target1"]

## ----predicted_dist 2---------------------------------------------------------
petra$predicted_dist[target == "target2"]

## ----predicted_dist 3---------------------------------------------------------
petra$predicted_dist[target == "target3"]

## ----MPD----------------------------------------------------------------------
MPD(x = petra)

## ----plot predicted, fig.width=5, fig.height=5, warning=FALSE-----------------
plot(x = petra,
     xlab = "MDS D1", 
     ylab = "MDS D2")
legend("bottomleft", 
       c("Predicted trajectories", "EDR trajectories"),
       lwd = 2, 
       col = c("red", "grey"), 
       cex = 0.8, 
       bty = "n")

## ----plot predicted col, fig.width=5, fig.height=5, warning=FALSE-------------
plot(x = petra, 
     petra.colors = grDevices::palette.colors(6, "Paired")[c(2, 4, 6)],
     xlab = "MDS D1", 
     ylab = "MDS D2")
legend("topleft", 
       c("Predicted trajectory 1", "Predicted trajectory 2",
         "Predicted trajectory 3", "EDR trajectories"),
       lwd = 2, 
       col = c(grDevices::palette.colors(6, "Paired")[c(2, 4, 6)], "grey"), 
       cex = 0.8, 
       bty = "n")

## ----plot targets, fig.width=5, fig.height=5, warning=FALSE-------------------
plot(x = petra, 
     traj.colors = grDevices::palette.colors(9, "Set 3")[9],
     petra.colors = grDevices::palette.colors(6, "Paired")[c(1, 3, 5)],
     target.colors = grDevices::palette.colors(6, "Paired")[c(2, 4, 6)],
     xlab = "MDS D1", 
     ylab = "MDS D2")

legend("bottomleft", 
       c("Target 1", "Predicted trajectory 1", 
         "Target 2", "Predicted trajectory 2",
         "Target 3", "Predicted trajectory 3", 
         "EDR trajectories"),
       col = c(grDevices::palette.colors(6, "Paired")[c(2, 1, 4, 3, 6, 5)],
               grDevices::palette.colors(9, "Paired")[9]),
       lwd = 2, 
       ncol = 2, 
       cex = 0.8, 
       bty = "n")

## ----plot mean_dist, fig.width=5, fig.height=5, warning=FALSE-----------------
plot(x = petra, 
     petra.colors = grDevices::hcl.colors(5, "Viridis")[3],
     target.colors = grDevices::hcl.colors(5, "Viridis")[1], 
     uncert.metric = "mean_dist", 
     uncert.colors = grDevices::hcl.colors(5, "Viridis"),
     xlab = "MDS D1", 
     ylab = "MDS D2")

legend("topleft", 
       legend = c(paste0("mean_dist = ", 
                         round(min(petra$predicted_dist$mean_dist), 2)),
                  rep(NA, 18), 
                  paste0("mean_dist = ",
                         round(max(petra$predicted_dist$mean_dist), 2))),
       fill = grDevices::hcl.colors(20, "Viridis"), 
       border = NA, 
       y.intersp = 0.2,
       cex = 0.8, 
       bty = "n")