# R/irl_graph.R In jsta/irlgraph: Generate Accumulated Cost Surfaces with Irregular Landscape Graphs

#### Documented in irl_graph

```#'@name irl_graph
#'@title Generate an irregular landscape graph
#'@description Subset the nodes in a matrix of cost values based on user defined point-of-interest coordinates, a user-defined importance cutoff, and a randomly selected proportion of non-null cells. Generate a graph object from this subset using delaunay triangulation.
#'@param dm matrix of 2d cost values
#'@param poicoords matrix 2 column point-of-interest coordinates
#'@param grainprop numeric proportion of cells selected as grain cells
#'@param cutoff numeric user-defined importance cutoff
#'@param irregular logical skip irregular node subsetting routine?
#'@import sp
#'@importFrom Matrix Matrix
#'@importFrom raster raster extract boundaries cellFromCol cellFromRow cellFromXY xyFromCell
#'@importFrom geometry delaunayn
#'@export
#'@examples \dontrun{
#'
#'dm <- matrix(
#'c(1,  1,  1,  1,  1,
#'  NA, NA, NA, NA, 1,
#'  1,  1,  1,  1,  1,
#'  1,  1,  1,  1,  1,
#'  1,  1,  1,  1,  1
#'  ), ncol = 5, byrow = TRUE)
#'
#'sp::plot(raster::raster(dm))
#'raster::text(raster::raster(dm), which(raster::raster(dm)[] == 1))
#'
#'graph <- irl_graph(dm)
#'plot(graph\$graph)
#'geometry::trimesh(graph\$tri\$tri, graph\$allcoords)
#'}

irl_graph <- function(dm, poicoords = NA, grainprop = 0.25, cutoff = 0, irregular = TRUE){
csurf <- raster::raster(nrows = dim(dm)[1], ncols = dim(dm)[2], resolution = 1, xmn = 0, xmx = dim(dm)[1], ymn = 0, ymx = dim(dm)[2])
csurf[] <- dm

get_cells <- function(dm, grainprop, poicoords, cutoff = 0){

csurf <- raster::raster(nrows = dim(dm)[1], ncols = dim(dm)[2], resolution = 1, xmn = 0, xmx = dim(dm)[1], ymn = 0, ymx = dim(dm)[2])
csurf[] <- dm

xmax <- dim(csurf)[1]
ymax <- dim(csurf)[2]

allcells <- raster::rasterToPoints(is.na(csurf))
nullcells <- which(is.na(raster::extract(csurf,allcells[,1:2])))
limitcells <- c(raster::cellFromRow(csurf, c(1, ymax)), raster::cellFromCol(csurf, c(1, xmax)))
limitcells <- limitcells[!(limitcells %in% nullcells)]
limitcells <- c(limitcells,which(raster::extract(raster::boundaries(csurf),allcells[,1:2]) == 1))
limitcells <- limitcells[!duplicated(limitcells)]

#allow for user-defined "importance cutoff"
focalsurf <- raster::focal(csurf, w = matrix(1/9, nrow = 3, ncol = 3), sd)
vipcells <- which(focalsurf[] > cutoff)
vipcells <- vipcells[!(vipcells %in% limitcells)]
vipcells <- vipcells[!(vipcells %in% nullcells)]

uncells <- raster::cellFromXY(csurf, allcells)[!raster::cellFromXY(csurf, allcells) %in% c(nullcells, vipcells, limitcells)]

graincells <- sample(uncells, length(uncells) * grainprop)

poicells<-NULL
if(length(poicoords)>1){
poicells <- apply(poicoords, 1, function(x) raster::cellFromXY(csurf, x))
}

cells <- c(limitcells, vipcells, graincells, poicells)
cells <- cells[!duplicated(cells)]
cells <- cells[order(cells)]

return(list(cells = cells, nullcells = nullcells, limitcells = limitcells, vipcells = vipcells, graincells = graincells, poicells = poicells, allcells = which(allcells[,3]==0)))
}

cells <- get_cells(dm = dm, poicoords = poicoords, grainprop = grainprop, cutoff = cutoff)

#cells contains all the cells in the following block
poicells <- cells\$poicells
limitcells <- cells\$limitcells
vipcells <- cells\$vipcells
graincells <- cells\$graincells

nullcells <- cells\$nullcells
nullcoords <- raster::xyFromCell(csurf, nullcells)[,1:2]

if(irregular == TRUE){
cells <- cells\$cells
}else{
cells <- cells\$allcells
}
cellcoords <- raster::xyFromCell(csurf, cells)[,1:2]

allcells <- c(cells, nullcells)
allcells <- allcells[order(allcells)]
allcoords <- rbind(cellcoords, nullcoords)
allcoords <- allcoords[order(-allcoords[,2], allcoords[,1]),]

#create graph====================================================#
create_tri <- function(cellcoords){
#deldir::deldir(cellcoords[,1], cellcoords[,2], suppressMsge = TRUE)

tri <- geometry::delaunayn(cellcoords)
#trimesh(tri, allcoords)

nodes <- rbind(tri[,-1], tri[,-2], tri[,-3])
nodes <- nodes[order(nodes[,1]),]
list(nodes = nodes[!duplicated(paste(nodes[,1], nodes[,2])),], tri = tri)
}

dtri <- create_tri(allcoords)
#trimesh(dtri\$tri, allcoords)

orig <- allcells[dtri\$nodes[,1]]
neigh <- allcells[dtri\$nodes[,2]]

keep <- !(orig %in% nullcells) & !(neigh %in% nullcells)
orig <- orig[keep]
neigh <- neigh[keep]
elist <- cbind(orig, neigh)

edgelengths<-rep(1, nrow(elist))
celldiff <- elist[,1] - elist [,2]
longedges<-which(!celldiff %in% c(1, dim(csurf)[1]))
longedgecoords_from<-raster::xyFromCell(csurf, elist[longedges,1])
longedgecoords_to<-raster::xyFromCell(csurf, elist[longedges,2])

elength_from_coords <- function(from, to){
a <- from[,1] - to[,1]
b <- from[,2] - to[,2]
sqrt(a^2 + b^2)
}

edgelengths[longedges] <- elength_from_coords(longedgecoords_from, longedgecoords_to)

#edgelengths[longedges] <- apply(cbind(longedgecoords_from, longedgecoords_to), 1, function(x) elength_from_coords(x[1:2], x[3:4]))

#if(length(unique(elist[,1]))!=length(cells)){
#  warning("Mismatch between edgelist and nonnull cells")
#}

datavals <- cbind(raster::getValues(csurf)[elist[,1]], raster::getValues(csurf)[elist[,2]])
datavals <- apply(datavals, 1, function(x) mean(x, na.rm = TRUE))

adjmat <- Matrix::Matrix(0, nrow = raster::ncell(csurf), ncol =  raster::ncell(csurf))
adjmat[elist] <- as.vector(1 / (datavals * edgelengths))