Implement changes as requested by JNCC

Parses.R

@@ -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("Saving model prior to network modelling")
-  print(paste0(inputText, edges))
+  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
@@ -339,12 +353,6 @@ parseSheet <- function(fName) {
         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(
           list(

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")
+layers <- c("Pressures to Functional Groups", "Functional Groups to Output Processes", "Output Processes to Ecosystem services")
-transitions <- c("Pressures to Bio-Assemblages", "Pressures to Output Processes", "Pressures 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)
@@ -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)
 
     }
   }

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")
+}
+
+
+

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)
+
+}