knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(cache.extra = knitr::rand_seed)
setwd("~/GitHub/SWF-molder/R")
library(terra)
library(devtools)
library(parallel)
library(igraph)
load_all()
## ℹ Loading SWFmolder
## Warning: Objects listed as exports, but not present in namespace:
## • ClosestDiagonalAgriCell
extractRandWindow <- function(x, p){
  firstIndex = sample(seq(length(x) - p + 1), 1)
  x[firstIndex:(firstIndex + p -1)]}
  color_map <- c('#ef8a62', '#f7f7f7', '#67a9cf')  
dims <- seq(10,170,20)
iters = 50
swfCat = 1
agriCat = 2
Q = 2 # 2 px allocated per time
NNeighbors = 3 # 3 NN
swfCover = 0.9 
np = detectCores()

1 5 % initial SWF cover

Eben.res <- lapply(dims, function(dim) {
    message("Dimension: ", dim)
    Sys.sleep(1)
    null.mt <- matrix(2, ncol=dim,nrow=dim)
    set.seed(123)
    # 5 % inititial SWF cover
    null.mt[sample((1:dim^2),floor(dim^2*0.05))] <-1
    # 5 % inititial Forest cover
    null.mt[as.vector(replicate(3,extractRandWindow(1:dim^2,floor(dim^2*0.05))))] <-3
    # plot(rast(t(null.mt)),col=color_map,cex=10)
    null.gr <- GfM(null.mt)

    matrix.time <- system.time(matrix.m <- swf.molder(Hmatrix=null.mt, swfCover=swfCover, swfCat=swfCat, agriCat=agriCat, Q=Q, ExpPriority="mixed", ExpDirection="mixed", reduceQTo=0, iterations = iters, NNeighbors=NNeighbors, maxDistance = 1, queensCase=TRUE, maxGDistance=1, np=np))

    Sys.sleep(1)
    
    graph.time <- system.time(graph.m <- swf_molderN(g=null.gr, swfCat, agriCat, iterations = iters, Q, NNeighbors, swfCover, maxGDistance=1, max_radius=1, np=np))

    matrix.m1 <- lapply(graph.m, MfG, nrow(null.mt), ncol(null.mt))
    check <- identical(t(matrix.m[[length(matrix.m)]]), matrix.m1[[length(matrix.m1)]])
    return(list(do.call(rbind,list(graph.time,matrix.time)),check,matrix.m,matrix.m1))
    })
timeG <- sapply(1:length(Eben.res), function(x) Eben.res[[x]][[1]][1,3])
timeM <- sapply(1:length(Eben.res), function(x) Eben.res[[x]][[1]][2,3])
# sapply(1:length(ben.results), function(x) ben.results[[x]][[2]])
timePlot <- cbind(timeG,timeM,dims)

1.1 Benchmark execution times

dims2 <- timePlot[,3]^2
plot(timeM~dims2,timePlot,type="l",col="red",ylim=c(0,max(timeM,timeG)), ylab="Seconds elapsed", xlab="matrix dimension (pixels)")
lines(timeG~dims2,timePlot,type="l",col="blue")
legend("topleft",lty=c(1,1),col=c("red","blue"),legend=c("matrix","graph"))
title("5% initial SWF cover")

1.2 Overview of matrices/networks used for testing

ncol=4
par(mfrow = c(ceiling(length(dims)/ncol),ncol) , mai = c(0.1, 0.1, 0.1, 0.1))
for (dim in dims) {
    null.mt <- matrix(2, ncol=dim,nrow=dim)
    set.seed(123)
    null.mt[sample((1:dim^2),floor(dim^2*0.15))] <-1
    null.mt[as.vector(replicate(3,extractRandWindow(1:dim^2,floor(dim^2*0.05))))] <-3
    smfTemp <- rast(t(null.mt))
    plot(smfTemp, col = color_map, legend=FALSE, axes = FALSE, box = FALSE, main = paste("Dimension ", dim^2," pixels"))
    # text(smfTemp, digits=1)
}

2 15 % initial SWF cover

dims <- seq(10,170,20)
iters = 50
swfCat = 1
agriCat = 2
Q = 1 # 2 px allocated per time
NNeighbors = 3 # 3 NN
swfCover = 0.9 
np = detectCores()
Hben.res <- lapply(dims, function(dim) {
    message("Dimension: ", dim)
    null.mt <- matrix(2, ncol=dim,nrow=dim)
    set.seed(123)
    # 15 % inititial SWF cover
    null.mt[sample((1:dim^2),floor(dim^2*0.15))] <-1
    # 5 % inititial SWF cover
    null.mt[as.vector(replicate(3,extractRandWindow(1:dim^2,floor(dim^2*0.05))))] <-3
    # plot(rast(t(null.mt)),col=color_map,cex=10)
    null.gr <- GfM(null.mt)
    Sys.sleep(0.5)
    
    matrix.time <- system.time(matrix.m <- swf.molder(Hmatrix=null.mt, swfCover=swfCover, swfCat=swfCat, agriCat=agriCat, Q=Q, ExpPriority="mixed", ExpDirection="mixed", reduceQTo=0, iterations = iters, NNeighbors=NNeighbors, maxDistance = 1, queensCase=TRUE, maxGDistance=1, np=np))
    
    Sys.sleep(1)
    
    graph.time <- system.time(graph.m <- swf_molderN(g=null.gr, swfCat, agriCat, iterations = iters, Q, NNeighbors, swfCover, maxGDistance=1, max_radius=1, np=np))

    matrix.m1 <- lapply(graph.m, MfG, nrow(null.mt), ncol(null.mt))
    check <- identical(t(matrix.m[[length(matrix.m)]]), matrix.m1[[length(matrix.m1)]])

    return(list(do.call(rbind,list(graph.time,matrix.time)),check,matrix.m,matrix.m1))
    })
timeG <- sapply(1:length(Hben.res), function(x) Hben.res[[x]][[1]][1,3])
timeM <- sapply(1:length(Hben.res), function(x) Hben.res[[x]][[1]][2,3])
# sapply(1:length(Hben.res), function(x) Hben.res[[x]][[2]])
timePlot <- cbind(timeG,timeM,dims)

2.1 Benchmark execution times

dims2 <- timePlot[,3]^2
plot(timeM~dims2,timePlot,type="l",col="red",ylim=c(0,max(timeM,timeG)), ylab="Seconds elapsed", xlab="matrix dimension (pixels)")
lines(timeG~dims2,timePlot,type="l",col="blue")
legend("topleft",lty=c(1,1),col=c("blue","red"),legend=c("graph","matrix"))
title("15% initial SWF cover")

2.2 Overview of matrices/networks used for testing

ncol=4
par(mfrow = c(ceiling(length(dims)/ncol),ncol) , mai = c(0.1, 0.1, 0.1, 0.1))
for (dim in dims) {
    null.mt <- matrix(2, ncol=dim,nrow=dim)
    set.seed(123)
    null.mt[sample((1:dim^2),floor(dim^2*0.15))] <-1
    null.mt[as.vector(replicate(3,extractRandWindow(1:dim^2,floor(dim^2*0.15))))] <-3
    smfTemp <- rast(t(null.mt))
    plot(smfTemp, col = color_map, legend=FALSE, axes = FALSE, box = FALSE, main = paste("Dimension ", dim^2," pixels"))
    # text(smfTemp, digits=1)
}