Merge pull request #232 from canmod/tv-stitch

* more trajectory methods
* parameterized trajectory
* optimized spec
* forecaster version of a simulator
* check what functions are used
* expose print out of parameterization info
* sim offset and much more robust sim offset/sim bounds
* insert log linear (experimental)
* insert basic transformations of trajectories
* expose objective function
* more consistently use best last parameter
* wrap more optimizers
* should we really reset best?
* current$data_frame() should give best.last.par
* make my own term column -- keep an eye on this
* only back transform if there are more than zero coefficientss
* check if opt_attempted
* potentiallly widespread effects of bind_rows fix (keep an eye on this)
* two new engine functions -- getting last item in a vector and checking if a matrix is finite

Author: stevencarlislewalker

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: macpan2
 Title: Fast and Flexible Compartmental Modelling
-Version: 1.15.3
+Version: 1.16.0
 Authors@R: c(
   person("Steve Walker", email="swalk@mcmaster.ca", role=c("cre", "aut")),
   person("Weiguang Guan", role="aut"),
@@ -37,6 +37,7 @@ Suggests:
     visNetwork,
     kableExtra,
     broom.mixed,
+    DEoptim,
     outbreaks,
     tmbstan,
     numDeriv,

Makefile

@@ -53,6 +53,7 @@ quick-test-all:
 	make run-vignette-code
 	make run-tests
 	make run-examples
+	make run-starter-readme-code
 
 quick-test:
 	make quick-doc-install

NAMESPACE

@@ -38,8 +38,10 @@ S3method(length,StructuredVector)
 S3method(macpan2::TMBPar,ParArg)
 S3method(macpan2::TMBPar,character)
 S3method(macpan2::TMBPar,list)
+S3method(mp_default,TMBCalibrator)
 S3method(mp_default,TMBModelSpec)
 S3method(mp_default,TMBSimulator)
+S3method(mp_default_list,TMBCalibrator)
 S3method(mp_default_list,TMBModelSpec)
 S3method(mp_default_list,TMBSimulator)
 S3method(mp_dynamic_simulator,DynamicModel)
@@ -53,16 +55,23 @@ S3method(mp_expand,TMBModelSpec)
 S3method(mp_extract,DynamicModel)
 S3method(mp_extract,Ledger)
 S3method(mp_extract,ModelDefRun)
+S3method(mp_final,TMBCalibrator)
 S3method(mp_final,TMBSimulator)
 S3method(mp_final_list,TMBSimulator)
 S3method(mp_fit,character)
 S3method(mp_fit,numeric)
+S3method(mp_functions_used,ExprList)
+S3method(mp_functions_used,TMBCalibrator)
+S3method(mp_functions_used,TMBModelSpec)
+S3method(mp_functions_used,TMBSimulator)
 S3method(mp_hazard,TMBModelSpec)
 S3method(mp_index,character)
 S3method(mp_index,data.frame)
 S3method(mp_index,numeric)
+S3method(mp_initial,TMBCalibrator)
 S3method(mp_initial,TMBModelSpec)
 S3method(mp_initial,TMBSimulator)
+S3method(mp_initial_list,TMBCalibrator)
 S3method(mp_initial_list,TMBModelSpec)
 S3method(mp_initial_list,TMBSimulator)
 S3method(mp_labels,Index)
@@ -71,7 +80,11 @@ S3method(mp_nofit,character)
 S3method(mp_nofit,numeric)
 S3method(mp_optimize,TMBCalibrator)
 S3method(mp_optimize,TMBSimulator)
+S3method(mp_optimized_spec,TMBCalibrator)
+S3method(mp_optimized_spec,TMBSimulator)
 S3method(mp_optimizer_output,TMBCalibrator)
+S3method(mp_parameterization,TMBCalibrator)
+S3method(mp_parameterization,TMBSimulator)
 S3method(mp_reduce,TMBModelSpec)
 S3method(mp_reference,Index)
 S3method(mp_reference,Ledger)
@@ -94,13 +107,18 @@ S3method(mp_tmb_coef,TMBCalibrator)
 S3method(mp_tmb_coef,TMBSimulator)
 S3method(mp_tmb_fixef_cov,TMBCalibrator)
 S3method(mp_tmb_fixef_cov,TMBSimulator)
+S3method(mp_tmb_objective,TMBCalibrator)
+S3method(mp_tmb_objective,TMBSimulator)
 S3method(mp_tmb_test,TMBModelSpec)
 S3method(mp_tmbstan_coef,TMBCalibrator)
 S3method(mp_tmbstan_coef,TMBSimulator)
 S3method(mp_trajectory,TMBCalibrator)
 S3method(mp_trajectory,TMBSimulator)
 S3method(mp_trajectory_ensemble,TMBCalibrator)
 S3method(mp_trajectory_ensemble,TMBSimulator)
+S3method(mp_trajectory_par,TMBCalibrator)
+S3method(mp_trajectory_par,TMBSimulator)
+S3method(mp_trajectory_replicate,TMBCalibrator)
 S3method(mp_trajectory_replicate,TMBSimulator)
 S3method(mp_trajectory_sd,TMBCalibrator)
 S3method(mp_trajectory_sd,TMBSimulator)
@@ -180,6 +198,8 @@ S3method(to_string,default)
 S3method(to_string,formula)
 S3method(to_values,data.frame)
 S3method(to_values,numeric)
+S3method(updated_param_vector,data.frame)
+S3method(updated_param_vector,list)
 export(BinaryOperator)
 export(CSVReader)
 export(JSONReader)
@@ -223,6 +243,9 @@ export(mp_final_list)
 export(mp_fit)
 export(mp_flow_frame)
 export(mp_flow_vars)
+export(mp_forecaster)
+export(mp_functions_used)
+export(mp_generates_randomness)
 export(mp_group)
 export(mp_hazard)
 export(mp_identity)
@@ -247,8 +270,10 @@ export(mp_neg_bin)
 export(mp_nofit)
 export(mp_normal)
 export(mp_optimize)
+export(mp_optimized_spec)
 export(mp_optimizer_output)
 export(mp_par)
+export(mp_parameterization)
 export(mp_per_capita_flow)
 export(mp_per_capita_inflow)
 export(mp_per_capita_outflow)
@@ -268,11 +293,13 @@ export(mp_set_numbers)
 export(mp_setdiff)
 export(mp_show_models)
 export(mp_sim_bounds)
+export(mp_sim_offset)
 export(mp_simulator)
 export(mp_slices)
 export(mp_sqrt)
 export(mp_square)
 export(mp_state_dependence_frame)
+export(mp_state_flow_vars)
 export(mp_state_vars)
 export(mp_structured_vector)
 export(mp_subset)
@@ -286,14 +313,18 @@ export(mp_tmb_entire_library)
 export(mp_tmb_expr_list)
 export(mp_tmb_fixef_cov)
 export(mp_tmb_insert)
+export(mp_tmb_insert_log_linear)
 export(mp_tmb_insert_reports)
+export(mp_tmb_insert_trans)
 export(mp_tmb_library)
 export(mp_tmb_model_spec)
+export(mp_tmb_objective)
 export(mp_tmb_update)
 export(mp_tmbstan_coef)
 export(mp_traj)
 export(mp_trajectory)
 export(mp_trajectory_ensemble)
+export(mp_trajectory_par)
 export(mp_trajectory_replicate)
 export(mp_trajectory_sd)
 export(mp_trajectory_sim)
@@ -340,6 +371,7 @@ importFrom(oor,return_object)
 importFrom(stats,as.formula)
 importFrom(stats,make.link)
 importFrom(stats,na.omit)
+importFrom(stats,plogis)
 importFrom(stats,qlogis)
 importFrom(stats,qnorm)
 importFrom(stats,quantile)

R/cal_utils.R

@@ -0,0 +1,6 @@
+get_last_best_par = function(ad_fun) {
+  best_par = ad_fun$env$last.par.best
+  ranef_indices = ad_fun$env$random
+  fixef_indices = setdiff(seq_along(best_par), ranef_indices)
+  best_par[fixef_indices] |> rep_name("params")
+}

R/calibration.R

@@ -17,14 +17,6 @@ TMBOptimizer = function(simulator) {
     force(arg_mappings)
     force(opt_func)
     force(opt_method_nm)
-    ## TODO: add package dependencies to assert using assert_dependencies
-    
-    get_last_best_par = function(ad_fun) {
-      best_par = ad_fun$env$last.par.best
-      ranef_indices = ad_fun$env$random
-      fixef_indices = setdiff(seq_along(best_par), ranef_indices)
-      best_par[fixef_indices]
-    }
 
     self[[opt_method_nm]] = function() {
 
@@ -63,8 +55,7 @@ TMBOptimizer = function(simulator) {
       ## invalidate the now out-of-date sdreport
       self$simulator$cache$sdreport$invalidate() 
 
-      #ad_fun$fn(opt_obj$par) ## probably this should be last.par.best
-      ad_fun$fn(get_last_best_par(ad_fun))
+      self$simulator$objective(get_last_best_par(ad_fun))
 
       opt_obj
     }
@@ -78,6 +69,32 @@ TMBOptimizer = function(simulator) {
       "nlminb", stats::nlminb
     , par = "start", fn = "objective", gr = "gradient", he = "hessian"
   )
+  wrap(
+      "DEoptim", DEoptim::DEoptim
+    , fn = "fn"
+  )
+  wrap(
+      "optimize", stats::optimize
+    , fn = "f"
+  )
+  wrap(
+      "optimise", stats::optimize
+    , fn = "f"
+  )
+  
+  self$extract_best = list(
+      nlminb = \(obj) obj$par
+    , optim = \(obj) obj$par
+    , optimize = \(obj) obj$minimum
+    , optimise = \(obj) obj$minimum
+    , DEoptim = \(obj) obj$optim$bestmem
+  )
+  self$reset_best = function(opt_obj, optimizer) {
+    ## this is useful to make sure that the names
+    ## of the parameter vector do not get mangled
+    ## by the optimizers
+    self$simulator$objective(self$extract_best[[optimizer]](opt_obj))
+  }
 
   return_object(self, "TMBOptimizer")
 }
@@ -94,7 +111,8 @@ TMBCurrentParams = function(simulator) { ## TMBSimulator
   }
   self$params_vector = function() {
     if (self$n_params() == 0L) return(numeric())
-    self$simulator$ad_fun()$env$parList()$params
+    get_last_best_par(self$simulator$ad_fun())
+    # $env$parList()$params
   }
   self$random_vector = function() {
     if (self$n_random() == 0L) return(numeric())
@@ -107,6 +125,58 @@ TMBCurrentParams = function(simulator) { ## TMBSimulator
   self$random_frame = function() {
     self$simulator$tmb_model$random$data_frame(current = self$random_vector())
   }
+  
+  ## update a matrix_list, which is a list of numeric matrices (or numeric 
+  ## vectors, which will be treated as n-by-1 matrices) so that the current 
+  ## fitted values of parameters are used to replace associated values in the 
+  ## matrix_list
+  self$update_matrix_list = function(matrix_list) {
+    current_frame = rbind(self$params_frame(), self$random_frame())
+    keepers = current_frame$mat %in% names(matrix_list)
+    current_frame = current_frame[keepers, , drop = FALSE]
+    if (nrow(current_frame) == 0L) return(matrix_list)
+    for (i in seq_len(nrow(current_frame))) {
+      mat = current_frame[i, "mat"]
+      mat = mat[]
+      
+      row = current_frame[i, "row"]
+      nr = nrow(matrix_list[[mat]])
+      if (is.null(nr)) nr = length(matrix_list[[mat]])
+      if (row > nr) {
+        stop(
+          sprintf(
+              "Matrix, %s, only has %s rows, but attempting to update row %s"
+            , mat, nr, row + 1L
+          )
+        )
+      }
+      
+      col = current_frame[i, "col"]
+      nc = ncol(matrix_list[[mat]])
+      if (is.null(nc)) nc = 1L
+      if (col > nc) {
+        stop(
+          sprintf(
+              "Matrix, %s, only has %s columns, but attempting to update column %s"
+            , mat, nc, col + 1L
+          )
+        )
+      }
+      
+      val = current_frame[i, "current"]
+      if (is.matrix(matrix_list[[mat]])) {
+        matrix_list[[mat]][row + 1L, col + 1L] = val
+      } else {
+        matrix_list[[mat]][row + 1L] = val
+      }
+    } # endfor
+    return(matrix_list)
+  }
+  
+  self$update_params_vector = function(matrix_list) {
+    
+  }
+  
   return_object(self, "OptimizedParams")
 }
 

R/calibrator_arg_constructors.R

@@ -5,15 +5,15 @@
 #' Define the prior distributions for parameters and random effects to be
 #' passed to `par` argument of the \code{\link{mp_tmb_calibrator}} function.
 #' 
-#' @param param Named list of distributional specifications for the
+#' @param params Named list of distributional specifications for the
 #' fixed effects.
 #' @param random Named list of distributional specifications for the random
 #' effects.
 #' @concept create-model-calibrator-args
 #' @export
-mp_par = function(param, random) {
+mp_par = function(params, random) {
   arg = list()
-  arg$param = param
+  arg$params = params
   arg$random = random
   structure(arg, class = "ParArg")
 }
@@ -75,6 +75,14 @@ mp_rbf = function(tv, dimension, initial_weights, seed, prior_sd = 1, fit_prior_
 }
 
 
+mp_piecewise = function(tv, data) {
+  arg = list()
+  arg$tv = tv
+  arg$data = data
+  structure(arg, class = "PiecewiseArg")
+}
+
+
 ## construct objects to pass to the traj argument of mp_tmb_calibrator
 
 #' Trajectory Specification
@@ -90,8 +98,8 @@ mp_rbf = function(tv, dimension, initial_weights, seed, prior_sd = 1, fit_prior_
 #' @concept create-model-calibrator-args
 #' @export
 mp_traj = function(
-      likelihood = list()
-    , condensation = list()
+      likelihood = empty_named_list()
+    , condensation = empty_named_list()
   ) {
   arg = list()
   arg$likelihood = likelihood