Skip to content

GitLab

Commit c77d8c22 authored by gross47's avatar gross47
Browse files

redo source

parent 2047d848
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 474 additions and 476 deletions
ui - old.R deleted 100644 → 0
View file @ 2047d848
#---#################################################################---#
#### Try to apply the stepwise linear approach on the Indonesia Data ####
#---#################################################################---#
# Thu Oct 15 15:33:34 2020 ------------------------------
# Volker und Kai
# See /rProgramming/indonesienOptimierung/Kopie von 08082018_Opti_Indon_Kai.xlsx
# Hinweis: Diese Optimierung hier unterscheided sich in der Definition des adjusted SEM vom einfachen Excel Bsp.
# Die optimistische Ausprägung ist gleich dem Land-Use Score (also mean) und nicht wie im einfachen Bsp. mean + sem
# (bzw. mean - sem bei "less is more").
# Deshalb habe ich die Funktion, die das adjSEM berechnet nochmal speziell für diese Optimierung angepasst. Besprechungsbedarf: Soll das
# generisch werden oder ist es ein Einzelfall? (Bei Indonesien ist es wie bei Gosling et al.)
#-------------------------------#
#### Load data and functions ####
#-------------------------------#
source("initScenario.R")
# source("calcDistanceToPerformanceScenario.R")
source("helper.R")
source("solveScenario.R")
source("functions.R")
library(lpSolveAPI)
library(tidyverse)
library(readxl)
library(shiny)
library(rsconnect)
library(shinyjs)
#-----------------#
#### Load data ####
#-----------------#
ui <- navbarPage(title = img(src="Logo_TUM_GOE.jpg", height = "40px", width = "250px"), id = "navBar",
theme = "test.css",
# inverse = TRUE,
windowTitle = "Robust multi-criterial optimization",
position = "fixed-top",
header = tags$style(
".navbar-right {
float: right !important;
}",
"body {padding-top: 75px;}"),
tabPanel("Home", value = "home",
shinyjs::useShinyjs(),
fluidRow(
HTML("
<section class='banner'>
<h2 class='parallax'>Robust multi-criterial optimization</h2>
<p class='parallax_description'></p>
</section>
")
),
fluidRow(
column(3),
column(6,
shiny::HTML("<br><br><center> <h1>What you'll find here</h1> </center><br>"),
shiny::HTML("<h5>Text</h5>")
),
column(3)
),
fluidRow(
style = "height:50px;"),
# PAGE BREAK
tags$hr(),
# HOW
fluidRow(
column(3),
column(6,
shiny::HTML("<br><br><center> <h1>How it can help you</h1> </center><br>"),
shiny::HTML("<h5>Text</h5>")
),
column(3)
),
fluidRow(
style = "height:50px;"),
# PAGE BREAK
tags$hr(),
# WHERE
fluidRow(
column(3),
column(6,
shiny::HTML("<br><br><center> <h1>Where it came from</h1> </center><br>"),
shiny::HTML("<h5>Text</h5>")
),
column(3)
),
fluidRow(
style = "height:50px;"),
# PAGE BREAK
tags$hr(),
# HOW TO START
fluidRow(
column(3),
column(6,
shiny::HTML("<br><br><center> <h1>How to get started</h1> </center><br>"),
shiny::HTML("<h5>Text.</h5>")
),
column(3)
)),
tabPanel("Upload file",
titlePanel("Uploading files"),
sidebarLayout(
sidebarPanel(
fileInput("file1", "Choose xlsx file", accept = ".xlsx"),
tags$br(),
checkboxInput('colnames', 'Colnames (experimental)', FALSE)
),
mainPanel(
tableOutput("contents")
)
)),
tabPanel(title = "Model",
headerPanel(""),
fluidRow(
column(4,
sidebarPanel(
checkboxGroupInput(inputId = "IndicatorGroups", "IndicatorGroups", "",
selected = "")
)),
column(6,
mainPanel(
plotOutput("plot1")
)),
column(2)),
tags$hr(),
fluidRow(
downloadButton('downloadPlot','Download Plot')
)
),
tabPanel(title = "About us")
)
ui - test.R 0 → 100644
View file @ c77d8c22
This diff is collapsed.
ui.R
View file @ c77d8c22
##--#####################################################--##
#### Transform the input table in an optimLanduse object ####
##--#####################################################--##
# Fri Jan 24 23:53:50 2020 ------------------------------
#' Initialize the robust optimization
#'
#' The function translates the indicators values and uncertainties for the land-use options into a solvable *optimLanduse* object.
#'
#' @param coefTable See the exemplary import folder.
#' @param uValue u Value.
#' @param optimisticRule Either *expectation* or *uncertaintyAdjustedExpectation*. It indicates weather the optimistic outcomes of an indicator are directly reflected by the *expectation* or if the indicator is *adjusted*.
#' @return An initialized landUse portfolio ready for optimization.
#' @import dplyr
#' @import tidyr
#' @export
initScenario <- function(coefTable, uValue = 3, optimisticRule = "expectation") {
#----------------------------#
#### Initialise the table ####
#----------------------------#
landUse <- as.character(unique(coefTable$landUse))
indicatorNames <- as.character(unique(coefTable$indicator))
expandList <- list()
expandList[landUse] <- list(c("High", "Low"))
expandMatrix1 <- as.matrix(expand.grid(expandList, stringsAsFactors = FALSE))
expandMatrix2 <- do.call(rbind, replicate(length(indicatorNames), expandMatrix1, simplify = FALSE))
scenarioTable <- tibble(indicator = rep(indicatorNames, each = dim(expandMatrix1)[1])) %>%
bind_cols(as_tibble(expandMatrix2))
# scenarioTable <- tibble(indicator = rep(indicatorNames, each = dim(expandMatrix1)[1])) %>%
# left_join(indicatorNames, by = "indicator") %>% bind_cols(as_tibble(expandMatrix2))
# Alter Version. Evtl relevant bei Fehlersuche. Ich weiß nicht mehr was ich mir bei dem left join gedacht habe.
# tbd. Tidy raus
scenarioTable <- scenarioTable %>% rename_at(.vars = vars(!!landUse[1] : !!landUse[length(landUse)]), .funs = funs(paste0("outcome", .)))
#--------------------#
## Attach direction ##
#--------------------#
scenarioTableTemp1 <- scenarioTable
scenarioTable <- merge(scenarioTable, unique(coefTable[, c("indicator","direction")]), by = "indicator")
if(!dim(scenarioTableTemp1)[1] == dim(scenarioTable)[1]) {cat("Error: Direction mising or wrong.")}
#---------------------------------------------#
## Attach indicator values and uncertainties ##
#---------------------------------------------#
scenarioTableTemp2 <- scenarioTable
# Das muss noch umgeschrieben werden, da funs "deprecated" ist: Am besten ganz ohne tidy ...
spread1 <- coefTable %>% select(-indicatorUncertainty) %>% spread(key = landUse, value = indicatorValue)
names(spread1)[names(spread1) %in% eval(landUse)] <- paste0("mean", names(spread1)[names(spread1) %in% eval(landUse)])
spread2 <- coefTable %>% select(-indicatorValue) %>% spread(key = landUse, value = indicatorUncertainty)
names(spread2)[names(spread2) %in% eval(landUse)] <- paste0("sem", names(spread2)[names(spread2) %in% eval(landUse)])
for(i in landUse) {
byIndicator <- c("indicator")
names(byIndicator) <- "indicator"
scenarioTable <- left_join(scenarioTable, spread1[, c("indicator", paste0("mean", i))], by = byIndicator)
scenarioTable <- left_join(scenarioTable, spread2[, c("indicator", paste0("sem", i))], by = byIndicator)
}
scenarioTable <- scenarioTable %>% select(-contains("mean"), everything()) # Order the variables, such that the means and uncertainties follow in direct succession
scenarioTable <- scenarioTable %>% select(-contains("sem"), everything()) # Alternatively, but slower, a second loop would be suitable
if(!dim(scenarioTableTemp1)[1] == dim(scenarioTable)[1]) {cat("Error: Attaching expectation or uncertainty failed.")}
#--------------------------------------------#
## Calculate indicator uncertainty adjusted ##
#--------------------------------------------#
scenarioTableTemp3 <- scenarioTable
newColumnNames <- paste0("adjSem", landUse)
scenarioTable[, newColumnNames] <- NA # Initialise empty
for(i in landUse) {
# Ugly. But fast and less error-prone
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("mean", i)] +
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("sem", i)] * uValue
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("mean", i)] -
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "Low", paste0("sem", i)] * uValue
if(optimisticRule == "uncertaintyAdjustedExpectation") {
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("mean", i)] -
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("sem", i)] * uValue
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("mean", i)] +
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("sem", i)] * uValue
}
if(optimisticRule == "expectation") {
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "less is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("mean", i)]
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("adjSem", i)] <-
scenarioTable[scenarioTable$direction == "more is better" & scenarioTable[, paste0("outcome", i)] == "High", paste0("mean", i)]
}
}
if(!optimisticRule %in% c("uncertaintyAdjustedExpectation", "expectation")) {cat("optimisticRule must be uncertaintyAdjustedExpectation or expectation")}
if(!dim(scenarioTableTemp3)[1] == dim(scenarioTable)[1] | any(is.na(scenarioTable))) {cat("Error: Calculation of adjusted uncertainty.")}
#--------------------------#
## calculate Min Max Diff ##
#--------------------------#
scenarioTable[, c("minAdjSem", "maxAdjSem", "diffAdjSem")] <-
apply(scenarioTable[, startsWith(names(scenarioTable), "adjSem")], 1, function(x) {c(min(x), max(x), (max(x) - min(x)))}) %>% t()
#-------------------------------------------------------------#
## Define the coefficients for the linear objective function ##
#-------------------------------------------------------------#
#and the restrictions.
coefObjective <- defineObjectiveCoefficients(scenarioTable)
#-------------------------------------#
#### Define the constraints matrix ####
#-------------------------------------#
constraintCoefficients <- defineConstraintCoefficients(scenarioTable)
retList <- list(scenarioSettings = data.frame(uValue = uValue,
optimisticRule = optimisticRule, stringsAsFactors = FALSE),
scenarioTable = scenarioTable,
coefObjective = coefObjective,
coefConstraint = constraintCoefficients,
status = "initialized",
beta = NA,
landUse = setNames(data.frame(matrix(rep(NA, length(landUse)), ncol = length(landUse), nrow = 1)), landUse),
optimDetails = list()
)
class(retList) <- "optimLanduse"
return(retList)
}
#----------------------------------------------------------#
#### Define the coefficients for the objective function ####
#----------------------------------------------------------#
# Calculate coefficients, one for each variable from the scenario table
#' @export
defineObjectiveCoefficients <- function(scenarioTable) {
# Set "less is better to negative" and divide by maximum difference
scenarioTable[scenarioTable$direction == "less is better", grep(c("^adjSem"), names(scenarioTable))] <-
-scenarioTable[scenarioTable$direction == "less is better", grep(c("^adjSem"), names(scenarioTable))]
scenarioTable[, grep(c("^adjSem"), names(scenarioTable))] <-
scenarioTable[, grep(c("^adjSem"), names(scenarioTable))] / scenarioTable$diffAdjSem * 100
coefObjective <- apply(scenarioTable[, grep(c("^adjSem"), names(scenarioTable))], 2, sum)
return(coefObjective)
}
#-------------------------------------#
#### Define the constraints matrix ####
#-------------------------------------#
#' @export
defineConstraintCoefficients <- function (scenarioTable) {
tempTableMore <- scenarioTable %>% filter(direction == "more is better") %>%
mutate_at(vars(starts_with("adjSem")), funs(modified = (. - minAdjSem) / diffAdjSem))
tempTableLess <- scenarioTable %>% filter(direction == "less is better") %>%
mutate_at(vars(starts_with("adjSem")), funs(modified = (maxAdjSem - .) / diffAdjSem))
tempTableMore %>% bind_rows(tempTableLess) %>% select(ends_with("modified")) %>% as.matrix() %>%
return()
}
##--#############################--##
#### Solve a optimLandUse object ####
##--#############################--##
# Mon Jan 27 10:35:54 2020 ------------------------------
#' Perform the robust optimization
#'
#' The function solves the initialized *optimLanduse* object.
#'
#' @param x The *optimLanduse* object.
#' @param digitsPrecision Precision.
#' @return A solved landUse portfolio.
#' @import lpSolveAPI
#' @export
solveScenario <- function (x, digitsPrecision = 4) {
# Bases on funDraft4 (rProgramming/uncertaintyOptimization/helperFunctions.R)
coefObjective <- x$coefObjective
piConstraintCoefficients <- x$coefConstraint
#tbd. Die Variablen sollte ich noch umbenennen. Von piConstraintCoefficients zu coefConstraint
precision <- 1 / 10^(digitsPrecision)
constraintCoef <- rbind(rep(1, length(coefObjective)), piConstraintCoefficients)
constraintDirection <- c("==", rep(">=", dim(piConstraintCoefficients)[1]))
piConstraintRhs <- c(0, .6, 1)
piConstraintRhsFeasible <- rep(FALSE, 3)
emergencyStop <- 1000
# Init lpa Object
lpaObj <- make.lp(nrow = 0, ncol = length(coefObjective))
set.objfn(lprec = lpaObj, obj = coefObjective)
add.constraint(lprec = lpaObj, xt = rep(1, length(coefObjective)), type = "=", rhs = 1)
apply(piConstraintCoefficients, 1, function(x) {add.constraint(lprec = lpaObj, xt =x, type = ">=", rhs = piConstraintRhs[2])})
lp.control(lprec = lpaObj, sense = "max")
counter <- 1 # 1 as the first iteration is outside the loop
set.rhs(lprec = lpaObj, b = c(1, rep(piConstraintRhs[2], dim(piConstraintCoefficients)[1])))
statusOpt <- lpSolveAPI::solve.lpExtPtr(lpaObj)
# Stepwise approximation loop
while (counter < emergencyStop) {
solutionFeasible <- TRUE
counter <- counter + 1
if (statusOpt == 0) {
piConstraintRhs <- c(piConstraintRhs[2], round((piConstraintRhs[2] + piConstraintRhs[3]) / 2, digitsPrecision), piConstraintRhs[3])
} else {
piConstraintRhs <- c(piConstraintRhs[1], round((piConstraintRhs[1] + piConstraintRhs[2]) / 2, digitsPrecision), piConstraintRhs[2])
}
set.rhs(lprec = lpaObj, b = c(1, rep(piConstraintRhs[2], dim(piConstraintCoefficients)[1])))
statusOpt <- lpSolveAPI::solve.lpExtPtr(lpaObj)
if(all(c(piConstraintRhs[3] - piConstraintRhs[2], piConstraintRhs[2] - piConstraintRhs[1]) <= precision)) {
break()
}
}
if(statusOpt == 2) {
set.rhs(lprec = lpaObj, b = c(1, rep(piConstraintRhs[1], dim(piConstraintCoefficients)[1])))
statusOpt <- lpSolveAPI::solve.lpExtPtr(lpaObj)
retPiConstraintRhs <- piConstraintRhs[1]
} else {
retPiConstraintRhs <- piConstraintRhs[2]
}
if(statusOpt != 0) {
cat("No optimum found.")
}
x$status <- "optimized"
x$beta <- round(retPiConstraintRhs, digitsPrecision)
x$landUse[1, ] <- get.variables(lpaObj)
return(x)
}
#### Excel-Manipulation function ####
library(readxl)
library(tidyverse)
library(janitor)
dataPreparation <- function(dat, uncertainty = "SE"){
## Filter all Rows with only NA ##
dat.final <- dat[rowSums(is.na(dat)) != ncol(dat), ]
dat.final <- dat.final %>% drop_na(...4)
## select landUse names ##
landUse <- dat[1, ]
landUse <- landUse[, colSums(is.na(landUse)) != nrow(landUse)]
landUse <- landUse %>% row_to_names(1)
landUse <- names(landUse)
## Create column names ##
dat.final <- dat.final %>% row_to_names(1)
## rename duplicated Columnnames ##
names(dat.final) <- make.unique(colnames(dat.final))
## rename first columns for initScenario function and define data structure ##
names(dat.final)[1:4] <- c("branch", "indicatorGroup", "indicator", "direction")
dat.final[, 5:ncol(dat.final)][is.na(dat.final[, 5:ncol(dat.final)])] <- 0
dat.final[5:ncol(dat.final)] <- lapply(dat.final[5:ncol(dat.final)], as.numeric)
## select mean values, rename columns and gather ##
importValues <- dat.final %>% select(branch, indicatorGroup, indicator, direction, starts_with("mean"))
colnames(importValues)[grepl("mean", colnames(importValues))] <- landUse
importValues <- importValues %>% gather(key = "landUse", value = "indicatorValue", Forest:`Oil palm plantation`)
## select uncertainty, rename columns and gather ##
importUnc <- dat.final %>%select(branch, indicatorGroup, indicator, direction, starts_with(uncertainty))
colnames(importUnc)[grepl(uncertainty, colnames(importUnc))] <- landUse
importUnc <- importUnc %>% gather(key = "landUse", value = "indicatorUncertainty", Forest:`Oil palm plantation`)
## combine mean and uncertainty ##
dataSource <- left_join(importValues, importUnc, by = c("branch", "indicatorGroup", "indicator", "direction", "landUse"))
return(dataSource)
}
#### Bray-Curtis dissimilarity ####
# BrayCurtis function - Knoke et al. 2015 (Compositional diversity of rehabilitated tropical lands supports multiple ecosystem services and buffers uncertainties)
# At the moment only usable for the exact same structure of data
BrayCurtis <- function(x){
# calculates the absolute difference between the different results of the total options and the desired (for example optimal) distribution
sum(abs(x$result - x$wantedResult)) / 200 * 100
}
#---#################################################################---#
#### Try to apply the stepwise linear approach on the Indonesia Data ####
......@@ -338,6 +20,10 @@ BrayCurtis <- function(x){
#-------------------------------#
#### Load data and functions ####
#-------------------------------#
source("initScenario.R")
source("solveScenario.R")
source("helper.R")
source("functions.R")
library(lpSolveAPI)
library(shiny)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment