cpdist integrated into model

Parses.R

@@ -30,6 +30,20 @@ delNA <- function(vec) {
   return(vec[!is.na(vec)])
 }
 
+buildExpr <- function(pressStatus) {
+  #pressStatus is a two column DF of name of pressure and status Ii.e. on or off)
+  MEANPRESS = 0
+  expr <- "("
+  for (p in 1:nrow(pressStatus)) {
+    if (pressStatus$status[p] == 'On') symbol='>=' else symbol='<='
+    expr <- paste0(expr, "(\"", pressStatus$code[p], "\"", symbol, MEANPRESS, ") & ")
+  }
+  expr<-substr(expr, 1, nchar(expr)-2)
+  expr<-paste0(expr, ')')
+  
+  return(expr)
+}
+
 parseScenario <- function(press, prefix = 'p') {
   pressNames <- colnames(press)[2:length(colnames(press))]
   coefs <- matrix(data=NA, nrow=length(pressNames), ncol=2, dimnames=list(NULL, c('growth', 'confidence')))

app.R

@@ -3,6 +3,7 @@ modules::import(shiny)
 modules::import(shinyBS)
 modules::import(shinyjs)
 modules::import(shinydashboard)
+modules::import(shinydashboardPlus)
 modules::import(htmltools)
 modules::import(DiagrammeR)
 modules::import(magrittr)
@@ -12,7 +13,6 @@ modules::import(Rgraphviz)
 modules::import(knitr)
 modules::import(shinycssloaders)
 modules::import(googleway)
-modules::import(Rgraphviz)
 modules::import(bnlearn)
 
 parser <- modules::use('Parses.R')
@@ -21,13 +21,21 @@ layers <- c("Pressures to Bio-Assemblages", "Bio-Assemblages to Output Processes
 transitions <- c("Pressures to Bio-Assemblages", "Pressures to Output Processes", "Pressures to Eco-system services")
 addResourcePath("js", "./www/js")
 
+
+    
+
 ui<-dashboardPage(
       dashboardHeader(title = "JNCC MESO online"),
+      
+      #tags$style(.times-circle {color:800000 }),
+      #tags$style(.check-square {color:008000 }),
+      
       dashboardSidebar(
         sidebarMenu(id = "tabs",
             menuItem("Pressure Test", tabName = "1", icon = icon("arrow-down")),
             menuItem("Bayesian Network", tabName = "2", icon = icon("atom")),
             menuItem("Habitats", tabName = "3", icon = icon("atlas")),
+            menuItem("Ingestion", tabName = "4", icon = icon("utensils")),
             selectInput("modelSelect", "Select MESO model", choices=c(""), selected=NULL, multiple=FALSE),
             selectInput("layerSelect", "Select Transition", 
                         choices=transitions, 
@@ -37,56 +45,54 @@ ui<-dashboardPage(
       dashboardBody(
         tabItems(
           tabItem(tabName = "1",
-                 fluidRow(
+             fluidRow(
-                  column(width=2,
+              column(width=2,
-                    h4('Pressure Test'),
+                h4('Pressure Test'),
-                    radioButtons("pressure1", 'Sediment type', choices=c('On', 'Off'), inline=TRUE),
+                uiOutput("pressureList"),
-                    radioButtons("pressure2", 'Seabed type', choices=c('On', 'Off'), inline=TRUE),
+                actionButton("calcAB", icon('calculator'))
-                    radioButtons("pressure3", 'Material extraction', choices=c('On', 'Off'), inline=TRUE),
+              ),
-                    radioButtons("pressure4", 'Abrasion of seabed', choices=c('On', 'Off'), inline=TRUE),
+              column(width=10, 
-                    radioButtons("pressure5", 'Penetration of seabed', choices=c('On', 'Off'), inline=TRUE),
+                h4('Effect on bio-assemblage'),
-                    radioButtons("pressure6", 'Siltation', choices=c('On', 'Off'), inline=TRUE),
+                plotlyOutput("layer1", height="270px") %>% withSpinner(),
-                    radioButtons("pressure7", 'Wave exposure', choices=c('On', 'Off'), inline=TRUE),
+                h4('Effect on Output Processes'),
-                    radioButtons("pressure8", 'Suspended sediment', choices=c('On', 'Off'), inline=TRUE),
+                plotlyOutput("layer2", height="270px") %>% withSpinner(),
-                    radioButtons("pressure9", 'Generic contamination', choices=c('On', 'Off'), inline=TRUE),
+                h4('Effect on Eco-system services'),
-                    radioButtons("pressure10", 'Deoxygenation', choices=c('On', 'Off'), inline=TRUE),
+                plotlyOutput("layer3", height="270px") %>% withSpinner()
-                    radioButtons("pressure11", 'Removal of target species', choices=c('On', 'Off'), inline=TRUE),
+              )
-                    actionButton("calcAB", "Calc")
+           )
-                  ),
-                  column(width=10, 
-                    h4('Effect on bio-assemblage'),
-                    plotlyOutput("layer1", height="270px") %>% withSpinner(),
-                    h4('Effect on Output Processes'),
-                    plotlyOutput("layer2", height="270px") %>% withSpinner(),
-                    h4('Effect on Eco-system services'),
-                    plotlyOutput("layer3", height="270px") %>% withSpinner()
-                  )
-               )
           ),
           tabItem(tabName = "2",h4("Bayesian Network"),
-              fluidPage(
+            fluidPage(
-                fluidRow(
+              fluidRow(
-                  plotOutput("bbnGraphPlot")
+                plotOutput("bbnGraphPlot")
+              ),
+              fluidRow(
+                column(
+                  width=6,
+                  h4('Ecoservice nodes'),
+                  DT::dataTableOutput('nodeTable')
                 ),
-                fluidRow(
+                column(
-                  column(
+                  width=6,
-                    width=6,
+                  h4('Ecoservice influences'),
-                    h4('Ecoservice nodes'),
+                  DT::dataTableOutput('edgeTable')
-                    DT::dataTableOutput('nodeTable')
-                  ),
-                  column(
-                    width=6,
-                    h4('Ecoservice influences'),
-                    DT::dataTableOutput('edgeTable')
-                  )
                 )
               )
+            )
           ),
           tabItem(tabName = "3",h4("Habitats"),
               fluidPage(
                 google_mapOutput(outputId = "map", width = "100%", height = "750px")
               )
-          )
+          ), 
+          tabItem(tabName = "4",h4("Ingestion"),
+             fluidPage(
+               p("Select a spreadsheet from your network for input into the JNCC Bayesian Network Analyser:"),
+               fileInput("fileSelect", "Choose Excel Spreadsheet File (xlsx format)", multiple = FALSE, accept = "xlsx"),
+               fluidRow(renderUI('status')),
+               actionButton('loadAB', 'Load') # icon='upload')
+             )
+           )
         )
       )
 )
@@ -100,6 +106,30 @@ server <- function(input, output, session) {
   dataStorage <- 'data/'
   
   models<-NULL
+  pressures <- NULL
+  
+  #disable(input$loadAb)
+  
+  .loadStatus <- reactiveValues(
+    valid = c(p=FALSE, ba=FALSE, op=FALSE, es=FALSE),
+    msgs = NULL
+  )
+  
+  .likelihoods <-reactiveValues(
+    p_ba = NULL,
+    ba_os = NULL,
+    os_es = NULL
+  )
+  
+  setPressures <- function(newPressures) {
+    pressures <<- newPressures
+  }
+  
+  validateSheets <- function() {
+    req(inputs$selectFile)
+    
+    ##TO DO - run parser on it and output the errors to 
+  }
   
   getAvailableModels <- function() {
     fileList <- list.files(dataStorage, pattern='.xlsx')
@@ -114,14 +144,11 @@ server <- function(input, output, session) {
       
       if (!is.null(tmp)) {
         modelList[[cnt]] <- tmp
-        
+
-        #tidy up the list for displaying
-        
         models <<- c(models, substr(fileList[idx], 1, (nchar(fileList[idx])-5)))
         print(paste('Model file successfully loaded', fileList[idx]))
         cnt=cnt+1
-        
+      
-        
       }
     }
     
@@ -133,91 +160,172 @@ server <- function(input, output, session) {
   
   .selections <- reactiveValues(model=1, layer=1)
   
-  #parse on load sheets in the input sheet folder
+  #parse on load sheets in the input sheet folder - replace with R Data
   modelList <- getAvailableModels()
   
-  calcLikelihood <- function(layer) {
+  calcLikelihood <- function(layer, pressStatus) {
     
-    isolate({
+    # isolate({
-      if (layer==1) layerStr='ba' else if (layer==2) layerStr='op' else layerStr ='es'
+    #   if (layer==1) layerStr='ba' else if (layer==2) layerStr='op' else layerStr ='es'
-      nodeList <- modelList[[.selections$model]][[.selections$layer]]$nodes
+    #   nodeList <- modelList[[.selections$model]][[.selections$layer]]$nodes
-      str(nodeList)
+    #   str(nodeList)
-      nodeNames <- nodeList$name[startsWith(nodeList$code, layerStr)]
+    #   nodeNames <- nodeList$name[startsWith(nodeList$code, layerStr)]
-      mean = runif(length(nodeNames), min=-1, max=1)
+    #   mean = runif(length(nodeNames), min=-1, max=1)
-      sd = runif(length(nodeNames), min=-0.25, max=0.25)
+    #   sd = runif(length(nodeNames), min=-0.25, max=0.25)
-      
-      df <- data.frame(
-        nodeNames = nodeNames,
-        range = c((mean - (3*sd)), (mean - (2*sd)), (mean - sd), mean,
-                  (mean + sd), (mean + (2*sd)), (mean + (3*sd))),
-        stringsAsFactors=FALSE
-      )
-      print(df)
-    })
-    return(
-      df
-    )
-    # isolate({   
-    #   
-    #   if (layer==1) layerStr='ba' else if (layer==2) layerStr='op' else if (layer==3) layerStr='es'
-    #     
-    #   layerRange <- which(startsWith(modelList[[.selections$model]][[layer]]$nodes, layerStr))
-    #   distList <- modelList[[.selections$model]][[layer]]$summDist[,layerRange]
-    #   nodeNames <- modelList[[.selections$model]][[layer]]$nodes$name[layerRange]
-    #     
-    #   }
-      
-    #   print(paste('Length of layer & node names',layer, length(nodeNames)))
-    #
-    #   distList <- modelList[[.selections$model]][[layer]]$summDist
-    #   colNames <- c('min', 'q1', 'q1', 'mean', 'q3', 'q3', 'max')
-    #   distM <- matrix(data=NA, nrow=ncol(distList), ncol=length(colNames))
-    #   
-    #   print(paste('Length of distributions',nrow(distM)))
-
-    #   for (col in 1:ncol(distList)) {
-    #     valsAsStrs <- unlist(strsplit(distList[,col], ":"))
-    #     valIdxs <- seq(from=2, to=length(valsAsStrs), by=2)
-    #     distVals <- as.numeric(valsAsStrs[valIdxs])
-    #     distM[col,] <- c(distVals[1], distVals[2], distVals[2], distVals[4], distVals[5], distVals[5], distVals[6])
-    #   }
-    # })
-    #
-    #  df <- data.frame(
-    #        nodeNames = nodeNames,
-    #        dist = distM,
-    #        stringsAsFactors=FALSE
-    #      )
-    #  print(df)
     # 
+    #   df <- data.frame(
+    #     nodeNames = nodeNames,
+    #     range = c((mean - (3*sd)), (mean - (2*sd)), (mean - sd), mean,
+    #               (mean + sd), (mean + (2*sd)), (mean + (3*sd))),
+    #     stringsAsFactors=FALSE
+    #   )
+    #   print(df)
+    # })
     # return(
-    #   df   
+    #   df
-    #    )
+    # )
+    isolate({
+
+      if (layer==1) layerStr='ba' else if (layer==2) layerStr='op' else if (layer==3) layerStr='es'
+
+      layerRange <- which(startsWith(modelList[[.selections$model]][[layer]]$nodes$code, layerStr))
+   
+      nodeCodes <- modelList[[.selections$model]][[layer]]$nodes$code[layerRange]
+      nodeNames <- modelList[[.selections$model]][[layer]]$nodes$name[layerRange]
+      
+      MEANPOS=1
+      MEANNEG=0
+      
+      # expr <- "("
+      # for (p in 1:nrow(pressStatus)) {
+      #   if (pressStatus$status[p] == 'On') {
+      #     threshold = MEANPOS
+      #   } else {
+      #     threshold = MEANNEG
+      #   }
+      #   
+      #   expr <- paste0(expr, "(\"", pressStatus$code[p], "\">=", threshold, ") & ")
+      # }
+      # expr <-substr(expr, 1, nchar(expr)-2)
+      # expr<-paste0(expr, ')')
+      # 
+      # print(expr)
+      
+      expr <- "list("
+      for (p in 1:nrow(pressStatus)) {
+        if (pressStatus$status[p] == 'On') {
+          threshold = MEANPOS
+        } else {
+          threshold = MEANNEG
+        }
+          
+        expr <- paste0(expr, "\"", pressStatus$code[p], "\"=", threshold, ", ")
+      }
+      expr <-substr(expr, 1, nchar(expr)-2)
+      expr<-paste0(expr, ')')
+      
+      print(expr)
+      
+      #txtStringWkg = "((\"p1\">=0.5) & (\"p10\">=0.5) & (\"p2\">=0.5))"  
+      
+      print(bnlearn::nodes(modelList[[.selections$model]][[layer]]$cfit))
+
+      sampleDists <- cpdist(
+        fitted = modelList[[.selections$model]][[layer]]$cfit,
+        nodes = bnlearn::nodes(modelList[[.selections$model]][[layer]]$cfit),
+        #evidence = eval(parse(text = expr)),
+        evidence = eval(parse(text = expr)),
+        method = "lw",
+        n = 10000,
+        debug=TRUE
+      )
+    })
+    
+    #print (sum(res[, 1] * attr(res, "weights")) / sum(attr(res, "weights")))
+    
+    print("Sample dists")
+    print(sampleDists)
+    
+    print("Weights")
+    print(unique(attr(sampleDists, "weights")))
+    displayCols <- match(nodeCodes, colnames(sampleDists))
+    sampleDists <- sampleDists[,displayCols]
+    means <- apply(sampleDists, 2, mean)
+    stdDev <- apply(sampleDists, 2, sd)
+    
+    print(modelList[[.selections$model]][[layer]]$nodes$name)
+      
+     return(data.frame(
+           nodeNames = nodeNames,
+           range = c(
+                  apply(sampleDists, 2, min),
+                  means - 2*stdDev,
+                  means - stdDev,
+                  means,
+                  means + stdDev,
+                  means + 2*stdDev,
+                  apply(sampleDists, 2, max)
+           ),
+           stringsAsFactors=FALSE
+         ))
+  }
+
+  renderStatus <- function(layer) {
+    isolate({
+      if (.loadStatus$valid[layer]) return('check-square') else return('times-circle')
+    })
   }
   
+  output$status <- renderUI({
+
+    tagList(
+      fluidRow(
+        column(width=3, h4('Pressures')),
+        column(width=3, h4('Bio-assemblages')),
+        column(width=3, h4('Output processes')),
+        column(width=3, h4('Eco-system services'))
+      ),
+     fluidRow(
+       column(width=3, icon(renderStatus(1))),
+       column(width=3, icon(renderStatus(2))),
+       column(width=3, icon(renderStatus(3))),
+       column(width=3, icon(renderStatus(4)))
+      )#,
+      #fluidRow(
+      #  verbatimTextOutput("msgBoard", .loadStatus$msg, placeholder=TRUE)
+      #)
+    )
+  })
   
-  .likelihoods <-reactiveValues(
+               
-    p_ba = calcLikelihood(1),
+  observeEvent(input$loadAB, {
-    ba_os = calcLikelihood(2),
+    #TO DO get spreadsheet
-    os_es = calcLikelihood(3)
+    #copy validated sheet into the data folder and either add or replace the sheet in the RData file
-  )
+    #reload the RData file
+    print('Load button pressed')
+  })             
  
   observeEvent(input$modelSelect, {
     .selections$model <<- match(input$modelSelect, models)
-    #print(.selections$model)
   })
   
   observeEvent(input$layerSelect, {
     .selections$layer <<- match(input$layerSelect, transitions)
-    #print(.selections$layer)
   })
   
+  
   observeEvent(input$calcAB, {
-    #print(paste('Action button pressed', input$calcAB))
+    #get the status of action buttons
+    isolate(myList <- reactiveValuesToList(input))
+    matches <- match(pressures$code,  names(myList))
+    status <-NULL
+    for (n in 1:length(matches))  status[n] = myList[[matches[n]]]
     
-      .likelihoods$p_ba <<- calcLikelihood(1)
+    pressStatus <- data.frame(code=pressures$code, status=status, stringsAsFactors = FALSE)
-      .likelihoods$ba_os <<- calcLikelihood(2)
+
-      .likelihoods$os_es <<- calcLikelihood(3)
+    .likelihoods$p_ba <<- calcLikelihood(1,  pressStatus)
+    .likelihoods$ba_os <<- calcLikelihood(2,  pressStatus)
+    .likelihoods$os_es <<- calcLikelihood(3, pressStatus)
 
   })
   
@@ -225,19 +333,32 @@ server <- function(input, output, session) {
     google_map(location = c(55, 0), zoom = 7)
   })
   
- 
+  makeRadioButtons <- function(row) {
+    radioButtons(row['code'], row['name'], choices=c('Off', 'On'), selected='Off', inline=TRUE)
+  }
+  
+  output$pressureList <- renderUI({
+    #isolate({
+    if (!is.null(modelList[[.selections$model]][[1]]$nodes)) {
+      pressCodes <- which(startsWith(modelList[[.selections$model]][[1]]$nodes$code, 'p'))
+      pressures <- data.frame(code = modelList[[.selections$model]][[1]]$nodes$code[pressCodes], 
+                     name = modelList[[.selections$model]][[1]]$nodes$name[pressCodes], stringsAsFactors=FALSE)
+      setPressures(pressures)
+      btnList <- apply(pressures, 1, makeRadioButtons)
+    }
+  })
   
   
   output$nodeTable <- DT::renderDataTable(
     
     modelList[[.selections$model]][[.selections$layer]]$nodes,
-      selection = 'single',options = list(searching = TRUE, pageLength = 10),server = TRUE, escape = FALSE,rownames= TRUE
+      selection = 'single',options = list(searching = TRUE, pageLength = 10, editable=TRUE),server = TRUE, escape = FALSE,rownames= TRUE
     )
 
   output$edgeTable <- DT::renderDataTable(
     
     modelList[[.selections$model]][[.selections$layer]]$edges,
-      selection = 'single',options = list(searching = TRUE, pageLength = 10),server = TRUE, escape = FALSE,rownames= TRUE
+      selection = 'single',options = list(searching = TRUE, pageLength = 10, editable=TRUE),server = TRUE, escape = FALSE,rownames= TRUE
   )
   
   output$bbnGraphPlot <- renderPlot({
@@ -245,21 +366,26 @@ server <- function(input, output, session) {
   })
   
   output$layer1 <- renderPlotly({
-    plot_ly(.likelihoods$p_ba, y = ~range, color = ~nodeNames, type = "box")
+    plot_ly(.likelihoods$p_ba, y = ~range, color = ~nodeNames, type = "box") %>%
+      layout(xaxis = list(zerolinewidth=2)) #%>%
+      #withSpinner()
+      
   })
 
   output$layer2 <- renderPlotly({
-    #print(.likelihoods)
-    
     if (.selections$layer>1) {
-      plot_ly(.likelihoods$ba_os, y = ~range, color = ~nodeNames, type = "box")
+      plot_ly(.likelihoods$ba_os, y = ~range, color = ~nodeNames, type = "box") %>%
+        layout(xaxis = list(zerolinewidth=2)) #%>%
+        #withSpinner()
     }
   })
 
   output$layer3 <- renderPlotly({
     
     if (.selections$layer>2) {
-      plot_ly(.likelihoods$os_es, y = ~range, color = ~nodeNames, type = "box")
+      plot_ly(.likelihoods$os_es, y = ~range, color = ~nodeNames, type = "box") %>%
+        layout(xaxis = list(zerolinewidth=2)) #%>%
+        #withSpinner()
     }
   })
 }