Optimisation #17
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -9,5 +9,6 @@ | |
| ^index\.md$ | ||
| ^README\.md$ | ||
| ^vignettes/index\.Rmd$ | ||
| ^vignettes/monotonicity_analysis$ | ||
| ^\.positai$ | ||
| ^\.claude$ | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,196 @@ | ||
| #' Smallest feasible parameter value via bisection (monotone threshold search) | ||
| #' | ||
| #' Core building block of the swale-design optimiser: finds the smallest | ||
| #' value of one design parameter for which the overflow target is met | ||
| #' (`n_overflows <= x_max`), assuming quasi-monotone feasibility (larger | ||
| #' value = never more overflows; verified for the RAINDROP model in the | ||
| #' `monotonicity_analysis` vignette). Each evaluation halves the search | ||
| #' interval, so `ceiling(log2(range / tol))` evaluations suffice. | ||
| #' | ||
| #' Two safety rules from the monotonicity analysis are built in: | ||
| #' \itemize{ | ||
| #' \item \strong{Edge guard}: if the upper bound is infeasible by no more | ||
| #' than `wobble` events (the +1 counting artefact of the 4-h event | ||
| #' separation), a descending ladder below the edge searches for a | ||
| #' feasible anchor before the branch is declared infeasible. | ||
| #' \item \strong{Volume referee}: whenever `n_overflows` increases with | ||
| #' the parameter (a counting flip), the overflow volume must have | ||
| #' decreased; if the volume increased as well, a warning is emitted and | ||
| #' `monotonicity_violation` is set (real non-monotonicity -- never | ||
| #' observed at the three validation sites). | ||
| #' } | ||
| #' | ||
| #' @param evaluate `function(value)` returning a list / one-row data.frame | ||
| #' with at least `n_overflows`; if it also contains `volume_column`, the | ||
| #' volume referee is active. Evaluations are memoised per value. | ||
| #' @param x_max Feasibility target: feasible iff `n_overflows <= x_max`. | ||
| #' @param lower,upper Numeric search bounds (continuous mode). | ||
| #' @param tol Resolution of the continuous search (same unit as the value). | ||
| #' @param levels Sorted numeric vector of discrete candidate values | ||
| #' (discrete mode, e.g. Sickerbox stack heights). If given, `lower`, | ||
| #' `upper` and `tol` are ignored and a binary search over the levels is | ||
| #' performed. | ||
| #' @param wobble Maximum counting-artefact size tolerated by the edge guard | ||
| #' (default 1, matching the observed +1 flips). | ||
| #' @param volume_column Name of the volume element in the `evaluate` result | ||
| #' used by the volume referee (default `"overflow_volume_m3"`). | ||
| #' @param verbose Print one line per evaluation. | ||
| #' | ||
| #' @return List with | ||
| #' \describe{ | ||
| #' \item{value}{smallest feasible value, or `NA` if infeasible} | ||
| #' \item{n_overflows}{overflow count at `value`} | ||
| #' \item{status}{`"ok"`, `"at_lower_bound"` (already feasible at the lower | ||
| #' end -- caller may widen the bracket) or `"infeasible"`} | ||
| #' \item{evaluations}{tibble of all evaluated values (value, n_overflows, | ||
| #' volume), sorted by value} | ||
| #' \item{n_evaluations}{number of distinct evaluations} | ||
| #' \item{monotonicity_violation}{`TRUE` if the volume referee fired} | ||
| #' } | ||
| #' | ||
| #' @examples | ||
| #' # synthetic monotone step function: feasible from 137.4 m2 on | ||
| #' f <- function(v) list(n_overflows = if (v >= 137.4) 0L else 10L) | ||
| #' find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2)$value | ||
| #' | ||
| #' @export | ||
| find_min_feasible <- function(evaluate, | ||
| x_max, | ||
| lower = NULL, | ||
| upper = NULL, | ||
| tol = 1, | ||
| levels = NULL, | ||
| wobble = 1L, | ||
| volume_column = "overflow_volume_m3", | ||
| verbose = FALSE) { | ||
|
|
||
| discrete <- !is.null(levels) | ||
| if (discrete) { | ||
| grid <- sort(unique(levels)) | ||
| axis_lo <- 0 # virtual: below the smallest level | ||
| axis_hi <- length(grid) | ||
| axis_tol <- 1 | ||
| to_value <- function(i) grid[[i]] | ||
| } else { | ||
| stopifnot(is.numeric(lower), is.numeric(upper), upper > lower, tol > 0) | ||
| axis_lo <- lower | ||
| axis_hi <- upper | ||
| axis_tol <- tol | ||
| to_value <- identity | ||
| } | ||
|
|
||
| evals <- new.env(parent = emptyenv()) | ||
| eval_at <- function(axis_pos) { | ||
| v <- to_value(axis_pos) | ||
| key <- format(v, digits = 15) | ||
| if (!is.null(evals[[key]])) return(evals[[key]]) | ||
| res <- as.list(evaluate(v)) | ||
| if (!"n_overflows" %in% names(res)) { | ||
| stop("find_min_feasible(): evaluate() must return an element 'n_overflows'") | ||
| } | ||
| row <- list( | ||
| value = v, | ||
| n_overflows = as.numeric(res$n_overflows), | ||
| volume = if (volume_column %in% names(res)) { | ||
| as.numeric(res[[volume_column]]) | ||
| } else { | ||
| NA_real_ | ||
| } | ||
| ) | ||
| assign(key, row, envir = evals) | ||
| if (isTRUE(verbose)) { | ||
| message(sprintf(" eval %s -> n_overflows = %s", | ||
| format(v), format(row$n_overflows))) | ||
| } | ||
| row | ||
| } | ||
| feasible <- function(axis_pos) { | ||
| n <- eval_at(axis_pos)$n_overflows | ||
| !is.na(n) && n <= x_max | ||
| } | ||
|
|
||
| status <- "ok" | ||
|
|
||
| # 1) Upper edge: with monotone feasibility the whole range is infeasible | ||
| # if the upper edge is -- unless the edge is a +1 counting wobble, in | ||
| # which case a descending ladder looks for a feasible anchor below. | ||
| if (!feasible(axis_hi)) { | ||
| n_hi <- eval_at(axis_hi)$n_overflows | ||
| anchor <- NA_real_ | ||
| if (!is.na(n_hi) && n_hi <= x_max + wobble) { | ||
| offset <- axis_tol | ||
| repeat { | ||
| p <- axis_hi - offset | ||
| if (discrete) p <- ceiling(p) | ||
| if (p <= axis_lo) break | ||
| if (feasible(p)) { | ||
| anchor <- p | ||
| break | ||
| } | ||
| offset <- offset * 2 | ||
| } | ||
| } | ||
| if (is.na(anchor)) { | ||
| status <- "infeasible" | ||
| } else { | ||
| axis_hi <- anchor | ||
| } | ||
| } | ||
|
|
||
| # 2) Bisection: invariant lo infeasible (or virtual/edge), hi feasible. | ||
| best <- NA_real_ | ||
| if (!identical(status, "infeasible")) { | ||
| lo <- axis_lo | ||
| hi <- axis_hi | ||
| if (!discrete && feasible(lo)) { | ||
| hi <- lo # optimum at (or below) the lower end | ||
| } | ||
| while ((hi - lo) > axis_tol) { | ||
| mid <- (lo + hi) / 2 | ||
| if (discrete) mid <- floor(mid) | ||
| if (mid <= lo || mid >= hi) break | ||
| if (feasible(mid)) hi <- mid else lo <- mid | ||
| } | ||
| best <- hi | ||
| at_lower <- if (discrete) best <= 1 else best <= axis_lo | ||
| if (at_lower) status <- "at_lower_bound" | ||
| } | ||
|
|
||
| # 3) Volume referee over all evaluations of this search | ||
| ev <- do.call(rbind, lapply(ls(evals), function(k) { | ||
| e <- get(k, envir = evals) | ||
| data.frame(value = e$value, n_overflows = e$n_overflows, | ||
| volume = e$volume) | ||
| })) | ||
| ev <- ev[order(ev$value), , drop = FALSE] | ||
| rownames(ev) <- NULL | ||
| monotonicity_violation <- FALSE | ||
| if (nrow(ev) >= 2) { | ||
| dn <- diff(ev$n_overflows) | ||
| dv <- diff(ev$volume) | ||
| flips <- which(!is.na(dn) & dn > 0) | ||
| real <- flips[!is.na(dv[flips]) & dv[flips] > 1e-9] | ||
| if (length(real) > 0) { | ||
| monotonicity_violation <- TRUE | ||
| warning(sprintf( | ||
| paste0("find_min_feasible(): n_overflows AND overflow volume ", | ||
| "increase between value %s and %s -- real non-monotonicity, ", | ||
| "bisection result unreliable for this branch."), | ||
| format(ev$value[real[1]]), format(ev$value[real[1] + 1]) | ||
| ), call. = FALSE) | ||
| } | ||
| } | ||
|
Comment on lines
+168
to
+182
There was a problem hiding this comment. Volume-referee blind spot when either endpoint's volume is Consider distinguishing "clean" from "unknown" — e.g., warn "volume unavailable at value X → referee disabled for this flip" and set |
||
|
|
||
| list( | ||
| value = if (identical(status, "infeasible")) NA_real_ else to_value(best), | ||
| n_overflows = if (identical(status, "infeasible")) { | ||
| NA_real_ | ||
| } else { | ||
| eval_at(best)$n_overflows | ||
| }, | ||
| status = status, | ||
| evaluations = tibble::as_tibble(ev), | ||
| n_evaluations = nrow(ev), | ||
| monotonicity_violation = monotonicity_violation | ||
| ) | ||
| } | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Edge-guard ladder never evaluates
axis_loitself. In continuous mode,if (p <= axis_lo) breakexits without probingaxis_lo. Under strict monotonicity this is safe (iffeasible(axis_lo)were true, the very first stepaxis_hi - tolwould already succeed), and the follow-upif (!discrete && feasible(lo))guard at line 145 catches the corner case after an anchor is found. But if the upper edge is truly infeasible-by-wobble AND no descending step lands on a feasible value (onlyaxis_loitself is feasible), the search returns"infeasible"when a lower-bound solution existed.This is an unlikely combination given the +1 wobble assumption, and I don't have a concrete failure scenario for the RAINDROP model — but the invariant is worth stating in the roxygen block, so future callers understand the assumption the guard rests on.