Commit ⋅ mrc-ide/SIMPLEGEN

  // objects for storing results
  daily_values = vector<vector<vector<double>>>(params->n_demes, vector<vector<double>>(params->max_time));
  // age distributions. Final level: 0 = Sh, 1 = Eh, 2 = Ah, 3 = Ch, 4 = Ph
  age_distributions = vector<vector<vector<vector<double>>>>(params->n_output_age_times, vector<vector<vector<double>>>(params->n_demes, vector<vector<double>>(params->n_life_table, vector<double>(5))));
  age_distributions = vector<vector<vector<vector<double>>>>(params->n_output_age_times,
                              vector<vector<vector<double>>>(params->n_demes,
                                      vector<vector<double>>(params->n_life_table, vector<double>(8))));
  
  // misc
  EIR = vector<double>(params->n_demes);
  prob_infectious_bite = vector<double>(params->n_demes);
  inc_infection = vector<double>(params->n_demes);
  inc_acute = vector<double>(params->n_demes);
  inc_chronic = vector<double>(params->n_demes);
  
}


@@ -192,7 +198,8 @@

      // get number of new infectious bites on humans
      EIR[k] = params->a * Iv[k] / double(H[k]);
     
      double prob_infectious_bite = 1 - exp(-EIR[k]);  // probability of new infectious bite on host
      // probability of new infectious bite on host
      prob_infectious_bite[k] = 1 - exp(-EIR[k]);
      
      // one method of drawing infections in humans would be to draw the total
      // number of infectious bites from Binomial(H[k], prob_infectious_bite),

@@ -206,17 +213,17 @@

      // host-specific and changes over inoculations. Therefore, as a
      // middleground, draw from Binomial(H[k],
      // max_prob_infection*prob_infectious_bite), where max_prob_infection is
      // the largest value that prob_infection could possibly take. Loop
      // the largest value that prob_infection could possibly take. Then loop
      // through these query infectious bites and draw from a Bernoulli
      // distribution relative to this value.
      //
      // For example, if prob_infection = {0.1, 0.05} then filter based on the
      // value 0.1, i.e. draw the number of query infected hosts from
      // Binomial(H[k], 0.1). Then loop these query hosts and draw from the
      // relative probability of infection, which is Bernoulli with
      // probability {0.1, 0.05}/0.1 = {1.0, 0.5}.
      // Binomial(H[k], 0.1*prob_infectious_bite). Then loop through these query
      // hosts and draw from the relative probability of infection, which is
      // Bernoulli with probability {0.1, 0.05}/0.1 = {1.0, 0.5}.
      
      int host_query_infection = rbinom1(H[k], params->max_prob_infection * prob_infectious_bite);
      int host_query_infection = rbinom1(H[k], params->max_prob_infection * prob_infectious_bite[k]);
      for (int i = 0; i < host_query_infection; ++i) {
        
        // choose host at random

@@ -344,11 +351,14 @@

    // update counts of each host status in each deme
    update_host_counts(t);
    
    // update incidence measures
    update_incidence(t);
    
    // store daily values
    for (int k = 0; k < params->n_demes; ++k) {
      daily_values[k][t] = {double(H[k]), double(Sh[k]), double(Eh[k]), double(Ah[k]), double(Ch[k]), double(Ph[k]),
                            double(Sv[k]), double(Ev[k]), double(Iv[k]),
                            EIR[k],
                            EIR[k], inc_infection[k], inc_acute[k], inc_chronic[k],
                            Ah_detectable_microscopy[k], Ch_detectable_microscopy[k],
                            Ah_detectable_PCR[k], Ch_detectable_PCR[k],n_inoc[k]};
    }

@@ -356,6 +366,7 @@

    // store age distributions
    if (params->output_age_distributions && (params->output_age_times[index_age_distributions] == t+1)) {
      
      // get age distribution
      get_age_distribution(index_age_distributions);
      
      // update index

@@ -371,9 +382,10 @@

        // get sampling parameters
        int this_deme = params->ss_deme[index_obtain_samples];
        int this_n = params->ss_n[index_obtain_samples];
        Diagnosis this_diagnosis = params->ss_diagnosis[index_obtain_samples];
        
        // obtain samples
        get_sample_details(t, this_deme, this_n);
        get_sample_details(t, this_deme, this_n, this_diagnosis);
        
        // increment index
        index_obtain_samples++;

@@ -447,12 +459,41 @@

    H[k] = Sh[k] + Eh[k] + Ah[k] + Ch[k] + Ph[k];
  }
  
}

//------------------------------------------------
// update incidence measures
void Dispatcher::update_incidence(int t) {
  
  // reset values
  fill(inc_infection.begin(), inc_infection.end(), 0.0);
  fill(inc_acute.begin(), inc_acute.end(), 0.0);
  fill(inc_chronic.begin(), inc_chronic.end(), 0.0);
  
  // loop through all hosts, update incidence in given deme
  for (unsigned int i = 0; i < host_pop.size(); ++i) {
    int this_deme = host_pop[i].deme;
    
    // incidence of infection
    double inc1 = prob_infectious_bite[this_deme] * host_pop[i].get_prob_infection() / double(H[this_deme]);
    inc_infection[this_deme] += inc1;
    
    // incidence of acute and chronic infection
    double inc2 = host_pop[i].get_prob_acute();
    inc_acute[this_deme] += inc1 * inc2;
    inc_chronic[this_deme] += inc1 * (1.0 - inc2);
  }
  
}

//------------------------------------------------
// draw sample from deme
void Dispatcher::get_sample_details(int t, int deme, int n) {
void Dispatcher::get_sample_details(int t, int deme, int n, Diagnosis diag) {
  
  // n cannot exceed human population size at this point in time
  if (n > H[deme]) {
    n = H[deme];
  }
  
  // draw vector by sampling without replacement
  vector<int> samp = sample4(n, 0, H[deme] - 1);

@@ -465,8 +506,31 @@

    int this_host_ID = host_pop[this_index].host_ID;
    
    // find if positive for malaria parasites
    bool test_positive =  false;
    Status_host this_host_status = host_pop[this_index].get_host_status();
    bool test_positive = (this_host_status == Host_Ah || this_host_status == Host_Ch);
    if (this_host_status == Host_Ah) {
      
      // get positive prob for acute microscopy vs PCR
      double prob_positive = 0.0;
      if (diag == microscopy) {
        prob_positive = host_pop[this_index].get_detectability_microscopy_acute(t);
      } else if (diag == PCR) {
        prob_positive = host_pop[this_index].get_detectability_PCR_acute(t);
      }
      test_positive = rbernoulli1(prob_positive);
      
    } else if (this_host_status == Host_Ch) {
      
      // get positive prob for chronic microscopy vs PCR
      double prob_positive = 0.0;
      if (diag == microscopy) {
        prob_positive = host_pop[this_index].get_detectability_microscopy_chronic(t);
      } else if (diag == PCR) {
        prob_positive = host_pop[this_index].get_detectability_PCR_chronic(t);
      }
      test_positive = rbernoulli1(prob_positive);
      
    }
    
    // save basic details
    vector<int> this_details = {t+1, deme+1, this_host_ID, test_positive};

@@ -495,7 +559,7 @@

  // get current time
  int t = params->output_age_times[t_index];
  
  // loop through all hosts, update age distribution
  // loop through all hosts, update age distributions
  for (unsigned int i = 0; i < host_pop.size(); ++i) {
    int this_deme = host_pop[i].deme;
    int this_age = host_pop[i].get_age(t);

@@ -519,6 +583,16 @@

    default:
      Rcpp::stop("invalid host status in get_age_distribution()");
    }
  }
    
    // incidence of infection
    double inc1 = prob_infectious_bite[this_deme] * host_pop[i].get_prob_infection() / double(H[this_deme]);
    age_distributions[t_index][this_deme][this_age][5] += inc1;
    
    // incidence of acute and chronic infection
    double inc2 = host_pop[i].get_prob_acute();
    age_distributions[t_index][this_deme][this_age][6] += inc1 * inc2;
    age_distributions[t_index][this_deme][this_age][7] += inc1 * (1.0 - inc2);
    
  }  // end loop through hosts
  
}

@@ -1,653 +1,653 @@


#### HELPER FUNCTIONS ####################################################################

#------------------------------------------------
# for single value, return value as string. For vector of values return string
# of comma-separated values enclosed in curly brackets
#' @noRd
nice_format <- function(x) {
  if (is.null(x)) {
    return("")
  }
  if (length(x)==1) {
    ret <- as.character(x)
  } else {
    ret <- paste0("{", paste(x, collapse = ", "), "}")
  }
  return(ret)
}

#### BASIC OBJECT TYPES ####################################################################

#------------------------------------------------
# x is NULL
#' @noRd
assert_null <- function(x, message = "%s must be null",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.null(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is not NULL
#' @noRd
assert_non_null <- function(x, message = "%s cannot be null",
                            name = paste(deparse(substitute(x)), collapse = "")) {
  if (is.null(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is atomic
#' @noRd
assert_atomic <- function(x, message = "%s must be atomic (see ?is.atomic)",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.atomic(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is atomic and single valued (has length 1)
#' @noRd
assert_single <- function(x, message = "%s must be a single value",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_atomic(x, name = name)
  assert_length(x, 1, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is character string
#' @noRd
assert_string <- function(x, message = "%s must be character string",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.character(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single character string
#' @noRd
assert_single_string <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_string(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is logical
#' @noRd
assert_logical <- function(x, message = "%s must be logical",
                           name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.logical(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single logical
#' @noRd
assert_single_logical <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_logical(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is numeric
#' @noRd
assert_numeric <- function(x, message = "%s must be numeric",
                           name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.numeric(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single numeric
#' @noRd
assert_single_numeric <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_numeric(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is integer
#' @noRd
assert_int <- function(x, message = "%s must be integer valued",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (!isTRUE(all.equal(x, as.integer(x), check.attributes = FALSE))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single integer
#' @noRd
assert_single_int <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_int(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is positive (with or without zero allowed)
#' @noRd
assert_pos <- function(x, zero_allowed = TRUE, message1 = "%s must be greater than or equal to zero",
                       message2 = "%s must be greater than zero",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (zero_allowed) {
    if (!all(x>=0)) {
      stop(sprintf(message1, name), call. = FALSE)
    }
  } else {
    if (!all(x>0)) {
      stop(sprintf(message2, name), call. = FALSE)
    }
  }
  return(TRUE)
}

#------------------------------------------------
# x is single positive (with or without zero allowed)
#' @noRd
assert_single_pos <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_pos(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is positive integer (with or without zero allowed)
#' @noRd
assert_pos_int <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_int(x, name = name)
  assert_pos(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is single positive integer (with or without zero allowed)
#' @noRd
assert_single_pos_int <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_pos_int(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is single value bounded between limits
#' @noRd
assert_single_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_bounded(x, left = left, right = right,
                 inclusive_left = inclusive_left, inclusive_right = inclusive_right,
                 name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector (and is not a list or another recursive type)
#' @noRd
assert_vector <- function(x, message = "%s must be a non-recursive vector",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.vector(x) || is.recursive(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a numeric vector
#' @noRd
assert_vector_numeric <- function(x, message = "%s must be a non-recursive vector of numeric values",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, message = message, name = name)
  assert_vector(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of integers
#' @noRd
assert_vector_int <- function(x, message = "%s must be a non-recursive vector of integers",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_int(x, message = message, name = name)
  assert_vector(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of positive values
#' @noRd
assert_vector_pos <- function(x, zero_allowed = TRUE,
                              message1 = "%s must be a non-recursive vector of positive values or zero",
                              message2 = "%s must be a non-recursive vector of positive values",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos(x, zero_allowed = zero_allowed, message1 = message1, message2 = message2, name = name)
  if (zero_allowed) {
    assert_vector(x, message = message1, name = name)
  } else {
    assert_vector(x, message = message2, name = name)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a vector of positive integers
#' @noRd
assert_vector_pos_int <- function(x, zero_allowed = TRUE,
                                  message1 = "%s must be a non-recursive vector of positive integers or zero",
                                  message2 = "%s must be a non-recursive vector of positive integers",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos(x, zero_allowed = zero_allowed, message1 = message1, message2 = message2, name = name)
  
  if (zero_allowed) {
    message_final = message1
  } else {
    message_final = message2
  }
  assert_vector(x, message = message_final, name = name)
  assert_int(x, message = message_final, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of bounded values
#' @noRd
assert_vector_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_vector(x, name = name)
  assert_bounded(x, left = left, right = right,
                 inclusive_left = inclusive_left, inclusive_right = inclusive_right,
                 name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a matrix
#' @noRd
assert_matrix <- function(x, message = "%s must be a matrix",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.matrix(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a matrix of numeric values
#' @noRd
assert_matrix_numeric <- function(x, message = "%s must be a numeric matrix",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_matrix(x, message = message, name = name)
  assert_numeric(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a list
assert_list <- function(x, message = "%s must be a list",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.list(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a data frame
assert_dataframe <- function(x, message = "%s must be a data frame",
                             name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.data.frame(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x inherits from custom class c
assert_custom_class <- function(x, c, message = "%s must inherit from class '%s'",
                                name = paste(deparse(substitute(x)), collapse = "")) {
  if (!inherits(x, c)) {
    stop(sprintf(message, name, c), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a plotting limit, i.e. contains two increasing values
assert_limit <- function(x, message = "%s must be a valid plotting limit, i.e. contain two increasing values",
                         name = paste(deparse(substitute(x)), collapse = "")) {
  assert_vector(x, name = name)
  assert_length(x, 2, name = name)
  assert_numeric(x, name = name)
  assert_increasing(x, name = name)
  return(TRUE)
}


#### VALUE COMPARISONS ####################################################################

#------------------------------------------------
# x and y are equal in all matched comparisons. x and y can be any type
#' @noRd
assert_eq <- function(x, y, message = "%s must equal %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  assert_same_length(x, y, name_x = name_x, name_y = name_y)
  if (!isTRUE(all.equal(x, y, check.attributes = FALSE))) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x and y are unequal in all matched comparisons. x and y can be any type
#' @noRd
assert_neq <- function(x, y, message = "%s cannot equal %s",
                       name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  assert_same_length(x, y, name_x = name_x, name_y = name_y)
  if (any(x == y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is greater than y in all matched comparisons
#' @noRd
assert_gr <- function(x, y, message = "%s must be greater than %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""),
                      name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x>y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is greater than or equal to y in all matched comparisons
#' @noRd
assert_greq <- function(x, y, message = "%s must be greater than or equal to %s",
                        name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x >= y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is less than y in all matched comparisons
#' @noRd
assert_le <- function(x, y, message = "%s must be less than %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x<y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is less than or equal to y in all matched comparisons
#' @noRd
assert_leq <- function(x, y, message = "%s must be less than or equal to %s",
                       name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x<=y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is between bounds (inclusive or exclusive)
#' @noRd
assert_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                           name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (inclusive_left) {
    if (!all(x>=left)) {
      stop(sprintf("%s must be greater than or equal to %s", name, left), call. = FALSE)
    }
  } else {
    if (!all(x>left)) {
      stop(sprintf("%s must be greater than %s", name, left), call. = FALSE)
    }
  }
  if (inclusive_right) {
    if (!all(x<=right)) {
      stop(sprintf("%s must be less than or equal to %s", name, right), call. = FALSE)
    }
  } else {
    if (!all(x<right)) {
      stop(sprintf("%s must be less than %s", name, right), call. = FALSE)
    }
  }
  return(TRUE)
}

#------------------------------------------------
# all x are in y
#' @noRd
assert_in <- function(x, y, message = "all %s must be in %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  if (!all(x %in% y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# none of x are in y
#' @noRd
assert_not_in <- function(x, y, message = "none of %s can be in %s",
                          name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  if (any(x %in% y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}


#### DIMENSIONS ####################################################################

#------------------------------------------------
# length(x) equals n
#' @noRd
assert_length <- function(x, n, message = "%s must be of length %s",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos_int(n)
  if (length(x) != n[1]) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x and y are same length
#' @noRd
assert_same_length <- function(x, y, message =  "%s and %s must be the same length",
                               name_x = paste(deparse(substitute(x)), collapse = ""),
                               name_y = paste(deparse(substitute(y)))) {
  if (length(x) != length(y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# multiple objects all same length
#' @noRd
assert_same_length_multiple <- function(...) {
  l <- mapply(length, list(...))
  if (length(unique(l)) != 1) {
    l_names <- sapply(match.call(expand.dots = FALSE)$..., deparse)
    stop(sprintf("variables %s must be the same length", nice_format(l_names)), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is two-dimensional
#' @noRd
assert_2d <- function(x, message = "%s must be two-dimensional",
                      name = paste(deparse(substitute(x)), collapse = "")) {
  is_2d <- FALSE
  if (!is.null(dim(x))) {
    if (length(dim(x)) == 2) {
      is_2d <- TRUE
    }
  }
  if (!is_2d) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# nrow(x) equals n
#' @noRd
assert_nrow <- function(x, n, message = "%s must have %s rows",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  if (nrow(x) != n) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# ncol(x) equals n
#' @noRd
assert_ncol <- function(x, n, message = "%s must have %s cols",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  if (ncol(x) != n) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# dim(x) equals y
#' @noRd
assert_dim <- function(x, y, message = "%s must have %s rows and %s columns",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  assert_pos_int(y, name = "y variable in assert_dim()")
  assert_length(y, 2, name = "y variable in assert_dim()")
  if (nrow(x) != y[1] | ncol(x) != y[2]) {
    stop(sprintf(message, name, y[1], y[2]), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is square matrix
#' @noRd
assert_square_matrix <- function(x, message = "%s must be a square matrix",
                                 name = paste(deparse(substitute(x)), collapse = "")) {
  assert_matrix(x, name = name)
  if (nrow(x) != ncol(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# is symmetric matrix
#' @noRd
assert_symmetric_matrix <- function(x, message = "%s must be a symmetric matrix",
                                    name = paste(deparse(substitute(x)), collapse = "")) {
  assert_square_matrix(x, name = name)
  if (!isSymmetric(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#### MISC ####################################################################

#------------------------------------------------
# x contains no duplicates
#' @noRd
assert_noduplicates <- function(x, message = "%s must contain no duplicates",
                                name = paste(deparse(substitute(x)), collapse = "")) {
  if (any(duplicated(x))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# file exists at chosen path
#' @noRd
assert_file_exists <- function(x, message = "file not found at path %s",
                               name = paste(deparse(substitute(x)), collapse = "")) {
  if (!file.exists(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is increasing
#' @noRd
assert_increasing <- function(x, message = "%s must be increasing",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_numeric(x, name = name)
  if (!all.equal(x, sort(x))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is decreasing
#' @noRd
assert_decreasing <- function(x, message = "%s must be decreasing",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_numeric(x, name = name)
  if (!all.equal(x, sort(x, decreasing = TRUE))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}


#### HELPER FUNCTIONS ####################################################################

#------------------------------------------------
# for single value, return value as string. For vector of values return string
# of comma-separated values enclosed in curly brackets
#' @noRd
nice_format <- function(x) {
  if (is.null(x)) {
    return("")
  }
  if (length(x)==1) {
    ret <- as.character(x)
  } else {
    ret <- paste0("{", paste(x, collapse = ", "), "}")
  }
  return(ret)
}

#### BASIC OBJECT TYPES ####################################################################

#------------------------------------------------
# x is NULL
#' @noRd
assert_null <- function(x, message = "%s must be null",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.null(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is not NULL
#' @noRd
assert_non_null <- function(x, message = "%s cannot be null",
                            name = paste(deparse(substitute(x)), collapse = "")) {
  if (is.null(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is atomic
#' @noRd
assert_atomic <- function(x, message = "%s must be atomic (see ?is.atomic)",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.atomic(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is atomic and single valued (has length 1)
#' @noRd
assert_single <- function(x, message = "%s must be a single value",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_atomic(x, name = name)
  assert_length(x, 1, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is character string
#' @noRd
assert_string <- function(x, message = "%s must be character string",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.character(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single character string
#' @noRd
assert_single_string <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_string(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is logical
#' @noRd
assert_logical <- function(x, message = "%s must be logical",
                           name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.logical(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single logical
#' @noRd
assert_single_logical <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_logical(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is numeric
#' @noRd
assert_numeric <- function(x, message = "%s must be numeric",
                           name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.numeric(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single numeric
#' @noRd
assert_single_numeric <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_numeric(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is integer
#' @noRd
assert_int <- function(x, message = "%s must be integer valued",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (!isTRUE(all.equal(x, as.integer(x), check.attributes = FALSE))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is single integer
#' @noRd
assert_single_int <- function(x, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_int(x, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is positive (with or without zero allowed)
#' @noRd
assert_pos <- function(x, zero_allowed = TRUE, message1 = "%s must be greater than or equal to zero",
                       message2 = "%s must be greater than zero",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (zero_allowed) {
    if (!all(x>=0)) {
      stop(sprintf(message1, name), call. = FALSE)
    }
  } else {
    if (!all(x>0)) {
      stop(sprintf(message2, name), call. = FALSE)
    }
  }
  return(TRUE)
}

#------------------------------------------------
# x is single positive (with or without zero allowed)
#' @noRd
assert_single_pos <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_pos(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is positive integer (with or without zero allowed)
#' @noRd
assert_pos_int <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_int(x, name = name)
  assert_pos(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is single positive integer (with or without zero allowed)
#' @noRd
assert_single_pos_int <- function(x, zero_allowed = TRUE, name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_pos_int(x, zero_allowed = zero_allowed, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is single value bounded between limits
#' @noRd
assert_single_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_length(x, n = 1, name = name)
  assert_bounded(x, left = left, right = right,
                 inclusive_left = inclusive_left, inclusive_right = inclusive_right,
                 name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector (and is not a list or another recursive type)
#' @noRd
assert_vector <- function(x, message = "%s must be a non-recursive vector",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.vector(x) || is.recursive(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a numeric vector
#' @noRd
assert_vector_numeric <- function(x, message = "%s must be a non-recursive vector of numeric values",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, message = message, name = name)
  assert_vector(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of integers
#' @noRd
assert_vector_int <- function(x, message = "%s must be a non-recursive vector of integers",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_int(x, message = message, name = name)
  assert_vector(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of positive values
#' @noRd
assert_vector_pos <- function(x, zero_allowed = TRUE,
                              message1 = "%s must be a non-recursive vector of positive values or zero",
                              message2 = "%s must be a non-recursive vector of positive values",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos(x, zero_allowed = zero_allowed, message1 = message1, message2 = message2, name = name)
  if (zero_allowed) {
    assert_vector(x, message = message1, name = name)
  } else {
    assert_vector(x, message = message2, name = name)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a vector of positive integers
#' @noRd
assert_vector_pos_int <- function(x, zero_allowed = TRUE,
                                  message1 = "%s must be a non-recursive vector of positive integers or zero",
                                  message2 = "%s must be a non-recursive vector of positive integers",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos(x, zero_allowed = zero_allowed, message1 = message1, message2 = message2, name = name)
  
  if (zero_allowed) {
    message_final = message1
  } else {
    message_final = message2
  }
  assert_vector(x, message = message_final, name = name)
  assert_int(x, message = message_final, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a vector of bounded values
#' @noRd
assert_vector_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_vector(x, name = name)
  assert_bounded(x, left = left, right = right,
                 inclusive_left = inclusive_left, inclusive_right = inclusive_right,
                 name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a matrix
#' @noRd
assert_matrix <- function(x, message = "%s must be a matrix",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.matrix(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a matrix of numeric values
#' @noRd
assert_matrix_numeric <- function(x, message = "%s must be a numeric matrix",
                                  name = paste(deparse(substitute(x)), collapse = "")) {
  assert_matrix(x, message = message, name = name)
  assert_numeric(x, message = message, name = name)
  return(TRUE)
}

#------------------------------------------------
# x is a list
assert_list <- function(x, message = "%s must be a list",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.list(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a data frame
assert_dataframe <- function(x, message = "%s must be a data frame",
                             name = paste(deparse(substitute(x)), collapse = "")) {
  if (!is.data.frame(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x inherits from custom class c
assert_custom_class <- function(x, c, message = "%s must inherit from class '%s'",
                                name = paste(deparse(substitute(x)), collapse = "")) {
  if (!inherits(x, c)) {
    stop(sprintf(message, name, c), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is a plotting limit, i.e. contains two increasing values
assert_limit <- function(x, message = "%s must be a valid plotting limit, i.e. contain two increasing values",
                         name = paste(deparse(substitute(x)), collapse = "")) {
  assert_vector(x, name = name)
  assert_length(x, 2, name = name)
  assert_numeric(x, name = name)
  assert_increasing(x, name = name)
  return(TRUE)
}


#### VALUE COMPARISONS ####################################################################

#------------------------------------------------
# x and y are equal in all matched comparisons. x and y can be any type
#' @noRd
assert_eq <- function(x, y, message = "%s must equal %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  assert_same_length(x, y, name_x = name_x, name_y = name_y)
  if (!isTRUE(all.equal(x, y, check.attributes = FALSE))) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x and y are unequal in all matched comparisons. x and y can be any type
#' @noRd
assert_neq <- function(x, y, message = "%s cannot equal %s",
                       name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  assert_same_length(x, y, name_x = name_x, name_y = name_y)
  if (any(x == y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is greater than y in all matched comparisons
#' @noRd
assert_gr <- function(x, y, message = "%s must be greater than %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""),
                      name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x>y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is greater than or equal to y in all matched comparisons
#' @noRd
assert_greq <- function(x, y, message = "%s must be greater than or equal to %s",
                        name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x >= y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is less than y in all matched comparisons
#' @noRd
assert_le <- function(x, y, message = "%s must be less than %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x<y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is less than or equal to y in all matched comparisons
#' @noRd
assert_leq <- function(x, y, message = "%s must be less than or equal to %s",
                       name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_numeric(x, name = name_x)
  assert_numeric(y, name = name_y)
  assert_in(length(y), c(1,length(x)))
  if (!all(x<=y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is between bounds (inclusive or exclusive)
#' @noRd
assert_bounded <- function(x, left = 0, right = 1, inclusive_left = TRUE, inclusive_right = TRUE,
                           name = paste(deparse(substitute(x)), collapse = "")) {
  assert_numeric(x, name = name)
  if (inclusive_left) {
    if (!all(x>=left)) {
      stop(sprintf("%s must be greater than or equal to %s", name, left), call. = FALSE)
    }
  } else {
    if (!all(x>left)) {
      stop(sprintf("%s must be greater than %s", name, left), call. = FALSE)
    }
  }
  if (inclusive_right) {
    if (!all(x<=right)) {
      stop(sprintf("%s must be less than or equal to %s", name, right), call. = FALSE)
    }
  } else {
    if (!all(x<right)) {
      stop(sprintf("%s must be less than %s", name, right), call. = FALSE)
    }
  }
  return(TRUE)
}

#------------------------------------------------
# all x are in y
#' @noRd
assert_in <- function(x, y, message = "all %s must be in %s",
                      name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  if (!all(x %in% y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# none of x are in y
#' @noRd
assert_not_in <- function(x, y, message = "none of %s can be in %s",
                          name_x = paste(deparse(substitute(x)), collapse = ""), name_y = nice_format(y)) {
  assert_non_null(x, name = name_x)
  assert_non_null(y, name = name_y)
  if (any(x %in% y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}


#### DIMENSIONS ####################################################################

#------------------------------------------------
# length(x) equals n
#' @noRd
assert_length <- function(x, n, message = "%s must be of length %s",
                          name = paste(deparse(substitute(x)), collapse = "")) {
  assert_pos_int(n)
  if (length(x) != n[1]) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x and y are same length
#' @noRd
assert_same_length <- function(x, y, message =  "%s and %s must be the same length",
                               name_x = paste(deparse(substitute(x)), collapse = ""),
                               name_y = paste(deparse(substitute(y)))) {
  if (length(x) != length(y)) {
    stop(sprintf(message, name_x, name_y), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# multiple objects all same length
#' @noRd
assert_same_length_multiple <- function(...) {
  l <- mapply(length, list(...))
  if (length(unique(l)) != 1) {
    l_names <- sapply(match.call(expand.dots = FALSE)$..., deparse)
    stop(sprintf("variables %s must be the same length", nice_format(l_names)), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is two-dimensional
#' @noRd
assert_2d <- function(x, message = "%s must be two-dimensional",
                      name = paste(deparse(substitute(x)), collapse = "")) {
  is_2d <- FALSE
  if (!is.null(dim(x))) {
    if (length(dim(x)) == 2) {
      is_2d <- TRUE
    }
  }
  if (!is_2d) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# nrow(x) equals n
#' @noRd
assert_nrow <- function(x, n, message = "%s must have %s rows",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  if (nrow(x) != n) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# ncol(x) equals n
#' @noRd
assert_ncol <- function(x, n, message = "%s must have %s cols",
                        name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  if (ncol(x) != n) {
    stop(sprintf(message, name, n), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# dim(x) equals y
#' @noRd
assert_dim <- function(x, y, message = "%s must have %s rows and %s columns",
                       name = paste(deparse(substitute(x)), collapse = "")) {
  assert_2d(x, name = name)
  assert_pos_int(y, name = "y variable in assert_dim()")
  assert_length(y, 2, name = "y variable in assert_dim()")
  if (nrow(x) != y[1] | ncol(x) != y[2]) {
    stop(sprintf(message, name, y[1], y[2]), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is square matrix
#' @noRd
assert_square_matrix <- function(x, message = "%s must be a square matrix",
                                 name = paste(deparse(substitute(x)), collapse = "")) {
  assert_matrix(x, name = name)
  if (nrow(x) != ncol(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# is symmetric matrix
#' @noRd
assert_symmetric_matrix <- function(x, message = "%s must be a symmetric matrix",
                                    name = paste(deparse(substitute(x)), collapse = "")) {
  assert_square_matrix(x, name = name)
  if (!isSymmetric(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#### MISC ####################################################################

#------------------------------------------------
# x contains no duplicates
#' @noRd
assert_noduplicates <- function(x, message = "%s must contain no duplicates",
                                name = paste(deparse(substitute(x)), collapse = "")) {
  if (any(duplicated(x))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# file exists at chosen path
#' @noRd
assert_file_exists <- function(x, message = "file not found at path %s",
                               name = paste(deparse(substitute(x)), collapse = "")) {
  if (!file.exists(x)) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is increasing
#' @noRd
assert_increasing <- function(x, message = "%s must be increasing",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_numeric(x, name = name)
  if (!isTRUE(all.equal(x, sort(x)))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}

#------------------------------------------------
# x is decreasing
#' @noRd
assert_decreasing <- function(x, message = "%s must be decreasing",
                              name = paste(deparse(substitute(x)), collapse = "")) {
  assert_non_null(x, name = name)
  assert_numeric(x, name = name)
  if (!isTRUE(all.equal(x, sort(x, decreasing = TRUE)))) {
    stop(sprintf(message, name), call. = FALSE)
  }
  return(TRUE)
}


@@ -378,7 +378,7 @@

#'
#' @param project a SIMPLEGEN project, as produced by the
#'   \code{simplegen_project()} function.
#' @param x a dataframe containing all of the following columns:
#' @param df_sample a dataframe containing all of the following columns:
#'   \itemize{
#'     \item time: the time (in days) at which samples are taken.
#'     \item deme: the deme from which samples are taken.

@@ -393,24 +393,25 @@

#'
#' @export

define_sampling_strategy <- function(project, x) {
define_sampling_strategy <- function(project, df_sample) {
  
  # check inputs
  assert_custom_class(project, "simplegen_project")
  assert_dataframe(x)
  assert_eq(names(x), c("time", "deme", "case_detection", "diagnosis", "n"))
  assert_pos_int(x$time, zero_allowed = FALSE)
  assert_pos_int(x$deme, zero_allowed = FALSE)
  assert_in(x$case_detection, c("active", "passive"))
  assert_in(x$diagnosis, c("microscopy", "PCR"))
  assert_pos_int(x$n, zero_allowed = FALSE)
  assert_dataframe(df_sample)
  assert_eq(names(df_sample), c("time", "deme", "case_detection", "diagnosis", "n"))
  assert_pos_int(df_sample$time, zero_allowed = FALSE)
  assert_increasing(df_sample$time)
  assert_pos_int(df_sample$deme, zero_allowed = FALSE)
  assert_in(df_sample$case_detection, c("active", "passive"))
  assert_in(df_sample$diagnosis, c("microscopy", "PCR"))
  assert_pos_int(df_sample$n, zero_allowed = FALSE)
  
  # specify formats
  x$case_detection <- as.character(x$case_detection)
  x$diagnosis <- as.character(x$diagnosis)
  df_sample$case_detection <- as.character(df_sample$case_detection)
  df_sample$diagnosis <- as.character(df_sample$diagnosis)
  
  # load into project
  project$sampling_strategy <- x
  project$sampling_strategy <- df_sample
  
  invisible(project)
}

@@ -468,8 +469,7 @@

  assert_single_logical(overwrite_transmission_record)
  assert_single_logical(output_daily_counts)
  assert_single_logical(output_age_distributions)
  assert_vector(output_age_times)
  assert_pos_int(output_age_times)
  assert_vector_pos_int(output_age_times, zero_allowed = FALSE)
  assert_leq(output_age_times, max_time)
  assert_single_logical(pb_markdown)
  assert_single_logical(silent)

@@ -486,6 +486,13 @@

    stop("no epi parameters defined. See ?define_epi_params")
  }
  
  # check that inputs are compatible with sampling strategy
  if (!is.null(project$sampling_strategy)) {
    ss <- project$sampling_strategy
    
    assert_greq(max_time, max(ss$time), message = "%s exceeded by sampling time (%s)")
    assert_greq(nrow(project$epi_parameters$mig_mat), max(ss$deme), "number of demes defined in sampling strategy exceeds number defined in epi parameters")
  }
  
  # ---------- define argument lists ----------
  

@@ -559,7 +566,7 @@

      ret <- as.data.frame(cbind(seq_len(nrow(ret)), i, ret))
      names(ret) <- c("time", "deme", "H", "S", "E", "A", "C", "P",
                      "Sv", "Ev", "Iv",
                      "EIR",
                      "EIR", "inc_infection", "inc_acute", "inc_chronic",
                      "A_detectable_microscopy", "C_detectable_microscopy",
                      "A_detectable_PCR", "C_detectable_PCR","n_inoc")
      return(ret)

@@ -572,7 +579,7 @@

    age_distributions <- do.call(rbind, mapply(function(j) {
      ret <- do.call(rbind, mapply(function(i) {
        ret <- do.call(rbind, output_raw$age_distributions[[j]][[i]])
        colnames(ret) <- c("S", "E", "A", "C", "P")
        colnames(ret) <- c("S", "E", "A", "C", "P", "inc_infection", "inc_acute", "inc_chronic")
        data.frame(cbind(deme = i, age = seq_len(nrow(ret)) - 1, ret))
      }, seq_along(output_raw$age_distributions[[j]]), SIMPLIFY = FALSE))
      cbind(sample_time = j, ret)

@@ -260,7 +260,7 @@

  
  # subset data
  dat <- project$epi_output$age_distributions
  dat <- dat[dat$sample_time == sample_time & dat$deme == deme, c("age", state)]
  dat <- dat[(dat$sample_time == sample_time) & (dat$deme == deme), c("age", state)]
  names(dat) <- c("age", "y")
  
  # produce plot

@@ -474,7 +474,7 @@

  for (int i = 0; i < params->max_inoculations; ++i) {
    
    // anything that is currently in the acute or chronic stages is cured
    if (inoc_status_asexual[i] == Acute_asexual || inoc_status_asexual[i] == Chronic_asexual) {
    if ((inoc_status_asexual[i] == Acute_asexual) || (inoc_status_asexual[i] == Chronic_asexual)) {
      
      // if due to become infective at a future timepoint then store this
      // timepoint

@@ -509,7 +509,8 @@

    // prophylactic period is cured immediately upon emergence
    if (inoc_status_asexual[i] == Liverstage_asexual) {
      
      if (inoc_events[i].count(Event_Eh_to_Ah) != 0 && inoc_events[i][Event_Eh_to_Ah] <= t2) {
      // acute
      if ((inoc_events[i].count(Event_Eh_to_Ah) != 0) && (inoc_events[i][Event_Eh_to_Ah] <= t2)) {
        
        // store time at which due to emerge
        int t_emerge = inoc_events[i][Event_Eh_to_Ah];

@@ -522,7 +523,8 @@

        new_inoc_event(t_emerge, Event_Ah_to_Sh, i);
        new_inoc_event(t_emerge, Event_end_infective, i);
        
      } else if (inoc_events[i].count(Event_Eh_to_Ch) != 0 && inoc_events[i][Event_Eh_to_Ch] <= t2) {
      // chronic
      } else if ((inoc_events[i].count(Event_Eh_to_Ch) != 0) && (inoc_events[i][Event_Eh_to_Ch] <= t2)) {
        
        // store time at which due to emerge
        int t_emerge = inoc_events[i][Event_Eh_to_Ch];

@@ -544,6 +546,8 @@

  t_next_inoc_event = params->max_time + 1;
  for (int i = 0; i < params->max_inoculations; ++i) {
    for (const auto & x : inoc_events[i]) {
      
      // catch to ensure that events cannot predate current time
      if (x.second < t) {
        print("time =", t);
        print_inoc_status();

@@ -814,6 +818,13 @@

//------------------------------------------------
// get current probability of infection
double Host::get_prob_infection() {
  
  // situations in which zero chance of infection
  if ((get_n_active_inoc() == params->max_inoculations) || prophylaxis_on) {
    return 0.0;
  }
  
  // get probability from flexible distribution
  int index = (infection_index < params->n_prob_infection) ? infection_index : params->n_prob_infection - 1;
  return params->prob_infection[index];
}
imilarity index 97%
ename from tests/testthat/test-assertions_v8.R
ename to tests/testthat/test-assertions_v9.R

@@ -78,6 +78,10 @@

  
  // misc
  std::vector<double> EIR;
  std::vector<double> prob_infectious_bite;
  std::vector<double> inc_infection;
  std::vector<double> inc_acute;
  std::vector<double> inc_chronic;
  
  // number of active inoculations
  std::vector<double> n_inoc;

@@ -91,7 +95,8 @@

  void init(Parameters &params);
  void run_simulation(Rcpp::List &args_functions, Rcpp::List &args_progress);
  void update_host_counts(int t);
  void get_sample_details(int t, int deme, int n);
  void update_incidence(int t);
  void get_sample_details(int t, int deme, int n, Diagnosis diag);
  void get_age_distribution(int t_index);
  
};

@@ -1,54 +1,81 @@


#------------------------------------------------
#' @title Define new SIMPLEGEN project
#'
#' @description Define a new SIMPLEGEN project. This project will hold all
#'   simulation inputs and outputs for a given analysis, and is initialised with
#'   the default values of all parameters.
#'
#' @export

simplegen_project <- function() {
  
  # create empty project
  project <- list(epi_parameters = NULL,
                  sampling_strategy = NULL,
                  epi_output = NULL,
                  sample_details = NULL,
                  relatedness = NULL,
                  true_genotypes = NULL,
                  observed_genotypes = NULL)
  
  class(project) <- "simplegen_project"
  
  # return
  invisible(project)
}

#------------------------------------------------
# overload print() function for simplegen_project
#' @method print simplegen_project
#' @export
print.simplegen_project <- function(x, ...) {
  
  # print summary
  summary(x)
  
  # return invisibly
  invisible(x)
}

#------------------------------------------------
# overload summary() function for simplegen_project
#' @method summary simplegen_project
#' @export
summary.simplegen_project <- function(object, ...) {
  
  
  
  message("TODO - some default print method")
  
  # return invisibly
  invisible(object)
}


#------------------------------------------------
#' @title Define new SIMPLEGEN project
#'
#' @description Define a new SIMPLEGEN project. This project will hold all
#'   simulation inputs and outputs for a given analysis, and is initialised with
#'   the default values of all parameters.
#'
#' @export

simplegen_project <- function() {
  
  # create empty project
  project <- list(epi_parameters = NULL,
                  sampling_strategy = NULL,
                  epi_output = NULL,
                  sample_details = NULL,
                  relatedness = NULL,
                  true_genotypes = NULL,
                  observed_genotypes = NULL)
  
  class(project) <- "simplegen_project"
  
  # return
  invisible(project)
}

#------------------------------------------------
# overload print() function for simplegen_project
#' @method print simplegen_project
#' @export
print.simplegen_project <- function(x, ...) {
  
  # print summary
  summary(x)
  
  # return invisibly
  invisible(x)
}

#------------------------------------------------
# overload summary() function for simplegen_project
#' @method summary simplegen_project
#' @export
summary.simplegen_project <- function(object, ...) {
  p <- object
  
  # if empty project
  if (all(mapply(is.null, p))) {
    message("(empty project)")
    invisible(object)
  }
  
  # print epi parameters
  if (!is.null(p$epi_parameters)) {
    message("Epidemiological model:")
    n_demes <- length(p$epi_parameters$H)
    message(sprintf("  demes: %s", n_demes))
    message(sprintf("  H:\t %s", paste(p$epi_parameters$H, collapse = ", ")))
    message(sprintf("  M:\t %s", paste(p$epi_parameters$M, collapse = ", ")))
    message(sprintf("  seed infections: %s", paste(p$epi_parameters$seed_infections, collapse = ", ")))
  }
  
  # print sampling strategy
  if (!is.null(p$sampling_strategy)) {
    message("Sampling strategy:")
    n_time <- unique(p$sampling_strategy$time)
    n_samp_time <- mapply(sum, split(p$sampling_strategy$n, f = p$sampling_strategy$time))
    if (length(n_time) <= 5) {
      message(sprintf("  time: %s", paste(n_time, collapse = ", ")))
      message(sprintf("  n: %s", paste(n_samp_time, collapse = ", ")))
    } else {
      message("  time: (more than 5)")
      message("  n: (more than 5)")
    }
    
  }
  
  invisible(object)
}


Files		•	•	Coverage
R	613	398	215	64.93%
src	1,397	741	656	53.04%
Project Totals (24 files)	2,010	1,139	871	56.67%

View CI Build Download

comment: false

coverage:
  status:
    project:
      default:
        target: auto
        threshold: 1%
    patch:
      default:
        target: auto
        threshold: 1%

Read the docs

Interact now

The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.

SVG

Learn more about the coverage sunburst graph

120	120		// objects for storing results
121	121		daily_values = vector<vector<vector<double>>>(params->n_demes, vector<vector<double>>(params->max_time));
122	122		// age distributions. Final level: 0 = Sh, 1 = Eh, 2 = Ah, 3 = Ch, 4 = Ph
123		-	age_distributions = vector<vector<vector<vector<double>>>>(params->n_output_age_times, vector<vector<vector<double>>>(params->n_demes, vector<vector<double>>(params->n_life_table, vector<double>(5))));
	123	+	age_distributions = vector<vector<vector<vector<double>>>>(params->n_output_age_times,
	124	+	vector<vector<vector<double>>>(params->n_demes,
	125	+	vector<vector<double>>(params->n_life_table, vector<double>(8))));
124	126
125	127		// misc
126	128		EIR = vector<double>(params->n_demes);
	129	+	prob_infectious_bite = vector<double>(params->n_demes);
	130	+	inc_infection = vector<double>(params->n_demes);
	131	+	inc_acute = vector<double>(params->n_demes);
	132	+	inc_chronic = vector<double>(params->n_demes);
127	133
128	134		}
129	135

206	213		// host-specific and changes over inoculations. Therefore, as a
207	214		// middleground, draw from Binomial(H[k],
208	215		// max_prob_infection*prob_infectious_bite), where max_prob_infection is
209		-	// the largest value that prob_infection could possibly take. Loop
	216	+	// the largest value that prob_infection could possibly take. Then loop
210	217		// through these query infectious bites and draw from a Bernoulli
211	218		// distribution relative to this value.
212	219		//
213	220		// For example, if prob_infection = {0.1, 0.05} then filter based on the
214	221		// value 0.1, i.e. draw the number of query infected hosts from
215		-	// Binomial(H[k], 0.1). Then loop these query hosts and draw from the
216		-	// relative probability of infection, which is Bernoulli with
217		-	// probability {0.1, 0.05}/0.1 = {1.0, 0.5}.
	222	+	// Binomial(H[k], 0.1*prob_infectious_bite). Then loop through these query
	223	+	// hosts and draw from the relative probability of infection, which is
	224	+	// Bernoulli with probability {0.1, 0.05}/0.1 = {1.0, 0.5}.
218	225
219		-	int host_query_infection = rbinom1(H[k], params->max_prob_infection * prob_infectious_bite);
	226	+	int host_query_infection = rbinom1(H[k], params->max_prob_infection * prob_infectious_bite[k]);
220	227		for (int i = 0; i < host_query_infection; ++i) {
221	228
222	229		// choose host at random

344	351		// update counts of each host status in each deme
345	352		update_host_counts(t);
346	353
	354	+	// update incidence measures
	355	+	update_incidence(t);
	356	+
347	357		// store daily values
348	358		for (int k = 0; k < params->n_demes; ++k) {
349	359		daily_values[k][t] = {double(H[k]), double(Sh[k]), double(Eh[k]), double(Ah[k]), double(Ch[k]), double(Ph[k]),
350	360		double(Sv[k]), double(Ev[k]), double(Iv[k]),
351		-	EIR[k],
	361	+	EIR[k], inc_infection[k], inc_acute[k], inc_chronic[k],
352	362		Ah_detectable_microscopy[k], Ch_detectable_microscopy[k],
353	363		Ah_detectable_PCR[k], Ch_detectable_PCR[k],n_inoc[k]};
354	364		}

356	366		// store age distributions
357	367		if (params->output_age_distributions && (params->output_age_times[index_age_distributions] == t+1)) {
358	368
	369	+	// get age distribution
359	370		get_age_distribution(index_age_distributions);
360	371
361	372		// update index

371	382		// get sampling parameters
372	383		int this_deme = params->ss_deme[index_obtain_samples];
373	384		int this_n = params->ss_n[index_obtain_samples];
	385	+	Diagnosis this_diagnosis = params->ss_diagnosis[index_obtain_samples];
374	386
375	387		// obtain samples
376		-	get_sample_details(t, this_deme, this_n);
	388	+	get_sample_details(t, this_deme, this_n, this_diagnosis);
377	389
378	390		// increment index
379	391		index_obtain_samples++;

192	198		// get number of new infectious bites on humans
193	199		EIR[k] = params->a * Iv[k] / double(H[k]);
194	200
195		-	double prob_infectious_bite = 1 - exp(-EIR[k]); // probability of new infectious bite on host
	201	+	// probability of new infectious bite on host
	202	+	prob_infectious_bite[k] = 1 - exp(-EIR[k]);
196	203
197	204		// one method of drawing infections in humans would be to draw the total
198	205		// number of infectious bites from Binomial(H[k], prob_infectious_bite),

447	459		H[k] = Sh[k] + Eh[k] + Ah[k] + Ch[k] + Ph[k];
448	460		}
449	461
	462	+	}
	463	+
	464	+	//------------------------------------------------
	465	+	// update incidence measures
	466	+	void Dispatcher::update_incidence(int t) {
	467	+
	468	+	// reset values
	469	+	fill(inc_infection.begin(), inc_infection.end(), 0.0);
	470	+	fill(inc_acute.begin(), inc_acute.end(), 0.0);
	471	+	fill(inc_chronic.begin(), inc_chronic.end(), 0.0);
	472	+
	473	+	// loop through all hosts, update incidence in given deme
	474	+	for (unsigned int i = 0; i < host_pop.size(); ++i) {
	475	+	int this_deme = host_pop[i].deme;
	476	+
	477	+	// incidence of infection
	478	+	double inc1 = prob_infectious_bite[this_deme] * host_pop[i].get_prob_infection() / double(H[this_deme]);
	479	+	inc_infection[this_deme] += inc1;
	480	+
	481	+	// incidence of acute and chronic infection
	482	+	double inc2 = host_pop[i].get_prob_acute();
	483	+	inc_acute[this_deme] += inc1 * inc2;
	484	+	inc_chronic[this_deme] += inc1 * (1.0 - inc2);
	485	+	}
450	486
451	487		}
452	488
453	489		//------------------------------------------------
454	490		// draw sample from deme
455		-	void Dispatcher::get_sample_details(int t, int deme, int n) {
	491	+	void Dispatcher::get_sample_details(int t, int deme, int n, Diagnosis diag) {
	492	+
	493	+	// n cannot exceed human population size at this point in time
	494	+	if (n > H[deme]) {
	495	+	n = H[deme];
	496	+	}
456	497
457	498		// draw vector by sampling without replacement
458	499		vector<int> samp = sample4(n, 0, H[deme] - 1);

465	506		int this_host_ID = host_pop[this_index].host_ID;
466	507
467	508		// find if positive for malaria parasites
	509	+	bool test_positive = false;
468	510		Status_host this_host_status = host_pop[this_index].get_host_status();
469		-	bool test_positive = (this_host_status == Host_Ah \|\| this_host_status == Host_Ch);
	511	+	if (this_host_status == Host_Ah) {
	512	+
	513	+	// get positive prob for acute microscopy vs PCR
	514	+	double prob_positive = 0.0;
	515	+	if (diag == microscopy) {
	516	+	prob_positive = host_pop[this_index].get_detectability_microscopy_acute(t);
	517	+	} else if (diag == PCR) {
	518	+	prob_positive = host_pop[this_index].get_detectability_PCR_acute(t);
	519	+	}
	520	+	test_positive = rbernoulli1(prob_positive);
	521	+
	522	+	} else if (this_host_status == Host_Ch) {
	523	+
	524	+	// get positive prob for chronic microscopy vs PCR
	525	+	double prob_positive = 0.0;
	526	+	if (diag == microscopy) {
	527	+	prob_positive = host_pop[this_index].get_detectability_microscopy_chronic(t);
	528	+	} else if (diag == PCR) {
	529	+	prob_positive = host_pop[this_index].get_detectability_PCR_chronic(t);
	530	+	}
	531	+	test_positive = rbernoulli1(prob_positive);
	532	+
	533	+	}
470	534
471	535		// save basic details
472	536		vector<int> this_details = {t+1, deme+1, this_host_ID, test_positive};

495	559		// get current time
496	560		int t = params->output_age_times[t_index];
497	561
498		-	// loop through all hosts, update age distribution
	562	+	// loop through all hosts, update age distributions
499	563		for (unsigned int i = 0; i < host_pop.size(); ++i) {
500	564		int this_deme = host_pop[i].deme;
501	565		int this_age = host_pop[i].get_age(t);

519	583		default:
520	584		Rcpp::stop("invalid host status in get_age_distribution()");
521	585		}
522		-	}
	586	+
	587	+	// incidence of infection
	588	+	double inc1 = prob_infectious_bite[this_deme] * host_pop[i].get_prob_infection() / double(H[this_deme]);
	589	+	age_distributions[t_index][this_deme][this_age][5] += inc1;
	590	+
	591	+	// incidence of acute and chronic infection
	592	+	double inc2 = host_pop[i].get_prob_acute();
	593	+	age_distributions[t_index][this_deme][this_age][6] += inc1 * inc2;
	594	+	age_distributions[t_index][this_deme][this_age][7] += inc1 * (1.0 - inc2);
	595	+
	596	+	} // end loop through hosts
523	597
524	598		}

378	378		#'
379	379		#' @param project a SIMPLEGEN project, as produced by the
380	380		#' \code{simplegen_project()} function.
381		-	#' @param x a dataframe containing all of the following columns:
	381	+	#' @param df_sample a dataframe containing all of the following columns:
382	382		#' \itemize{
383	383		#' \item time: the time (in days) at which samples are taken.
384	384		#' \item deme: the deme from which samples are taken.

393	393		#'
394	394		#' @export
395	395
396		-	define_sampling_strategy <- function(project, x) {
	396	+	define_sampling_strategy <- function(project, df_sample) {
397	397
398	398		# check inputs
399	399		assert_custom_class(project, "simplegen_project")
400		-	assert_dataframe(x)
401		-	assert_eq(names(x), c("time", "deme", "case_detection", "diagnosis", "n"))
402		-	assert_pos_int(x$time, zero_allowed = FALSE)
403		-	assert_pos_int(x$deme, zero_allowed = FALSE)
404		-	assert_in(x$case_detection, c("active", "passive"))
405		-	assert_in(x$diagnosis, c("microscopy", "PCR"))
406		-	assert_pos_int(x$n, zero_allowed = FALSE)
	400	+	assert_dataframe(df_sample)
	401	+	assert_eq(names(df_sample), c("time", "deme", "case_detection", "diagnosis", "n"))
	402	+	assert_pos_int(df_sample$time, zero_allowed = FALSE)
	403	+	assert_increasing(df_sample$time)
	404	+	assert_pos_int(df_sample$deme, zero_allowed = FALSE)
	405	+	assert_in(df_sample$case_detection, c("active", "passive"))
	406	+	assert_in(df_sample$diagnosis, c("microscopy", "PCR"))
	407	+	assert_pos_int(df_sample$n, zero_allowed = FALSE)
407	408
408	409		# specify formats
409		-	x$case_detection <- as.character(x$case_detection)
410		-	x$diagnosis <- as.character(x$diagnosis)
	410	+	df_sample$case_detection <- as.character(df_sample$case_detection)
	411	+	df_sample$diagnosis <- as.character(df_sample$diagnosis)
411	412
412	413		# load into project
413		-	project$sampling_strategy <- x
	414	+	project$sampling_strategy <- df_sample
414	415
415	416		invisible(project)
416	417		}

468	469		assert_single_logical(overwrite_transmission_record)
469	470		assert_single_logical(output_daily_counts)
470	471		assert_single_logical(output_age_distributions)
471		-	assert_vector(output_age_times)
472		-	assert_pos_int(output_age_times)
	472	+	assert_vector_pos_int(output_age_times, zero_allowed = FALSE)
473	473		assert_leq(output_age_times, max_time)
474	474		assert_single_logical(pb_markdown)
475	475		assert_single_logical(silent)

486	486		stop("no epi parameters defined. See ?define_epi_params")
487	487		}
488	488
	489	+	# check that inputs are compatible with sampling strategy
	490	+	if (!is.null(project$sampling_strategy)) {
	491	+	ss <- project$sampling_strategy
	492	+
	493	+	assert_greq(max_time, max(ss$time), message = "%s exceeded by sampling time (%s)")
	494	+	assert_greq(nrow(project$epi_parameters$mig_mat), max(ss$deme), "number of demes defined in sampling strategy exceeds number defined in epi parameters")
	495	+	}
489	496
490	497		# ---------- define argument lists ----------
491	498

559	566		ret <- as.data.frame(cbind(seq_len(nrow(ret)), i, ret))
560	567		names(ret) <- c("time", "deme", "H", "S", "E", "A", "C", "P",
561	568		"Sv", "Ev", "Iv",
562		-	"EIR",
	569	+	"EIR", "inc_infection", "inc_acute", "inc_chronic",
563	570		"A_detectable_microscopy", "C_detectable_microscopy",
564	571		"A_detectable_PCR", "C_detectable_PCR","n_inoc")
565	572		return(ret)

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

572	579		age_distributions <- do.call(rbind, mapply(function(j) {
573	580		ret <- do.call(rbind, mapply(function(i) {
574	581		ret <- do.call(rbind, output_raw$age_distributions[[j]][[i]])
575		-	colnames(ret) <- c("S", "E", "A", "C", "P")
	582	+	colnames(ret) <- c("S", "E", "A", "C", "P", "inc_infection", "inc_acute", "inc_chronic")
576	583		data.frame(cbind(deme = i, age = seq_len(nrow(ret)) - 1, ret))
577	584		}, seq_along(output_raw$age_distributions[[j]]), SIMPLIFY = FALSE))
578	585		cbind(sample_time = j, ret)

260	260
261	261		# subset data
262	262		dat <- project$epi_output$age_distributions
263		-	dat <- dat[dat$sample_time == sample_time & dat$deme == deme, c("age", state)]
	263	+	dat <- dat[(dat$sample_time == sample_time) & (dat$deme == deme), c("age", state)]
264	264		names(dat) <- c("age", "y")
265	265
266	266		# produce plot

474	474		for (int i = 0; i < params->max_inoculations; ++i) {
475	475
476	476		// anything that is currently in the acute or chronic stages is cured
477		-	if (inoc_status_asexual[i] == Acute_asexual \|\| inoc_status_asexual[i] == Chronic_asexual) {
	477	+	if ((inoc_status_asexual[i] == Acute_asexual) \|\| (inoc_status_asexual[i] == Chronic_asexual)) {
478	478
479	479		// if due to become infective at a future timepoint then store this
480	480		// timepoint

509	509		// prophylactic period is cured immediately upon emergence
510	510		if (inoc_status_asexual[i] == Liverstage_asexual) {
511	511
512		-	if (inoc_events[i].count(Event_Eh_to_Ah) != 0 && inoc_events[i][Event_Eh_to_Ah] <= t2) {
	512	+	// acute
	513	+	if ((inoc_events[i].count(Event_Eh_to_Ah) != 0) && (inoc_events[i][Event_Eh_to_Ah] <= t2)) {
513	514
514	515		// store time at which due to emerge
515	516		int t_emerge = inoc_events[i][Event_Eh_to_Ah];

522	523		new_inoc_event(t_emerge, Event_Ah_to_Sh, i);
523	524		new_inoc_event(t_emerge, Event_end_infective, i);
524	525
525		-	} else if (inoc_events[i].count(Event_Eh_to_Ch) != 0 && inoc_events[i][Event_Eh_to_Ch] <= t2) {
	526	+	// chronic
	527	+	} else if ((inoc_events[i].count(Event_Eh_to_Ch) != 0) && (inoc_events[i][Event_Eh_to_Ch] <= t2)) {
526	528
527	529		// store time at which due to emerge
528	530		int t_emerge = inoc_events[i][Event_Eh_to_Ch];

544	546		t_next_inoc_event = params->max_time + 1;
545	547		for (int i = 0; i < params->max_inoculations; ++i) {
546	548		for (const auto & x : inoc_events[i]) {
	549	+
	550	+	// catch to ensure that events cannot predate current time
547	551		if (x.second < t) {
548	552		print("time =", t);
549	553		print_inoc_status();

814	818		//------------------------------------------------
815	819		// get current probability of infection
816	820		double Host::get_prob_infection() {
	821	+
	822	+	// situations in which zero chance of infection
	823	+	if ((get_n_active_inoc() == params->max_inoculations) \|\| prophylaxis_on) {
	824	+	return 0.0;
	825	+	}
	826	+
	827	+	// get probability from flexible distribution
817	828		int index = (infection_index < params->n_prob_infection) ? infection_index : params->n_prob_infection - 1;
818	829		return params->prob_infection[index];
819	830		}
820	831		imilarity index 97%
821	832		ename from tests/testthat/test-assertions_v8.R
822	833		ename to tests/testthat/test-assertions_v9.R

78	78
79	79		// misc
80	80		std::vector<double> EIR;
	81	+	std::vector<double> prob_infectious_bite;
	82	+	std::vector<double> inc_infection;
	83	+	std::vector<double> inc_acute;
	84	+	std::vector<double> inc_chronic;
81	85
82	86		// number of active inoculations
83	87		std::vector<double> n_inoc;

91	95		void init(Parameters &params);
92	96		void run_simulation(Rcpp::List &args_functions, Rcpp::List &args_progress);
93	97		void update_host_counts(int t);
94		-	void get_sample_details(int t, int deme, int n);
	98	+	void update_incidence(int t);
	99	+	void get_sample_details(int t, int deme, int n, Diagnosis diag);
95	100		void get_age_distribution(int t_index);
96	101
97	102		};

1		-
2		-	#------------------------------------------------
3		-	#' @title Define new SIMPLEGEN project
4		-	#'
5		-	#' @description Define a new SIMPLEGEN project. This project will hold all
6		-	#' simulation inputs and outputs for a given analysis, and is initialised with
7		-	#' the default values of all parameters.
8		-	#'
9		-	#' @export
10		-
11		-	simplegen_project <- function() {
12		-
13		-	# create empty project
14		-	project <- list(epi_parameters = NULL,
15		-	sampling_strategy = NULL,
16		-	epi_output = NULL,
17		-	sample_details = NULL,
18		-	relatedness = NULL,
19		-	true_genotypes = NULL,
20		-	observed_genotypes = NULL)
21		-
22		-	class(project) <- "simplegen_project"
23		-
24		-	# return
25		-	invisible(project)
26		-	}
27		-
28		-	#------------------------------------------------
29		-	# overload print() function for simplegen_project
30		-	#' @method print simplegen_project
31		-	#' @export
32		-	print.simplegen_project <- function(x, ...) {
33		-
34		-	# print summary
35		-	summary(x)
36		-
37		-	# return invisibly
38		-	invisible(x)
39		-	}
40		-
41		-	#------------------------------------------------
42		-	# overload summary() function for simplegen_project
43		-	#' @method summary simplegen_project
44		-	#' @export
45		-	summary.simplegen_project <- function(object, ...) {
46		-
47		-
48		-
49		-	message("TODO - some default print method")
50		-
51		-	# return invisibly
52		-	invisible(object)
53		-	}
54		-
	1	+
	2	+	#------------------------------------------------
	3	+	#' @title Define new SIMPLEGEN project
	4	+	#'
	5	+	#' @description Define a new SIMPLEGEN project. This project will hold all
	6	+	#' simulation inputs and outputs for a given analysis, and is initialised with
	7	+	#' the default values of all parameters.
	8	+	#'
	9	+	#' @export
	10	+
	11	+	simplegen_project <- function() {
	12	+
	13	+	# create empty project
	14	+	project <- list(epi_parameters = NULL,
	15	+	sampling_strategy = NULL,
	16	+	epi_output = NULL,
	17	+	sample_details = NULL,
	18	+	relatedness = NULL,
	19	+	true_genotypes = NULL,
	20	+	observed_genotypes = NULL)
	21	+
	22	+	class(project) <- "simplegen_project"
	23	+
	24	+	# return
	25	+	invisible(project)
	26	+	}
	27	+
	28	+	#------------------------------------------------
	29	+	# overload print() function for simplegen_project
	30	+	#' @method print simplegen_project
	31	+	#' @export
	32	+	print.simplegen_project <- function(x, ...) {
	33	+
	34	+	# print summary
	35	+	summary(x)
	36	+
	37	+	# return invisibly
	38	+	invisible(x)
	39	+	}
	40	+
	41	+	#------------------------------------------------
	42	+	# overload summary() function for simplegen_project
	43	+	#' @method summary simplegen_project
	44	+	#' @export
	45	+	summary.simplegen_project <- function(object, ...) {
	46	+	p <- object
	47	+
	48	+	# if empty project
	49	+	if (all(mapply(is.null, p))) {
	50	+	message("(empty project)")
	51	+	invisible(object)
	52	+	}
	53	+
	54	+	# print epi parameters
	55	+	if (!is.null(p$epi_parameters)) {
	56	+	message("Epidemiological model:")
	57	+	n_demes <- length(p$epi_parameters$H)
	58	+	message(sprintf(" demes: %s", n_demes))
	59	+	message(sprintf(" H:\t %s", paste(p$epi_parameters$H, collapse = ", ")))
	60	+	message(sprintf(" M:\t %s", paste(p$epi_parameters$M, collapse = ", ")))
	61	+	message(sprintf(" seed infections: %s", paste(p$epi_parameters$seed_infections, collapse = ", ")))
	62	+	}
	63	+
	64	+	# print sampling strategy
	65	+	if (!is.null(p$sampling_strategy)) {
	66	+	message("Sampling strategy:")
	67	+	n_time <- unique(p$sampling_strategy$time)
	68	+	n_samp_time <- mapply(sum, split(p$sampling_strategy$n, f = p$sampling_strategy$time))
	69	+	if (length(n_time) <= 5) {
	70	+	message(sprintf(" time: %s", paste(n_time, collapse = ", ")))
	71	+	message(sprintf(" n: %s", paste(n_samp_time, collapse = ", ")))
	72	+	} else {
	73	+	message(" time: (more than 5)")
	74	+	message(" n: (more than 5)")
	75	+	}
	76	+
	77	+	}
	78	+
	79	+	invisible(object)
	80	+	}
	81	+

Updating our web app

1	comment: false
2
3	coverage:
4	status:
5	project:
6	default:
7	target: auto
8	threshold: 1%
9	patch:
10	default:
11	target: auto
12	threshold: 1%