diff --git a/MBA_MESO_Nodes.xlsx b/MBA_MESO_Nodes.xlsx new file mode 100644 index 0000000..5567bc8 Binary files /dev/null and b/MBA_MESO_Nodes.xlsx differ diff --git a/Parses.R b/Parses.R index d7a98c9..dd8172c 100644 --- a/Parses.R +++ b/Parses.R @@ -85,7 +85,7 @@ getInitial <- function(string, letter) { } split <- function(cell) { - + params <- unlist(strsplit(cell, ",")) values <- rep(0, length(states)) @@ -126,6 +126,7 @@ buildGraph <- function(model, desc) { #inputCode - the top layer of the model #outputCodes - all subsequent layers to be included in the model + inputNodes <- model$nodes$code[which(startsWith(model$nodes$code, desc$inputCode))] inputText <- paste0("[", inputNodes, "]", collapse = "") @@ -151,8 +152,15 @@ buildGraph <- function(model, desc) { outDist[[idx]] <- list(coef = coefVal, sd = model$nodes$confidence[nodeRef]) } - print("about to build network") - print(paste0(inputText, edges)) + print("Saving model prior to network modelling") + modelDefn <- paste0(inputText, edges) + save(modelDefn, file="buildGraph.RData") + + + #print("about to build network") + #print(paste0(inputText, edges)) + + net <- model2network(paste0(inputText, edges), debug = FALSE) @@ -167,6 +175,8 @@ buildGraph <- function(model, desc) { } allDists <- as.list(setNames(c(inDist, outDist), c(inputNodes, outNodes))) + + #print(allDists) cfit <- custom.fit(net, allDists) cat("about to calculate sample distributions") @@ -264,6 +274,8 @@ getCode <- function(name, nodeDF) { getValidEdges <- function(mapping, nodeDF, prevEdge = NULL, prefix) { #utils::str(nodeDF) + #save(mapping, nodeDF, prevEdge, prefix, file="validEdges.RData") + edgeCols <- c("inputNode", "outputNode", "impact") edgeM <- matrix(data = NA, nrow = 0, ncol = length(edgeCols), dimnames = list(NULL, edgeCols)) @@ -309,6 +321,8 @@ parseMapping <- function(mapping, prevOutputs, prefix) { nodeDF <- getValidNodes(mapping, prevOutputs$nodes, prefix) edgeDF <- getValidEdges(mapping, nodeDF, prevEdge = prevOutputs$edges, prefix) + #save(nodeDF, edgeDF, file="mapping.RData") + return(list( #New structure nodes = nodeDF, @@ -329,7 +343,7 @@ parseSheet <- function(fName) { sheets <- sort(delNA(match(names, mappings))) cat("starting sheet parse") - print(sheets) + #print(sheets) if (sum(sheets == refs) == length(refs)) { #read all mapping tables @@ -338,12 +352,6 @@ parseSheet <- function(fName) { p_op <- parseMapping(readXL(fName,mappings[3], startRow = 1), p_ba, prefix = "op") p_es <- parseMapping(readXL(fName,mappings[4], startRow = 1), p_op, prefix = "es") legend <- readXL(fName,mappings[5], startRow = 1) - - #print("building graphs") - - #p_baNet <- buildGraph(p_ba, desc = list(inputCode = "p", outputCodes = "ba")) - #p_opNet <- buildGraph(p_op, desc = list(inputCode = "p", outputCodes = c("ba", "op"))) - #p_esNet <- buildGraph(p_es, desc = list(inputCode = "p", outputCodes = c("ba", "op", "es"))) print("sheet load completed") return( diff --git a/app.R b/app.R index 9578432..d866c99 100644 --- a/app.R +++ b/app.R @@ -3,6 +3,7 @@ modules::import(shinydashboard) modules::import(shinydashboardPlus) modules::import(shinycssloaders) modules::import(shinyjs) +modules::import(shinyBS) modules::import(bnlearn) modules::import(visNetwork) @@ -12,21 +13,25 @@ modules::import(openxlsx) modules::import(zip) modules::import(DT) modules::import(plyr) +modules::import(magrittr) parser <- modules::use("Parses.R") +rw <- modules::use("reWeight.R") + addResourcePath("js", "./www/js") -layers <- c("Pressures to Bio-Assemblages", "Bio-Assemblages to Output Processes", "Output Processes to Ecosystem services") -transitions <- c("Pressures to Bio-Assemblages", "Pressures to Output Processes", "Pressures to Ecosystem services") +layers <- c("Pressures to Functional Groups", "Functional Groups to Output Processes", "Output Processes to Ecosystem services") +transitions <- c("Pressures to Functional Groups", "Pressures to Output Processes", "Pressures to Ecosystem services") impacts <- c("Very High", ">= High", ">= Medium", ">= Low", "All") thresholds <- c(0.97, 0.9, 0.45, 0.17, 0) impLabels <- c("Very High", "High", "Medium", "Low", "Very Low") ui <- dashboardPage( + dashboardHeader(title = "JNCC MESO online", tags$li( id = "dropdownHelp", @@ -82,6 +87,7 @@ ui <- dashboardPage( #menuItem("Habitats", tabName = "3", icon = icon("atlas")), #menuItem("Ingestion", tabName = "3", icon = icon("utensils")), selectInput("modelSelect", "Select MESO model", choices = c(""), selected = NULL, multiple = FALSE), + #downloadButton("download", "", icon = icon("download")), uiOutput("pressureList") ) @@ -143,7 +149,7 @@ ui <- dashboardPage( fluidRow( column( width = 6, - h4("Effect on bio-assemblage") + h4("Effect on Functional Groups") ), column( width = 1, @@ -155,7 +161,8 @@ ui <- dashboardPage( ), column( width = 1, - downloadButton("download", "", icon = icon("download")) + downloadButton("download", "", icon = icon("download")), + shinyBS::bsTooltip("download", "Template provides for decimal values in degs column OR degs:mins:secs. Longitude west of meridian must be negative.") ), column( width = 2, @@ -174,11 +181,13 @@ ui <- dashboardPage( fluidRow( column( width = 4, - checkboxInput("bbnDisplayNames", "Display Node names", value = FALSE) + checkboxInput("bbnDisplayNames", "Display Node names", value = FALSE), + shinyBS::bsTooltip("bbnDisplayNames", "Four MESO models have been defined thus far") ), column( width = 4, - checkboxInput("bbnDisplayEdges", "Display edge status", value = FALSE) + checkboxInput("bbnDisplayEdges", "Display edge status", value = FALSE), + shinyBS::bsTooltip("bbnDisplayEdges", "Edges are removed") ), column( width = 4, @@ -261,6 +270,42 @@ server <- function(input, output, session) { as.numeric(v) } + newNameMap <- openxlsx::read.xlsx("MBA_MESO_Nodes.xlsx") %>% + dplyr::select(hab, nodeType, Suggestion, node, newname) + #save(newNameMap, file="nameMap.RData") + + stripStr <- function(nodeStr) { + nodeStr %>% stringr::str_replace_all("\\.", "") %>% + stringr::str_replace_all(" ", "") %>% + stringr::str_replace_all("\\(", "") %>% + stringr::str_replace_all("\\)", "") %>% + stringr::str_replace_all("\\/", "") %>% + tolower() + } + + setNewNames <- function(wb, habName) { + + #habName <- substr(fileList[idx], 1, (nchar(fileList[idx])-5)) + + print(habName) + possNames <- newNameMap %>% + dplyr::filter(hab==habName) %>% + dplyr::mutate(node=stripStr(node)) + + newNodes <- wb$p_es$nodes %>% dplyr::mutate(node=stripStr(name)) + + print(possNames$node) + print(newNodes$node) + newNames <- apply(newNodes, 1, function(row) { + id <- match(row["node"], possNames$node) + print(paste(id, row["node"])) + possNames$newname[id] + }) + print(newNames) + wb$p_es$nodes$name <- newNames + return(wb) + } + getAvailableModels <- function() { fileList <- list.files(dataStorage, pattern = ".xlsx") @@ -276,13 +321,20 @@ server <- function(input, output, session) { wb$p_es$edges$values <- sapply(wb$p_es$edges$impact, getImpact) if (!is.null(wb)) { - modelList[[cnt]] <- wb - models <<- c(models, substr(fileList[idx], 1, (nchar(fileList[idx])-5))) + + habName <- substr(fileList[idx], 1, (nchar(fileList[idx])-5)) + + wb2 <- setNewNames(wb, habName) + + modelList[[cnt]] <- wb2 + models <<- c(models, habName) print(paste("Model file successfully loaded", fileList[idx])) + #save(tmp, file = "tmp.RData") cnt <- cnt+1 } } + #save(modelList, file="models.RData") updateSelectInput(session, "modelSelect", choices = models) return(modelList) } @@ -290,6 +342,10 @@ server <- function(input, output, session) { #parse on load sheets in the input sheet folder - replace with R Data modelList <- getAvailableModels() + save(modelList, file="model.RData") + + #print(load("modelList.RData")) + calcLikelihood <- function(layer, pressStatus, forPlotly) { @@ -301,7 +357,6 @@ server <- function(input, output, session) { thisModel <- modelList[[.selections$model]] - MEANPOS <- 1 MEANNEG <- 0 @@ -318,6 +373,25 @@ server <- function(input, output, session) { expr <- substr(expr, 1, nchar(expr)-2) expr <- paste0(expr, ")") + print(names(thisModel)) + + #Now do it in stages with one assessment per stage + + + + thisModel$p_es$nodes$confidence <- 0.1 * thisModel$p_es$nodes$confidence + + + #save(pressStatus, thisModel, file="beforeWeight.RData") + + + + if (sum(pressStatus$status=="On")>0) { + thisModel$p_es <- rw$reWeightModel(thisModel$p_es, pressStatus) + } #else nothing to do + + #save(pressStatus, thisModel, file="afterWeight.RData") + thisNet <- parser$buildGraph(thisModel$p_es, desc = list(inputCode = "p", outputCodes = c("ba", "op", "es"))) sampleDists <- cpdist( @@ -333,7 +407,7 @@ server <- function(input, output, session) { #print(sampleDists) #displayCols <- match(nodeCodes, colnames(sampleDists)) - sampleDists <- sampleDists[,match(thisModel$p_es$nodes$code, colnames(sampleDists))] + sampleDists <- round(sampleDists[,match(thisModel$p_es$nodes$code, colnames(sampleDists))], digits=2) means <- apply(sampleDists, 2, mean) stdDev <- apply(sampleDists, 2, sd) @@ -341,7 +415,7 @@ server <- function(input, output, session) { quantiles <- t(apply(sampleDists, 2, quantile, c(0.01, 0.25, 0.5, 0.75, 0.99))) print(paste("Building likelihoods from model, sample dists", length(thisModel$p_es$nodes$name), length(sampleDists))) #str(quantiles) - + if (forPlotly) { return(data.frame( name = thisModel$p_es$nodes$name, @@ -371,7 +445,7 @@ server <- function(input, output, session) { maxes = apply(sampleDists, 2, max), stringsAsFactors = FALSE )) - + } } @@ -397,7 +471,7 @@ server <- function(input, output, session) { #.selections$runOnce = FALSE print("Running calc") .likelihoods$p_es <<- calcLikelihood(0, newStatus, TRUE) - + .selections$pressStatus <<- newStatus } @@ -420,7 +494,7 @@ server <- function(input, output, session) { #status = status, stringsAsFactors = FALSE ) - + #This assumes all pressures are the same... setPressures(pressures) @@ -466,7 +540,7 @@ server <- function(input, output, session) { }) observeEvent(input$modalOK, { - + .resistanceScores["nr"] <<- -input$l1VH .resistanceScores["lr"] <<- -input$l1H @@ -476,7 +550,7 @@ server <- function(input, output, session) { .resistanceScores["ssgr"] <<- input$ssgr .resistanceScores["pressSD"] <<- input$l1PressSD - + .likelihoods$p_es <<- calcLikelihood(0, .selections$pressStatus, TRUE) removeModal() @@ -570,7 +644,7 @@ server <- function(input, output, session) { } else { edgeNet <- edges } - + print(paste(nrow(model$legend), length(palette))) legendDF <- data.frame( @@ -619,7 +693,7 @@ server <- function(input, output, session) { zerolinewidth = 10) # plot_ly(boxPlot, x = boxPlot[,1], y = ~Range, color = as.character(boxPlot$Group), colors = palette, type = "box") %>% - layout(xaxis = xform, showlegend = FALSE, title = title) + layout(xaxis = xform, yaxis=list(range=c(-1.2, 1.2)), showlegend = FALSE, title = title) } } @@ -668,10 +742,10 @@ server <- function(input, output, session) { }, contentType = "application/xlsx" ) - + makeLikelihoods <- function() { - - + + likeliTab <- as.data.frame( cbind( .likelihoods$p_es, codeVal = sapply( @@ -682,15 +756,15 @@ server <- function(input, output, session) { )), stringsAsFactors=FALSE ) - + likeliTab <- arrange(likeliTab, layer, codeVal) - + outputRows <- trunc(nrow(likeliTab)/7) outputTab <- NULL - + for (idx in 1:outputRows) { elementRow <- (idx - 1) * 7 + 1 - + tabRow <-c( name = likeliTab$name[elementRow], code = likeliTab$code[elementRow], @@ -702,9 +776,9 @@ server <- function(input, output, session) { max =likeliTab$range[elementRow+6] ) outputTab <- rbind(outputTab, tabRow) - + } - + likelihoods <- data.frame( name = outputTab[,1], code = outputTab[,2], @@ -720,10 +794,10 @@ server <- function(input, output, session) { } output$download <- downloadHandler( - + filename = function() { paste0("MESO-", format(Sys.time(), "%m%d_%H%M"), ".xlsx") }, content = function(file) { - + showModal( modalDialog( fluidRow( @@ -734,13 +808,13 @@ server <- function(input, output, session) { ) oldDir <- getwd() - + tmp <- tempfile("") dir.create(tmp) setwd(tmp) - - - + + + l <- list( pressures = .selections$pressStatus, nodes = modelList[[.selections$model]]$p_es$nodes, @@ -751,7 +825,7 @@ server <- function(input, output, session) { xl <- write.xlsx(l, "dataset.xlsx") #zipFile <- zipr(file, c("dataset.xlsx")) - + file.copy("dataset.xlsx", file) #print(paste("zip file complete", zipFile)) diff --git a/extract.R b/extract.R new file mode 100644 index 0000000..7515eeb --- /dev/null +++ b/extract.R @@ -0,0 +1,111 @@ +#R script to upload the existing spreadsheets and homologise them +library(magrittr) +fList <- list.files("data", pattern="*.xlsx") + +#Objective to create data tables with +linkCheck <- function(nodeType, nodeString, nodeStringCheck) { + nodeString <- stringr::str_replace_all(nodeString, "\\.", " ") + res <- sapply(nodeString, match, nodeStringCheck$Nodes) %>% is.na() %>% which() + if (length(res)>0) print(paste("Clean up error found in", nodeType, "mapping at", names(res))) +} + +getNodeVals <- function(nodeStr) { + params <- stringr::str_split(nodeStr, ",") %>% unlist() %>% trimws() + paramVals <- stringr::str_split(params, "=") + vals <- c() + lapply(paramVals, function(l) { + val <- l[2] + names(val) <- l[1] + vals <<- c(vals, val) + }) + vals +} + +#We want to build a node table and an impact table. +#Colnames of the node table will be +#Hab, Node Type, Node, Node Layer, Growth, .... + +#The edges table will be +#Hab, In Node, Out Node, Params, .... + + +sheetNames <- c("TestScenario", "Map_P_BA", "Map_BA_OP", "Map_OP_ES", "Legend") + +cleanNames <- function(namVec) { + stringr::str_replace_all(namVec, "\\.", " ") %>% trimws() %>% tolower() +} + +nodeTable <- tibble::tibble() + +for (wbIdx in 1:length(fList)) { + wb <- openxlsx::loadWorkbook(paste0("data/", fList[wbIdx])) + hab <- stringr::str_split(fList[wbIdx], "\\.")[[1]][1] + #get pressure names + + #Drop the time column no use at all.... + sheet <- openxlsx::readWorkbook(wb, sheet=sheetNames[1])[ ,-1] + pressures <- cleanNames(colnames(sheet)) + pressure_nodes <- sheet[1,] + + + sheet <- openxlsx::readWorkbook(wb, sheet=sheetNames[2])[ ,-1] + pressure_check <- na.omit(sheet[,1:2]) + sheet2 <- na.omit(sheet[, -c(1,2)]) + ba <- cleanNames(colnames(sheet2)) + ba_nodes <- sheet2[1,] + pressImpact <- sheet2[-1,] + + #linkCheck("pressures", pressures, pressure_check) + + + sheet <- openxlsx::readWorkbook(wb, sheet=sheetNames[3])[ ,-1] + ba_check <- na.omit(sheet[,1:2]) + sheet2 <- na.omit(sheet[, -c(1,2)]) + op <- cleanNames(colnames(sheet2)) + op_nodes <- sheet2[1,] + baImpact <- sheet2[-1,] + + #linkCheck("bioassemblages", ba, ba_check) + + sheet <- openxlsx::readWorkbook(wb, sheet=sheetNames[4])[ ,-1] + op_check <- na.omit(sheet[,1:2]) + sheet2 <- na.omit(sheet[, -c(1,2)]) + es <- cleanNames(colnames(sheet2)) + es_nodes <- sheet2[1,] + opImpact <- sheet2[-1,] + + #linkCheck("outputprocesses", op, op_check) + + legend <- openxlsx::readWorkbook(wb, sheet=sheetNames[5]) + + nodeType <- c( + rep("pressure", length(pressures)), + rep("bioassemblage", length(ba)), + rep("outputprocess", length(op)), + rep("ecosystemservice", length(es)) + ) + + + + res <- t(sapply(es_nodes[1,], getNodeVals)) %>% as.data.frame() + names(res) <- cleanNames(names(res)) + res <- res %>% mutate(nodeName=names(res)) + + nodeTable <- nodeTable %>% dplyr::bind_rows( + tibble::tibble( + hab=hab, + nodeType=nodeType, + res + ) + ) + +} + +mapNewNames <- function() { + newNameMap <- openxlsx::read.xlsx("MBA_MESO_Nodes.xlsx") %>% + dplyr::select(hab, nodeType, Suggestion, node, newname) + save(newNameMap, file="nameMap.RData") +} + + + diff --git a/reWeight.R b/reWeight.R new file mode 100644 index 0000000..b4bd9c7 --- /dev/null +++ b/reWeight.R @@ -0,0 +1,132 @@ +modules::import(magrittr) + +reWeightLayer <- function(nestedLayerTib, fudge=1) { + + for (idx in 1:nrow(nestedLayerTib)) { + #print(nestedLayerTib$data[idx]) + thisData <- nestedLayerTib$data[idx][[1]] + + #Calculate the overall depletion rate + #depRate <- ifelse(thisData$values<0, -thisData$values, 0) + #Re-adjust those weightings in line with the number applied + survived <- 1 + grown <- 1 + for (depIdx in 1:nrow(thisData)) { + if (thisData$values[depIdx]<0) survived <- survived * (1 + thisData$values[depIdx]) else + grown <- (1-thisData$values[depIdx]) * grown + } + #Update the edge weightings to reflect the combined depletion on the BA from each of the edges + + effDepRate <- survived - 1 + effGrowthRate <- 1-grown + #print(effDepRate) + if (sum(thisData$values)==0) newValues <- rep(0, length(thisData$values)) else + newValues <- round(thisData$values/sum(thisData$values)*(effDepRate+effGrowthRate), digits=3) + #print(paste(idx, paste(newValues, collapse=","))) + nestedLayerTib$data[idx][[1]]$values <- newValues / fudge + } + + return(nestedLayerTib %>% tidyr::unnest(cols=c(data))) +} + +assignWeights <- function( + edgesTib, + incode, + outcode, + value) { + for (idx in 1:length(incode)) { + ref <- intersect(which(edgesTib$input == incode[idx]), + which(edgesTib$output == outcode[idx])) + + utils::str(ref) + + if (length(ref)>1) stop("Error has occurred with multiple edges between two nodes") + print(paste(ref, edgesTib$values[ref], value[idx])) + edgesTib$values[ref] <- value[idx] + #Set the appropriate values + + } + return(edgesTib) +} + +reWeightModel <- function(thisNet, pressStatus) { + + print("About to recalc p - ba") + + #what is the depletion factor for each of the pressures applied to the BA? + p_on <- pressStatus %>% + dplyr::filter(status=="On") %>% + dplyr::left_join(thisNet$nodes, by=c("code"="code")) %>% + dplyr::left_join(thisNet$edges, by=c("code"="input")) %>% + dplyr::mutate(values=values * 0.9) + + print("before") + print(sum(p_on$values)) + + p_on <- p_on %>% + dplyr::rename(presscode=code) %>% + dplyr::rename(ba_code=output) %>% + dplyr::select(presscode, layer, ba_code, values) %>% + tidyr::nest(data=c(presscode, values)) + + newP <- reWeightLayer(p_on, fudge=1) + + + + print("About to recalc ba - op") + + #Repeat for the linkage between ba and op + bas <- unique(newP$ba_code) + ba_impacted <- thisNet$nodes %>% + dplyr::filter(code %in% bas) %>% + dplyr::left_join(thisNet$edges, by=c("code"="input")) %>% + tidyr::drop_na() %>% + dplyr::rename(ba_code=code) %>% + dplyr::select(layer, output, ba_code, values) %>% + dplyr::rename(op_code=output) %>% + tidyr::nest(data=c(ba_code, values)) + + newBA <- reWeightLayer(ba_impacted, fudge=4) + + print("About to recalc op - es") + + #Repeat for the linkage between op and es + ops <- unique(newBA$op_code) + op_impacted <- thisNet$nodes %>% + dplyr::filter(code %in% ops) %>% + dplyr::left_join(thisNet$edges, by=c("code"="input")) %>% + dplyr::rename(op_code=code) %>% + tidyr::drop_na() %>% + dplyr::select(layer, output, op_code, values) %>% + dplyr::rename(es_code=output) %>% + tidyr::nest(data=c(op_code, values)) + + newOP <- reWeightLayer(op_impacted, fudge=4) + + #Check for any more links through the system + print("About to recalc es - es") + + + ess <- unique(newOP$es_code) + es_impacted <- thisNet$nodes %>% + dplyr::filter(code %in% ess) %>% + dplyr::left_join(thisNet$edges, by=c("code"="input")) %>% + dplyr::rename(es_code=code) %>% + tidyr::drop_na() %>% + dplyr::select(layer, output, es_code, values) %>% + dplyr::rename(lo_code=output) %>% + tidyr::nest(data=c(lo_code, values)) + + newES <- reWeightLayer(es_impacted, fudge=2) + + incode <- c(newP$presscode, newBA$ba_code, newOP$op_code, newES$es_code) + outcode <- c(newP$ba_code, newBA$op_code, newOP$es_code, newES$lo_code) + value <- c(newP$values, newBA$values, newOP$values, newES$values) + + thisNet$edges <- assignWeights(thisNet$edges, incode, outcode, value) + + print("exitting reweighting process") + + return(thisNet) + +}