From aa480f906d3cbf897b7987b31e140f93a3041d91 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 11 Oct 2026 09:30:19 +0000 Subject: [PATCH] fix(fit): start gamma/lambda from the prior run's final values With prior = or a PM_result, estimated gamma/lambda restarted at the model's initial value, so continued runs and cycles = 0 posteriors used the wrong error model. Fixed error models are left unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01GTyUgJEd5uCEmv9qYrv8dC --- R/PM_model.R | 46 +++++++++++++++++-- tests/testthat/test-pm-fit-prior-gamma.R | 57 ++++++++++++++++++++++++ 2 files changed, 100 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-pm-fit-prior-gamma.R diff --git a/R/PM_model.R b/R/PM_model.R index 9b9003f6..17c150d2 100644 --- a/R/PM_model.R +++ b/R/PM_model.R @@ -1301,7 +1301,9 @@ PM_model <- R6::R6Class( } # set prior + prior_gamlam <- NULL # final gamma/lambda of a prior run, named by outeq if (inherits(prior, "PM_result")) { + prior_gamlam <- gamlam_from_cycle(prior$cycle$data$gamlam) prior <- prior$final$popPoints if (!is.data.frame(prior)) { msg <- c(msg, "{.arg prior} is a {.cls PM_result} without final population support points.") @@ -1316,6 +1318,7 @@ PM_model <- R6::R6Class( msg <- c(msg, "{.arg prior} file does not exist.", "i" = "Check the file path.") run_error <- run_error + 1 } + prior_gamlam <- gamlam_from_cycles_csv(glue::glue("{path}/{prior}/outputs/cycles.csv")) prior_file <- make_prior_file() file.copy(glue::glue("{path}/{prior}/outputs/theta.csv"), prior_file, overwrite = TRUE) prior_res <- normalize_prior_csv(prior_file, self$model_list$parameters) @@ -1454,6 +1457,21 @@ PM_model <- R6::R6Class( } + # start estimated gamma/lambda from the prior run's final values + error_models <- lapply(seq_along(self$model_list$err), function(i) { + e <- pm_err_flatten(self$model_list$err[[i]], i) + oq <- self$model_list$err[[i]]$outeq + if (length(oq) != 1 || is.na(oq)) oq <- i + start <- if (length(prior_gamlam) == 1) prior_gamlam[[1]] else prior_gamlam[as.character(oq)] + if (!isTRUE(e$fixed) && length(start) == 1 && !is.na(start)) { + e$initial <- unname(start) + } + e + }) + if (!is.null(prior_gamlam)) { + msg <- c(msg, "Initial gamma/lambda set to the final values of the prior run.") + } + if (intern) { ### CALL RUST out_path <- normalizePath(file.path(path_run, "outputs"), mustWork = FALSE) @@ -1466,9 +1484,7 @@ PM_model <- R6::R6Class( params = list( ranges = ranges, # not important but needed for POSTPROB algorithm = algorithm, - error_models = lapply(seq_along(self$model_list$err), function(i) { - pm_err_flatten(self$model_list$err[[i]], i) - }), + error_models = error_models, idelta = as.numeric(idelta), tad = as.numeric(tad), max_cycles = as.numeric(cycles), # will be hardcoded in Rust to 0 for POSTPROB @@ -2318,6 +2334,30 @@ pm_err_flatten <- function(e, position) { ) } +# Final gamma/lambda of a previous run, named by outeq. `gamlam` is the +# `$cycle$data$gamlam` tibble of a PM_result. +gamlam_from_cycle <- function(gamlam) { + if (!is.data.frame(gamlam) || nrow(gamlam) == 0) { + return(NULL) + } + last <- gamlam[gamlam$cycle == max(gamlam$cycle), ] + stats::setNames(last$value, last$outeq) +} + +# As gamlam_from_cycle(), read from a run's outputs/cycles.csv. +gamlam_from_cycles_csv <- function(file) { + if (!file.exists(file)) { + return(NULL) + } + cyc <- utils::read.csv(file, check.names = FALSE) + cols <- grep("^gamlam\\.", names(cyc), value = TRUE) + if (nrow(cyc) == 0 || length(cols) == 0) { + return(NULL) + } + outeq <- normalize_engine_index(as.numeric(sub("^gamlam\\.", "", cols))) + stats::setNames(unlist(cyc[nrow(cyc), cols]), outeq) +} + # Canonical output labels for a list of error models, in declaration order. pm_err_labels <- function(err) { vapply(seq_along(err), function(i) pm_err_label(err[[i]], i), character(1)) diff --git a/tests/testthat/test-pm-fit-prior-gamma.R b/tests/testthat/test-pm-fit-prior-gamma.R new file mode 100644 index 00000000..fcef800a --- /dev/null +++ b/tests/testthat/test-pm-fit-prior-gamma.R @@ -0,0 +1,57 @@ +library(Pmetrics) + +fit_capturing_error_models <- function(model, prior, path) { + captured <- NULL + testthat::with_mocked_bindings( + fit = function(model_source, data, params, output_path, solver = NULL) { + captured <<- params$error_models + dir.create(output_path, recursive = TRUE) + invisible(NULL) + }, + PM_parse = function(path) invisible(NULL), + PM_load = function(path, file, ...) structure(list(path = path), class = c("PM_result", "R6")), + model$fit( + data = NPex$data, path = path, run = 2, cycles = 0, prior = prior, + overwrite = TRUE, report = "none", quiet = TRUE + ) + ) + captured +} + +test_that("PM_model$fit starts gamma from the final value of a prior run number", { + run_root <- withr::local_tempdir(pattern = "prior-gamma-") + params <- tolower(names(NPex$model$model_list$pri)) + out_dir <- file.path(run_root, "1", "outputs") + dir.create(out_dir, recursive = TRUE) + theta <- data.frame(matrix(0.5, nrow = 2, ncol = length(params), dimnames = list(NULL, params)), prob = 0.5) + utils::write.csv(theta, file.path(out_dir, "theta.csv"), row.names = FALSE) + utils::write.csv( + data.frame(cycle = 1:2, neg2ll = c(10, 5), gamlam.0 = c(4, 2.5)), + file.path(out_dir, "cycles.csv"), + row.names = FALSE + ) + + model <- PM_model$new(x = NPex$model, compile = FALSE) + err <- fit_capturing_error_models(model, 1, run_root) + + expect_equal(err[[1]]$initial, 2.5) + # the model itself is unchanged + expect_false(isTRUE(all.equal(model$model_list$err[[1]]$initial, 2.5))) +}) + +test_that("PM_model$fit starts gamma from the final value of a PM_result prior", { + params <- tolower(names(NPex$model$model_list$pri)) + theta <- data.frame(matrix(0.5, nrow = 2, ncol = length(params), dimnames = list(NULL, params)), prob = 0.5) + prior <- structure( + list( + final = list(popPoints = theta), + cycle = list(data = list(gamlam = tibble::tibble(cycle = c(1, 2), value = c(4, 3.5), outeq = 1, type = "proportional"))) + ), + class = "PM_result" + ) + + model <- PM_model$new(x = NPex$model, compile = FALSE) + err <- fit_capturing_error_models(model, prior, withr::local_tempdir(pattern = "prior-gamma-res-")) + + expect_equal(err[[1]]$initial, 3.5) +})