fix(fit): start gamma/lambda from the prior run

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))

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)
+})