From b848fb360aa865c37298ee7526344b5280c700cc Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 17 Oct 2022 10:28:54 +0200
Subject: SFORB in saem, update for mhmkin and multistart

---
 R/mhmkin.R     |  4 ++++
 R/multistart.R | 29 +++++++++++++++++++++++++++++
 R/saem.R       | 37 ++++++++++++++++++++++++++++++++++++-
 3 files changed, 69 insertions(+), 1 deletion(-)

(limited to 'R')

diff --git a/R/mhmkin.R b/R/mhmkin.R
index 2df01f0c..542bddf7 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -207,6 +207,10 @@ BIC.mhmkin <- function(object, ...) {
 #' @export
 update.mhmkin <- function(object, ..., evaluate = TRUE) {
   call <- attr(object, "call")
+  # For some reason we get mhkin.list in call[[1]] when using mhmkin from the
+  # loaded package so we need to fix this so we do not have to export
+  # mhmkin.list in addition to the S3 method mhmkin
+  call[[1]] <- mhmkin
 
   update_arguments <- match.call(expand.dots = FALSE)$...
 
diff --git a/R/multistart.R b/R/multistart.R
index a788953e..11736670 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -71,6 +71,7 @@ multistart <- function(object, n = 50,
 #' @export
 multistart.saem.mmkin <- function(object, n = 50, cores = 1,
   cluster = NULL, ...) {
+  call <- match.call()
   if (n <= 1) stop("Please specify an n of at least 2")
 
   mmkin_parms <- parms(object$mmkin, errparms = FALSE,
@@ -90,6 +91,7 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1,
   }
   attr(res, "orig") <- object
   attr(res, "start_parms") <- start_parms
+  attr(res, "call") <- call
   class(res) <- c("multistart.saem.mmkin", "multistart")
   return(res)
 }
@@ -178,3 +180,30 @@ which.best.default <- function(object, ...)
   ll <- sapply(object, llfunc)
   return(which.max(ll))
 }
+
+#' @export
+update.multistart <- function(object, ..., evaluate = TRUE) {
+  call <- attr(object, "call")
+  # For some reason we get multistart.saem.mmkin in call[[1]] when using multistart
+  # from the loaded package so we need to fix this so we do not have to export
+  # multistart.saem.mmkin
+  call[[1]] <- multistart
+
+  update_arguments <- match.call(expand.dots = FALSE)$...
+
+  if (length(update_arguments) > 0) {
+    update_arguments_in_call <- !is.na(match(names(update_arguments), names(call)))
+  }
+
+  for (a in names(update_arguments)[update_arguments_in_call]) {
+    call[[a]] <- update_arguments[[a]]
+  }
+
+  update_arguments_not_in_call <- !update_arguments_in_call
+  if(any(update_arguments_not_in_call)) {
+    call <- c(as.list(call), update_arguments[update_arguments_not_in_call])
+    call <- as.call(call)
+  }
+  if(evaluate) eval(call, parent.frame())
+  else call
+}
diff --git a/R/saem.R b/R/saem.R
index 5256f6b5..dbc19dec 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -313,7 +313,7 @@ saemix_model <- function(object, solution_type = "auto",
     # Parent only
     if (length(mkin_model$spec) == 1) {
       parent_type <- mkin_model$spec[[1]]$type
-      if (length(odeini_fixed) == 1) {
+      if (length(odeini_fixed) == 1 && !grepl("_bound$", names(odeini_fixed))) {
         if (transformations == "saemix") {
           stop("saemix transformations are not supported for parent fits with fixed initial parent value")
         }
@@ -344,6 +344,9 @@ saemix_model <- function(object, solution_type = "auto",
           }
         }
       } else {
+        if (length(odeini_fixed) == 2) {
+          stop("SFORB with fixed initial parent value is not supported")
+        }
         if (parent_type == "SFO") {
           if (transformations == "mkin") {
             model_function <- function(psi, id, xidep) {
@@ -386,6 +389,38 @@ saemix_model <- function(object, solution_type = "auto",
             transform.par = c(0, 1, 1, 3)
           }
         }
+        if (parent_type == "SFORB") {
+          if (transformations == "mkin") {
+            model_function <- function(psi, id, xidep) {
+              k_12 <- exp(psi[id, 3])
+              k_21 <- exp(psi[id, 4])
+              k_1output <- exp(psi[id, 2])
+              t <- xidep[, "time"]
+
+              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
+              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
+              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
+
+              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
+                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
+            }
+          } else {
+            model_function <- function(psi, id, xidep) {
+              k_12 <- psi[id, 3]
+              k_21 <- psi[id, 4]
+              k_1output <- psi[id, 2]
+              t <- xidep[, "time"]
+
+              sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + k_1output) * k_21)
+              b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
+              b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
+
+              psi[id, 1] * (((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
+                ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t))
+            }
+            transform.par = c(0, 1, 1, 1)
+          }
+        }
         if (parent_type == "HS") {
           if (transformations == "mkin") {
             model_function <- function(psi, id, xidep) {
-- 
cgit v1.2.1