<broken> debugging exploding derivatives

R/RcppR6.R

@@ -1,10 +1,6 @@
 ## Generated by RcppR6: do not edit by hand
 ## Version: 0.2.4
-<<<<<<< HEAD
-## Hash: f58c99ad15f3f665237bff0862aa5349
-=======
-## Hash: 584c07ca68ef8e920c7e89821952cc98
->>>>>>> develop
+## Hash: 4804ac8c26f229fdaad40fc82af29dec
 
 ##' @importFrom Rcpp evalCpp
 ##' @importFrom R6 R6Class

R/build_schedule.R

@@ -32,12 +32,15 @@ build_schedule <- function(p, env = make_environment(parameters = p),
   complete = FALSE
 
   # generate coarse initial size distribution
-  if(!is.null(splines))
+  if(!is.null(splines)) {
     state = lapply(splines, partition_spline, n = n_init)
-
+    p$initial_state = unlist(lapply(state, as.vector))
+    p$n_initial_cohorts = unlist(lapply(state, ncol))
+  }
+
   # the refine cohorts
   for (i in seq_len(ctrl$schedule_nsteps)) {
-    res <- run_scm_error(p, env, ctrl, state)
+    res <- run_scm_error(p, env, ctrl)
     net_reproduction_ratios <- res[["net_reproduction_ratios"]]
     split <- lapply(res$err$total, function(x) x > eps)
 

R/scm_support.R

@@ -185,36 +185,6 @@ make_scm <- function(p, env, ctrl, state=NULL) {
   types <- extract_RcppR6_template_types(p, "Parameters")
   scm <- do.call('SCM', types)(p, env, ctrl)
 
-  if(!is.null(state)) {
-    n_str <- length(p$strategies)
-    n_spp = length(state)
-
-    if(n_spp != n_str)
-      stop("State object has more species than strategies defined in Parameters")
-
-    times <- scm$cohort_schedule$all_times
-
-    initial_cohorts <- sapply(times, function(t) sum(t == 0), simplify = F)
-    new_cohorts <- mapply(function(s, i) max(0, ncol(s) - i),
-                          state, initial_cohorts, SIMPLIFY = F)
-
-    new_times <- mapply(function(i, t) c(rep(0, i), t),
-                        new_cohorts, times, SIMPLIFY = F)
-
-    # this introduces one more ind. than necessary, but if we
-    # overwrite the oldest cohorts first then we can just start at t1
-    scm$set_cohort_schedule_times(new_times)
-    scm$run_next()
-
-    # next step starts from first non-zero time
-    start_time <- sapply(times, function(t) min(t[t>0]))
-
-    # add as many new cohorts as required to fit `state` object
-    scm$set_state(min(start_time), unlist(state),
-                  n = mapply(`+`, new_cohorts, initial_cohorts))
-
-  }
-
   return(scm)
 }
 
@@ -239,8 +209,8 @@ scm_patch <- function(i, x) {
 }
 
 run_scm_error <- function(p, env = make_environment(parameters = p),
-                          ctrl = scm_base_control(), state = NULL) {
-  scm <- make_scm(p, env, ctrl, state)
+                          ctrl = scm_base_control()) {
+  scm <- make_scm(p, env, ctrl)
   n_spp <- length(p$strategies)
 
   lai_error <- rep(list(NULL), n_spp)

inst/RcppR6_classes.yml

@@ -360,6 +360,8 @@ Parameters:
     - strategies: "std::vector<T>"
     - birth_rate: "std::vector<double>"
     - is_resident: "std::vector<bool>"
+    - initial_state: "std::vector<double>"
+    - n_initial_cohorts: "std::vector<size_t>"
     - strategy_default: "T"
     - cohort_schedule_times_default: "std::vector<double>"
     - cohort_schedule_times: "std::vector<std::vector<double> >"

inst/include/plant/RcppR6_post.hpp

@@ -554,6 +554,8 @@ template <> inline SEXP wrap(const plant::Parameters<plant::FF16_Strategy,plant:
   ret["strategies"] = Rcpp::wrap(x.strategies);
   ret["birth_rate"] = Rcpp::wrap(x.birth_rate);
   ret["is_resident"] = Rcpp::wrap(x.is_resident);
+  ret["initial_state"] = Rcpp::wrap(x.initial_state);
+  ret["n_initial_cohorts"] = Rcpp::wrap(x.n_initial_cohorts);
   ret["strategy_default"] = Rcpp::wrap(x.strategy_default);
   ret["cohort_schedule_times_default"] = Rcpp::wrap(x.cohort_schedule_times_default);
   ret["cohort_schedule_times"] = Rcpp::wrap(x.cohort_schedule_times);
@@ -584,6 +586,10 @@ template <> inline plant::Parameters<plant::FF16_Strategy,plant::FF16_Environmen
   ret.birth_rate = Rcpp::as<std::vector<double> >(xl["birth_rate"]);
   // ret.is_resident = Rcpp::as<decltype(retis_resident) >(xl["is_resident"]);
   ret.is_resident = Rcpp::as<std::vector<bool> >(xl["is_resident"]);
+  // ret.initial_state = Rcpp::as<decltype(retinitial_state) >(xl["initial_state"]);
+  ret.initial_state = Rcpp::as<std::vector<double> >(xl["initial_state"]);
+  // ret.n_initial_cohorts = Rcpp::as<decltype(retn_initial_cohorts) >(xl["n_initial_cohorts"]);
+  ret.n_initial_cohorts = Rcpp::as<std::vector<size_t> >(xl["n_initial_cohorts"]);
   // ret.strategy_default = Rcpp::as<decltype(retstrategy_default) >(xl["strategy_default"]);
   ret.strategy_default = Rcpp::as<plant::FF16_Strategy >(xl["strategy_default"]);
   // ret.cohort_schedule_times_default = Rcpp::as<decltype(retcohort_schedule_times_default) >(xl["cohort_schedule_times_default"]);
@@ -605,6 +611,8 @@ template <> inline SEXP wrap(const plant::Parameters<plant::FF16w_Strategy,plant
   ret["strategies"] = Rcpp::wrap(x.strategies);
   ret["birth_rate"] = Rcpp::wrap(x.birth_rate);
   ret["is_resident"] = Rcpp::wrap(x.is_resident);
+  ret["initial_state"] = Rcpp::wrap(x.initial_state);
+  ret["n_initial_cohorts"] = Rcpp::wrap(x.n_initial_cohorts);
   ret["strategy_default"] = Rcpp::wrap(x.strategy_default);
   ret["cohort_schedule_times_default"] = Rcpp::wrap(x.cohort_schedule_times_default);
   ret["cohort_schedule_times"] = Rcpp::wrap(x.cohort_schedule_times);
@@ -635,6 +643,10 @@ template <> inline plant::Parameters<plant::FF16w_Strategy,plant::FF16_Environme
   ret.birth_rate = Rcpp::as<std::vector<double> >(xl["birth_rate"]);
   // ret.is_resident = Rcpp::as<decltype(retis_resident) >(xl["is_resident"]);
   ret.is_resident = Rcpp::as<std::vector<bool> >(xl["is_resident"]);
+  // ret.initial_state = Rcpp::as<decltype(retinitial_state) >(xl["initial_state"]);
+  ret.initial_state = Rcpp::as<std::vector<double> >(xl["initial_state"]);
+  // ret.n_initial_cohorts = Rcpp::as<decltype(retn_initial_cohorts) >(xl["n_initial_cohorts"]);
+  ret.n_initial_cohorts = Rcpp::as<std::vector<size_t> >(xl["n_initial_cohorts"]);
   // ret.strategy_default = Rcpp::as<decltype(retstrategy_default) >(xl["strategy_default"]);
   ret.strategy_default = Rcpp::as<plant::FF16w_Strategy >(xl["strategy_default"]);
   // ret.cohort_schedule_times_default = Rcpp::as<decltype(retcohort_schedule_times_default) >(xl["cohort_schedule_times_default"]);
@@ -656,6 +668,8 @@ template <> inline SEXP wrap(const plant::Parameters<plant::FF16r_Strategy,plant
   ret["strategies"] = Rcpp::wrap(x.strategies);
   ret["birth_rate"] = Rcpp::wrap(x.birth_rate);
   ret["is_resident"] = Rcpp::wrap(x.is_resident);
+  ret["initial_state"] = Rcpp::wrap(x.initial_state);
+  ret["n_initial_cohorts"] = Rcpp::wrap(x.n_initial_cohorts);
   ret["strategy_default"] = Rcpp::wrap(x.strategy_default);
   ret["cohort_schedule_times_default"] = Rcpp::wrap(x.cohort_schedule_times_default);
   ret["cohort_schedule_times"] = Rcpp::wrap(x.cohort_schedule_times);
@@ -686,6 +700,10 @@ template <> inline plant::Parameters<plant::FF16r_Strategy,plant::FF16_Environme
   ret.birth_rate = Rcpp::as<std::vector<double> >(xl["birth_rate"]);
   // ret.is_resident = Rcpp::as<decltype(retis_resident) >(xl["is_resident"]);
   ret.is_resident = Rcpp::as<std::vector<bool> >(xl["is_resident"]);
+  // ret.initial_state = Rcpp::as<decltype(retinitial_state) >(xl["initial_state"]);
+  ret.initial_state = Rcpp::as<std::vector<double> >(xl["initial_state"]);
+  // ret.n_initial_cohorts = Rcpp::as<decltype(retn_initial_cohorts) >(xl["n_initial_cohorts"]);
+  ret.n_initial_cohorts = Rcpp::as<std::vector<size_t> >(xl["n_initial_cohorts"]);
   // ret.strategy_default = Rcpp::as<decltype(retstrategy_default) >(xl["strategy_default"]);
   ret.strategy_default = Rcpp::as<plant::FF16r_Strategy >(xl["strategy_default"]);
   // ret.cohort_schedule_times_default = Rcpp::as<decltype(retcohort_schedule_times_default) >(xl["cohort_schedule_times_default"]);
@@ -707,6 +725,8 @@ template <> inline SEXP wrap(const plant::Parameters<plant::K93_Strategy,plant::
   ret["strategies"] = Rcpp::wrap(x.strategies);
   ret["birth_rate"] = Rcpp::wrap(x.birth_rate);
   ret["is_resident"] = Rcpp::wrap(x.is_resident);
+  ret["initial_state"] = Rcpp::wrap(x.initial_state);
+  ret["n_initial_cohorts"] = Rcpp::wrap(x.n_initial_cohorts);
   ret["strategy_default"] = Rcpp::wrap(x.strategy_default);
   ret["cohort_schedule_times_default"] = Rcpp::wrap(x.cohort_schedule_times_default);
   ret["cohort_schedule_times"] = Rcpp::wrap(x.cohort_schedule_times);
@@ -737,6 +757,10 @@ template <> inline plant::Parameters<plant::K93_Strategy,plant::K93_Environment>
   ret.birth_rate = Rcpp::as<std::vector<double> >(xl["birth_rate"]);
   // ret.is_resident = Rcpp::as<decltype(retis_resident) >(xl["is_resident"]);
   ret.is_resident = Rcpp::as<std::vector<bool> >(xl["is_resident"]);
+  // ret.initial_state = Rcpp::as<decltype(retinitial_state) >(xl["initial_state"]);
+  ret.initial_state = Rcpp::as<std::vector<double> >(xl["initial_state"]);
+  // ret.n_initial_cohorts = Rcpp::as<decltype(retn_initial_cohorts) >(xl["n_initial_cohorts"]);
+  ret.n_initial_cohorts = Rcpp::as<std::vector<size_t> >(xl["n_initial_cohorts"]);
   // ret.strategy_default = Rcpp::as<decltype(retstrategy_default) >(xl["strategy_default"]);
   ret.strategy_default = Rcpp::as<plant::K93_Strategy >(xl["strategy_default"]);
   // ret.cohort_schedule_times_default = Rcpp::as<decltype(retcohort_schedule_times_default) >(xl["cohort_schedule_times_default"]);

inst/include/plant/ode_solver.h

@@ -147,6 +147,7 @@ void Solver<System>::step(System& system) {
   const size_t size = y.size();
 
   // Compute the derivatives at the beginning.
+  std::cout << "Computing derivatives" << std::endl;
   setup_dydt_in(system);
 
   while (true) {
@@ -157,9 +158,8 @@ void Solver<System>::step(System& system) {
     }
     stepper.step(system, time, step_size, y, yerr, dydt_in, dydt_out);
 
-    const double step_size_next =
-      control.adjust_step_size(size, stepper.order(), step_size,
-			       y, yerr, dydt_out);
+    const double step_size_next = control.adjust_step_size(
+        size, stepper.order(), step_size, y, yerr, dydt_out);
 
     if (control.step_size_shrank()) {
       // GSL checks that the step size is actually decreased.
@@ -168,15 +168,15 @@ void Solver<System>::step(System& system) {
       // hmin << t
       const double time_next = time + step_size_next;
       if (step_size_next < step_size && time_next > time_orig) {
-	// Step was decreased. Undo step (resetting the state y and
-	// time), and try again with the new step_size.
-	y         = y_orig;
-	time      = time_orig;
-	step_size = step_size_next;
+        // Step was decreased. Undo step (resetting the state y and
+        // time), and try again with the new step_size.
+        y = y_orig;
+        time = time_orig;
+        step_size = step_size_next;
       } else {
-	// We've reached limits of machine accuracy in differences of
-	// step sizes or time (or both).
-	util::stop("Cannot achive the desired accuracy");
+        // We've reached limits of machine accuracy in differences of
+        // step sizes or time (or both).
+        util::stop("Cannot achive the desired accuracy");
       }
     } else {
       // We have successfully taken a step and will return.  Update
@@ -187,10 +187,10 @@ void Solver<System>::step(System& system) {
       //  suggested in the final step, because that step can be very
       //  small compared to previous step, to reach time_max.
       if (final_step) {
-	time = time_max;
+        time = time_max;
       } else {
-	time += step_size;
-	step_size_last = step_size_next;
+        time += step_size;
+        step_size_last = step_size_next;
       }
       prev_times.push_back(time);
       save_dydt_out_as_in();

inst/include/plant/ode_step.h

@@ -14,14 +14,11 @@ class Step {
 public:
   void resize(size_t size_);
   size_t order() const;
-  void step(System& system,
-            double time, double step_size,
-	    state_type &y,
-	    state_type &yerr,
-	    const state_type &dydt_in,
-	    state_type &dydt_out);
-  void derivs(System& system,
-              const state_type& y, state_type& dydt, double t) {
+  void step(System &system, double time, double step_size, state_type &y,
+            state_type &yerr, const state_type &dydt_in, state_type &dydt_out);
+
+  void derivs(System &system, const state_type &y, state_type &dydt, double t) {
+    std::cout << "Derivs. at solver time: " << t << std::endl;
     return ode::derivs(system, y, dydt, t);
   }
 
@@ -80,6 +77,7 @@ void Step<System>::step(System& system,
                         state_type &dydt_out) {
   const double h = step_size; // Historical reasons.
 
+  std::cout << "STEPPING: " << time << std::endl;
   // k1 step:
   std::copy(dydt_in.begin(), dydt_in.end(), k1.begin());
   for (size_t i = 0; i < size; ++i) {
@@ -128,6 +126,8 @@ void Step<System>::step(System& system,
     yerr[i] = h * (ec[1] * k1[i] + ec[3] * k3[i] + ec[4] * k4[i] +
 		   ec[5] * k5[i] + ec[6] * k6[i]);
   }
+
+  std::cout << "FINISHED STEP" << std::endl << std::endl;
 }
 
 // RKCK coefficients, from GSL

inst/include/plant/parameters.h

@@ -56,6 +56,10 @@ struct Parameters {
   std::vector<std::vector<double> > cohort_schedule_times;
   std::vector<double> cohort_schedule_ode_times;
 
+  // Initial size distribution
+  std::vector<double> initial_state;
+  std::vector<size_t> n_initial_cohorts;
+
   // Some little query functions for use on the C side:
   size_t size() const;
   size_t n_residents() const;