Fix error in distribution of measures

Keep original behaviour of run_rabimo_with_measures()
by setting new arg "old_version" to TRUE

Author: hsonne

R/apply_measures_to_blocks.R

@@ -0,0 +1,204 @@
+# @param measures list with elements green_roof, unpaved, to_swale representing 
+#   the target percentages of the total areas corresponding to each measure.
+#   A value of NA means that the corresponding measure column is not touched.
+apply_measures_to_blocks <- function(blocks, measures, dbg = FALSE, check = FALSE)
+{
+  #dbg = FALSE; check = FALSE
+  
+  # Define helper functions
+  {
+    report_problem <- function(...) {
+      problem <- paste(unlist(list(...)), collapse = "")
+      warning(problem, call. = FALSE)
+    }
+    
+    debug <- function(...) {
+      if (dbg) {
+        writeLines(...)
+      }
+    }
+    
+    share_of_sum <- function(x) {
+      # Do not divide by zero, return vector of zeros instead.
+      if ((s <- sum(x)) == 0) {
+        rep(0, length(x))
+      } else {
+        x/s  
+      }
+    }
+    
+    check_if_target_was_reached <- function(blocks, measure) {
+      obtained <- kwb.rabimo::get_measure_stats(blocks)[[measure]]$mean
+      target <- measures[[measure]]
+      if (!isTRUE(all.equal(obtained, target))) {
+        warning(
+          sprintf("Target value %0.2f for '%s' ", target, measure),
+          sprintf("could not be achieved. Actual value: %0.2f", obtained), 
+          call. = FALSE
+        )
+      }
+    }
+    
+    check_for_negative_values <- function(blocks, measure) {
+      is_negative <- blocks[[measure]] < 0
+      if (any(is_negative)) {
+        warning(call. = FALSE, sprintf(
+          "There are %d negative values in column '%s'", 
+          sum(is_negative), measure
+        ))
+      }
+    }
+  }
+  
+  # The prefix "a_" refers to absolute area (in square metres)
+  
+  # Provide the total areas and roof areas in advance. They are not changed by 
+  # the measures.
+  a_total <- blocks$total_area
+  a_roof <- blocks$total_area * blocks$roof
+  
+  # 1. Handle measure "green roof"
+  if (!is.na(measures$green_roof)) {
+    
+    a_green_roof <- a_roof * blocks$green_roof
+    
+    # Total green roof area to add (if value >= 0) or to remove (if value < 0)
+    a_green_roof_change <- sum(measures$green_roof * a_total) - sum(a_green_roof)
+    
+    # Roof area that can be converted to green roof area
+    if (a_green_roof_change >= 0) {
+      
+      # increase green roof area
+      a_green_roof_potential <- a_roof - a_green_roof
+      
+      if (a_green_roof_change > sum(a_green_roof_potential)) {
+        report_problem(sprintf(
+          "Not enough (non-green) roof area available (%0.2f m2 missing)",
+          a_green_roof_change - sum(a_green_roof_potential)
+        ))
+      }
+      
+    } else {
+      
+      # decrease green roof area
+      a_green_roof_potential <- a_green_roof
+      
+      if (- a_green_roof_change > sum(a_green_roof_potential)) {
+        report_problem(sprintf(
+          "Not enough green roof area available (%0.2f m2 missing)",
+          - a_green_roof_change - sum(a_green_roof_potential)
+        ))
+      }    
+    }
+    
+    # Distribute change in green roof area to the different blocks    
+    a_green_roof_new <- a_green_roof + share_of_sum(a_green_roof_potential) * a_green_roof_change
+    
+    # Update column "green_roof" (as fraction of roof area)
+    blocks$green_roof <- ifelse(a_roof == 0, 0, a_green_roof_new / a_roof)
+  }
+  
+  # 2. Handle measure "Unsealing"
+  if (!is.na(measures$unpaved)) {
+    
+    # current paved/unpaved areas
+    a_paved <- a_total * blocks$pvd
+    a_unpaved <- a_total - a_roof - a_paved
+
+    # Required increase/decrease in unpaved area
+    a_unpaved_change <- measures$unpaved * sum(a_total) - sum(a_unpaved)
+    unpave <- a_unpaved_change >= 0
+    
+    debug(sprintf(
+      "%s area to be %s: %0.2f m2", 
+      ifelse(unpave, "Paved", "Unpaved"),
+      ifelse(unpave, "unpaved", "paved"),
+      abs(a_unpaved_change)
+    ))
+    
+    if (unpave) {
+      
+      a_potential <- a_paved
+      
+      if (a_unpaved_change > sum(a_potential)) {
+        report_problem(sprintf(
+          "Not enough paved area available to be unpaved (%0.2f m2 missing)",
+          a_unpaved_change - sum(a_potential)
+        ))
+      }
+      
+    } else {
+      
+      # actually pave instead of unpave!
+      a_potential <- a_unpaved
+      
+      # a_unpaved_change is negative here
+      if (- a_unpaved_change > sum(a_potential)) {
+        report_problem(sprintf(
+          "Not enough unpaved area available to be paved (%0.2f m2 missing)",
+          - a_unpaved_change - sum(a_potential)
+        ))
+      }
+    }
+    
+    # Distribute change in paved/unpaved area to the different blocks
+    a_paved_new <- a_paved - share_of_sum(a_potential) * a_unpaved_change
+    
+    blocks$pvd <- ifelse(a_total == 0, 0, a_paved_new / a_total)
+  }
+  
+  # 3. Handle measure "Connection to swales"
+  if (!is.na(measures$to_swale)) {
+    a_sealed <- a_roof + a_total * blocks$pvd
+    a_to_swale <- blocks$to_swale * a_sealed
+    a_to_swale_change <- sum(measures$to_swale * a_total) - sum(a_to_swale)
+    
+    increase <- a_to_swale_change >= 0
+    
+    debug(sprintf(
+      "Sealed area to be %s swales: %0.2f m2",
+      ifelse(increase, "connected to", "disconnected from"),
+      abs(a_to_swale_change)
+    ))
+    
+    if (increase) {
+      
+      a_to_swale_potential <- a_sealed - a_to_swale
+      
+      if (a_to_swale_change > sum(a_to_swale_potential)) {
+        report_problem(sprintf(
+          "Not enough sealed area available to be connected to swales (%0.2f m2 missing)",
+          a_to_swale_change - sum(a_to_swale_potential)
+        ))
+      }
+      
+    } else {
+      
+      # a_to_swale_change is negative here
+      a_to_swale_potential <- a_to_swale
+      
+      if (- a_to_swale_change > sum(a_to_swale_potential)) {
+        report_problem(sprintf(
+          "Not enough swale-connected sealed area available to be disconnected (%0.2f m2 missing)",
+          abs(a_to_swale_change) - sum(a_to_swale_potential)
+        ))
+      }
+    }
+    
+    # distribute
+    a_to_swale_new <- a_to_swale + share_of_sum(a_to_swale_potential) * a_to_swale_change
+    
+    # Update column "to_swale"
+    blocks$to_swale <- ifelse(a_sealed == 0, 0, a_to_swale_new / a_sealed)
+  }
+  
+  # Targets reached?
+  if (check) {
+    for (measure in names(measures)[!is.na(measures)]) {
+      check_if_target_was_reached(blocks, measure)  
+      check_for_negative_values(blocks, measure)
+    }
+  }
+  
+  blocks
+}

R/run_rabimo_with_measures.R

@@ -7,21 +7,24 @@
 #'   corresponding to each measure
 #' @param config configuration object, default:
 #'   \code{\link{rabimo_inputs_2020}$config}
+#' @param old_version if \code{TRUE} the old, erroneous version of this function 
+#'   is used (not correctly considering the updated pvd value before calculating 
+#'   the new to_swale values). The default is \code{FALSE}.
 #' @export
 run_rabimo_with_measures <- function(
     blocks,
     measures,
-    config = kwb.rabimo::rabimo_inputs_2020$config
+    config = kwb.rabimo::rabimo_inputs_2020$config,
+    old_version = FALSE
 )
 {
   #kwb.utils::assignPackageObjects("kwb.rabimo")
-  rescaled_targets <- rescale_target_values(
-    new_targets = measures,
-    blocks = blocks
-  )
-
-  run_rabimo(
-    distribute_measures(blocks = blocks, targets = rescaled_targets),
-    config = config
-  )
+  
+  new_blocks <- if (old_version) {
+    distribute_measures(blocks, rescale_target_values(measures, blocks))
+  } else {
+    apply_measures_to_blocks(blocks, measures)
+  }
+  
+  run_rabimo(new_blocks, config = config)
 }

man/run_rabimo_with_measures.Rd

@@ -0,0 +1,88 @@
+#library(testthat)
+test_that("apply_measures_to_blocks() works", {
+  
+  apply_measures_to_blocks <- kwb.rabimo:::apply_measures_to_blocks
+  
+  test_me <- function(data) {
+    
+    blocks <- data[sample(seq_len(nrow(data)), 10L), ]
+    stats <- kwb.rabimo:::get_measure_stats(blocks)
+    safety_factor <- 0.999
+    
+    measures_max <- list(
+      green_roof = safety_factor * stats$green_roof$max,
+      unpaved_max = safety_factor * stats$unpaved$max,
+      to_swale = NA
+    )
+    
+    measures_too_big_green_roof <- list(
+      green_roof = measures_max$green_roof + 0.01,
+      unpaved = measures_max$unpaved,
+      to_swale = NA
+    )
+    
+    measures_too_big_unpaved <- list(
+      green_roof = measures_max$green_roof,
+      unpaved = measures_max$unpaved + 0.01,
+      to_swale = NA
+    )
+    
+    expect_silent(b1 <- apply_measures_to_blocks(
+      blocks, 
+      measures = measures_max
+    ))
+    
+    expect_warning(b2 <- apply_measures_to_blocks(
+      blocks, 
+      measures = measures_too_big_green_roof
+    ))
+    
+    expect_warning(b3 <- apply_measures_to_blocks(
+      blocks, 
+      measures = measures_too_big_unpaved
+    ))
+  }
+  
+  test_me(data = kwb.rabimo::rabimo_inputs_2020$data)
+  test_me(data = kwb.rabimo::rabimo_inputs_2025$data)
+  
+})
+
+if (FALSE) {
+  
+  blocks <- rbind(
+    data.frame(total_area = 100, roof = 0.1, green_roof = 0, pvd = 0.7, to_swale = 1)
+    , data.frame(total_area = 200, roof = 0.2, green_roof = 0.8, pvd = 0.5, to_swale = 1)
+    , data.frame(total_area =  50, roof = 0.9, green_roof = 0.2, pvd = 0.1, to_swale = 1)
+  )
+  
+  stats <- kwb.rabimo::get_measure_stats(blocks)
+  str(stats)
+  
+  n <- 10L
+  
+  m <- as.data.frame(matrix(
+    runif(3*n), 
+    ncol = 3L, 
+    dimnames = list(NULL, c("green_roof", "unpaved", "to_swale"))
+  ))
+  
+  combinations <- split(m, seq_len(n))
+  
+  result <- lapply(combinations, function(measures) {
+    apply_measures_to_blocks(blocks, measures, dbg = TRUE, check = TRUE)
+  })
+  
+  measures <- list(
+    green_roof = 0, 
+    unpaved = 0,
+    to_swale = 0.1
+  )
+  
+  new_blocks <- apply_measures_to_blocks(blocks, measures)
+  
+  new_stats <- kwb.rabimo::get_measure_stats(new_blocks)
+  new_stats$green_roof$mean
+  new_stats$unpaved$mean
+  new_stats$to_swale$mean
+}