First picked up on here

A rate can mean a lot of things, but in epi. it's most often interpreted as the number of events per unit time. Instead, a probability refers to the likelihood of an event happening to an individual.

At the very least we should document Bitset$sample to have a prob argument instead of a rate argument. But we should change the implementation variable names to make that less confusing for development.

DoubleVariables can only store one number per individual. But sometimes you want several numbers to represent a variable. So that you can access and update them together. e.g.

# there are 4 aspects of diarrhoea we want to model
diarrhoea_variables <- list(
  DoubleVariable$new(rep(0, parameters$population)), # diarrhoea_bacteria_prior 
  DoubleVariable$new(rep(0, parameters$population)), # diarrhoea_virus_prior
  DoubleVariable$new(rep(0, parameters$population)), # diarrhoea_parasite_prior
  DoubleVariable$new(rep(0, parameters$population)) # diarrhoea_rotavirus_prior
)

# in process...
for (v in diarrhoea_variables) {
  prior <- v$get_values()
  # do something with this aspect of diarrhoea (or perhaps combine with another aspect)
  v$queue_update(new_prior, changed_individuals)
}

It would be much easier to write:

diarrhoea_priors <- DoubleVectorVariable$new(matrix(0, ncol=4, nrow=parameters$population))

# in process...
priors <- diarrhoea_priors$get_values(
  individual_index = NULL, # get all individuals
  column_index = NULL # get all columns
)
# process your priors
diarrhoea_priors$queue_update(
  values = new_priors # matrix of updated priors
  individual_index = NULL, # update all individuals
  column_index = NULL # update all columns
) # will update the priors matrix at the end of the timestep

Clients want to mock Variables so that they can put values and test updates. Like so:

my_variable <- mock_double(rep(0, 4))
mockery::expect_args(
    my_variable$queue_update,
    1,
    expected_value,
    expected_index
)

Insertions and iterations to hasmaps are two of the main bottlenecks at the moment. See the latest profile

Create state intermittently fails. It's very hard to debug/recreate.

A version without Rcpp::Function or Rcpp::Environment would be better.

@giovannic I have an implementation of the shadow integer variable here. The way it works is pretty straightforward, there are just 2 vectors it contains, each time step an enum is used to track which one is "active" and which is the "shadow". Queuing updates overwrites the shadow. The update sets the shadow to be the new active vector by swapping the enum and copies the old shadow to the new shadow vector.

The result of one of the benchmarks is here (the pattern is similar for the rest). As you can see, for large variables when only a few elements are being updated, the current implementation is faster. For larger variables when a high proportion of elements are being updated however, the shadow variable is faster, see limit (variable size) of 1e6 with updating 9e5 elements in the lower right, where the shadow variable is quite a bit faster.

vv_bi.pdf

So it's a tradeoff = ) Of course the double variable's shadow implementation should look almost identical. It may be less of a tradeoff for the double variable, which is probably updated much more frequently, with more of the population being updated.

Currently Bitset$remove accepts a vector of elements, not indices to remove. For certain purposes, it would be useful to have another method which accepts either a vector of indices or perhaps a Bitset of indices to flip to zero. This occurs in cases where you have a Bitset, say x which is used to get values y and z from DoubleVariable objects, and then you only want to keep things in x fulfilling some logical condition, like y > x (added a use case example in case there's some good obvious solution I missed!)

This issue is related to openjournals/joss-reviews#3539.

A statement of need is missing in the software documentation (readme).

Currently, the only way to get the number of individuals in a state (or number of states) is

length(api$get_state(human, severe_covid, very_severe_covid))
# will get the combined index for individuals in the severe_covid and very_severe_covid state
# and then in R it will get the length of that vector

But I would like to do that without requesting the whole state vector

api$get_state_size(human, severe_covid, very_severe_covid)
# will quickly return the number of individuals in the severe_covid and very_severe_covid states

If the user misspecifies an individual, forgetting for example an event, the resulting error when running simulate() (which calls processes including the missing event) is ambiguous. I think the function that errors has access to the event$name, so it would be helpful to the user for the error message to state that a named avent is missing (if possible?). The following reprex produces this error:

# Load packages
library(individual)

# Produce outputs
render_state_sizes <- function(api) {
  api$render('S', length(api$get_state(human, S)))
  api$render('I', length(api$get_state(human, S)))
}
# Infection process
infection <- function() {
  function(api) {
     # Isolate susceptible individuals
    from_state <- api$get_state(human, S)
    # Infection with given probability
    api$queue_state_update(human, I, from_state[runif(length(from_state)) < 0.5])
  }
}
# Recovery process
recovery <- function(api) {
  function(api){
    # Isolate susceptible individuals
    from_state <- api$get_state(human, I)
    # Infection with given probability
    api$queue_state_update(human, S, from_state[runif(length(from_state)) < 0.2])
  }
}
# Create death event for scheduled natural death
death_event <- Event$new('death')
# Add the action associated with death (int his case "birth" = new susceptible)
death_event$add_listener(function(api, target) {
  api$queue_state_update(human, S, target)
})
# Schedule deaths for any indivduals where death is not scheduled
death <- function(lifespan){
  function(api) {
    # All individuals
    alive <- api$get_state(human, S, I)
    # Those with death scheduled
    death_already_scheduled <- api$get_scheduled(death_event)
    # Those (newborn) individuals without death scheduled
    death_not_scheduled <- setdiff(alive, death_already_scheduled)
    # Draw random lifespan for newborns
    lifespan <- round(rexp(length(death_not_scheduled), 1 / lifespan))
    # Schedule death
    api$schedule(death_event, death_not_scheduled, lifespan)
  }
}
# Define processes
processes <- list(
  infection(),
  death(lifespan = 1000),
  render_state_sizes
)
# Set population and states
population <- 1000
S <- State$new('S', population - 10)
I <- State$new('I', 10)
human <- Individual$new('human', list(S, I), events = list(death_event))

# Works fine
t1 <- simulate(human, processes, 100)


# Define but with missing event
human2 <- Individual$new('human', list(S, I))
# Ambiguous error:
t2 <- simulate(human2, processes, 100)

The API adds a lot of boilerplate to adding new variables (like categorical, vector-based, multi-category). The API is also complicated to mock.

So lets refactor it out. The high level changes will include:

• No more SimAPI or api objects. Clients will instead make their variables and pass them into their processes.
• No more State object. This is just another variable, it's not special. There could be several state-like variables.

So the basic SIR example goes from:

population <- 4
S <- State$new('S', population)
I <- State$new('I', 0)
R <- State$new('R', 0)
human <- Individual$new('human', list(S, I, R))

transition <- function(from, to, rate) {
  return(function(api) {
    from_state <- api$get_state(human, from)
    api$queue_state_update(
      human,
      to,
      from_state[runif(length(from_state), 0, 1) < rate]
    )
  })
}

processes <- list(
  transition(S, I, .2),
  transition(I, R, .1),
  transition(R, S, .05),
  state_count_renderer_process(individual$name, c('S', 'I', 'R'))
)

simulate(human, processes, 5)

To...

population <- 4
timesteps <- 5
state <- CategoricalVariable$new(c('S', 'I', 'R'), rep('S', population))
renderer <- Render$new(timesteps)

transition <- function(from, to, rate) {
  return(function(api) {
    from_state <- state$get_index_of(from)
    state$queue_update(
      to,
      from_state$sample(rate)
    )
  })
}

processes <- list(
  transition('S', 'I', .2),
  transition('I', 'R', .1),
  transition('R', 'S', .05),
  categorical_count_renderer_process(renderer, state, c('S', 'I', 'R'))
)

simulation_loop(variables=list(state), processes=processes, timesteps=timesteps)
renderer$to_dataframe()

Going forward, this decouples the variables and allows us to improve them and create new ones in isolation

Still needs to be specced with a concrete, reduced test case.

But api$queue_state_update should error if you try to update an index that's out of range

Long term thoughts, opening here for discussion/keeping a record of wishlist items.

Consider a situation where contact is structured by a network between individuals (different from metapopulation, which can be handled fairly easily by just keeping track of which node each individual is in at any time), an epidemic on a network. To simulate infection, we'd want to find discordant pairs, if we had a bitset for S and a bitset for I individuals, plus some sort of data structure keeping track of edges, we'd want to (in increasing complexity):

1. Simple infection: find for each S the number of I edges (if infection probability is identical for each S-I edge)
2. Heterogeneous infection: find for each S the bitset of I persons (if infection probability depended on some additional attribute of those individuals)

For the simple infection situation, once we have the number of I's connecting to each S, we can just use multi_probability_bernoulli_process to sample infections, or sample and schedule a future infection time if infectious contacts don't happen at the points of a Poisson process.

For the heterogeneous infection type, one would need to get, for each S, a subset of a DoubleVariable or IntegerVariable corresponding to the I's connected to it, which would then be used to sample infection, as the simple infection situation.

Depending on the relative number of S's and I's it may be better to find those I's connected to each S, or those S's connected to each I.

For keeping track of edges a few options are:

1. a bitset for each individual storing their contacts, memory hog. To get discordant pairs just AND this with I (or S).
2. a bitset with buckets (sparse) for each individual storing their contacts. To get discordant pairs just AND this with I (or S).
3. a sparse matrix, one could check if S > I or S < I and then decide which rows to check to get discordant pairs; that would involve somehow intersecting that row with I (or S).
4. something more complicated?

Hi! I have a question about how to properly use the Bitset class, in particular the to_vector member function seems to be giving odd results for me. An example of unexpected behavior is here.

> new <- Bitset$new(4)
> new$insert(15:18)
> new$to_vector()
[1] 0 0 0 0

What is the proper way to initialize a bitset when the number of desired integers to insert is known? Thanks!

Currently, something like the following (where variable is a CategoricalVariable isn't possible:

variable$queue_update(value = c("a", "b", "c", "c", "a"), index = 1:5)

and returns the Error:

Error: Expecting a single string value: [type=character; extent=5].

Would it be possible to modify this functionality so that CategoricalVariable queue_update could support this?

Noting here to include in the PR #139

> c <- CategoricalVariable$new(categories = c("A","B"),initial_values = rep(c("A","B"),each=10))
> c$get_index_of("A")$to_vector()
 [1]  1  2  3  4  5  6  7  8  9 10
> c$get_index_of("B")$to_vector()
 [1] 11 12 13 14 15 16 17 18 19 20
> c$queue_update(value = "A",index = c(15,25,50))
> c$.update()
> c$get_index_of("B")$to_vector()
[1] 11 12 13 14 16 17 18 19 20
> c$get_index_of("A")$to_vector()
 [1]  1  2  3  4  5  6  7  8  9 10 15  0  0

Just need to change insert to insert_safe in categorical_variable.cpp

Models often include binary variables that are not mutually exclusive (like State objects). These are currently modelled with Variables using 1 and 0. Investigate and integrate performant alternatives.

Hi individual team! I had a question/ request about the output formats

Instead of generating a data frame that contains the counts of each compartment at each time step, I'm interested in creating a data frame that shows each individual's status at each time point. Is there a way to do this within the Render class?

@giovannic I don't know if this is contributing to the issues coming up in the CRAN checks but I am seeing this coming up from building on Rhub. Was on this platform: Windows Server 2008 R2 SP1, R-release, 32/64 bit

In file included from ../inst/include/common_types.h:15,
                 from bitset.cpp:10:
../inst/include/IterableBitset.h: In instantiation of 'void bitset_sample_internal(IterableBitset<A>&, double) [with A = long long unsigned int]':
bitset.cpp:99:42:   required from here
../inst/include/IterableBitset.h:485:22: warning: comparison of integer expressions of different signedness: 'unsigned int' and 'int' [-Wsign-compare]
       while(bitset_i != i) {
             ~~~~~~~~~^~~~
../inst/include/IterableBitset.h: In instantiation of 'void bitset_choose_internal(IterableBitset<A>&, size_t) [with A = long long unsigned int; size_t = unsigned int]':
bitset.cpp:164:33:   required from here
../inst/include/IterableBitset.h:455:20: warning: comparison of integer expressions of different signedness: 'unsigned int' and 'int' [-Wsign-compare]
     while(bitset_i != i) {
           ~~~~~~~~~^~~~

It's coming from bitset_sample_internal and bitset_choose_internal. My guess is that the object returned from Rcpp::sample in this line isn't the right type to be compared with bitset_i, and because the iterator variable i "inherits" (sorry for loose language, I have no idea how auto type deduction works) from the return type of Rcpp::sample, we see the comparison between signed and unsigned types. That's my hypothesis at least.

The full build log is here: 00install.txt

We have a prefab for fixed transition rates:

Rcpp::XPtr<process_t> fixed_probability_state_change_process(
    const std::string individual,
    const std::string from_state,
    const std::string to_state,
    double rate) { ...

Would be really helpful to have a very similar one but for different rates for each individual:

Rcpp::XPtr<process_t> multi_probability_state_change_process(
    const std::string individual,
    const std::string from_state,
    const std::string to_state,
    std::vector<double> rate) { ...

Tilde can be confusing (relation to destructors)

A common IBM transition I can see being really useful would be for transitions where individuals transition at a rate (fixed or different per individual) into one of two states determined by an additional probability (maybe templating fixed or different per individual would be great).

Rcpp::XPtr<process_t> fixed_probability_forked_state_change_process(
    const std::string individual,
    const std::string from_state,
    const std::string to_state_1,
    const std::string to_state_2,
    double rate,
    double fork_probability) { 

    // Draw which individual will transition based on rate

    // Allocate drawn individuals to either to_state_1 or to_state_2 based on for_probability

}

Involves making a proper C++ interface for parameters

Invoking std::function is a major bottleneck see latest profile

If a client wants to test a process function, they have to write mocks for the R or C++ api. It would save lots of time if these mocks could be provided by the package.

An example of an R API mock can be found here and C++ API can be found here.

It would also be helpful to provide a vignette describing how you would test a process function.

Create an optimised variable to represent individuals who could be in a combination of states.

So we can have something like:

vaccinated_with = MulitCategory$new(
  c('AstraZeneca', 'Pfizer', 'Moderna'),
  matrix(FALSE, ncol=3, nrow=population) # initialise with a boolean matrix of shape (category, population)
)

vaccination_process <- function(timestep) {
  ...
  # add a new vaccine without overwriting existing vaccines
  vaccinated_with$queue_update(vaccine_type, vaccinated_population)
  ...
}

waning_process <- function(timestep) {
  ...
  # remove vaccine after it's no longer effective
  vaccinated_with$queue_removal(vaccine_type, vaccinated_population)
  ...
}

immunity_process <- function(timestep) {
  ...
  # do something with patients who have been vaccinated with a combination
  az <- vaccinated_with$get_index_of('AstraZeneca')
  pf <- vaccinated_with$get_index_of('Pfizer')
  both <- az$and(pf) # bitset and operator
  ...
}

This kind of behaviour requires bitset. Which is currently WIP. Issue to track

TODO:

• A MultiCategory C++ class, you could draw inspiration from CategoricalVariable
• queue_update logic allows for a client to have a new category without removing their existing one
• add queue_removal logic so clients can explicitly remove a category
• update logic executes updates and removals preserving FIFO queue order
• expose the C++ class to R in src/variable.cpp
• expose an R6 class to the client, you could draw inspiration from CategoricalVariable

When following the instructions in the docs (https://mrc-ide.github.io/individual/articles/Performance.html#c-prefabs) the resulting package cannot compile because it gets a "[type from individual] not declared in this scope" error. For some reason the header "individual.h" is not being included in the auto-generated RcppExports.cpp, so it doesn't see the type declarations. After finding out how to do it we will need to update the docs.

Everyone should know by the next version

The get_scheduled method in Event calls event_get_scheduled which wants an XPtr to a TargetedEvent but gets an Event passed to it, causing a crash when it tries to access things that don't exist.

Currently it is possible to queue an update to a Catergorical variable of a category that was not defined when the variable was initialised. This is very bug prone (e.g. with a misspelling of a queued update). Most of the other CategoricalVariable functionality seems to check and return a helpful error message if an unspecified catergory is used.

Simple reprex:

library(individual)
my_pets <- CategoricalVariable$new(categories = c("dog", "cat"), initial_values = "dog")
# Queue an update with a un-specified category (easy to do with a typo)
target <- Bitset$new(1)$insert(1)
my_pets$queue_update("hat", target)
my_pets$.update()

# No cats
my_pets$get_size_of("cat")
# No dogs either ... Where have my pets gone!!
my_pets$get_size_of("dog")

I'm working on this in this branch (along with some other improvements for events) https://github.com/slwu89/individual/tree/feat/events, creating an issue to remind myself to write tests to cover these cases.

example:

>   event <- TargetedEvent$new(10)
>   target <-1:5
>   delay <- NaN
> event$schedule(target = target,delay = delay)
target_timestep: 9223372036854775809---

I know this is documented as a feature, but it's highly prone to error.

The expected case is great:

# I want to reset everyone's immunity to 0:
api$queue_variable_update(human, immunity, 0)
# Everyone is set to 0

But it leads to unexpected behaviour if I give an index:

# I want to reset a random subset of human immunities to 0
n_died <- rbinom(10, .5)
died <- sample.int(10, n_died, replace=FALSE)
api$queue_variable_update(human, immunity, 0, died)
# If n_died == 0 everyone will be reset to 0!
# I really want no one to be set to 0

I propose a change to the R and C++ api to make variable updates with empty indices do nothing.

@giovannic for something in safir I realized it would be nice to add a method for Bitset objects that would let you randomly sample a N number of indices to retain. It can be done somewhat awkwardly right now like this:

n <- 100
bset <- Bitset$new(n)$insert(1:n)
num_to_keep <- 30
to_remove <- sample.int(n = n,size = n - num_to_keep,replace = FALSE)
bset$remove(to_remove)
bset$size() == num_to_keep

What do you think about making some function that acts on the IterableBitset where you tell it N, the number of items to retain in the bitset, and then it will randomly sample the indices of the ones to keep? (and accompanying R interface)? I imagine it could be virtually identical to bitset_sample_internal except the 2nd argument to Rcpp::sample would be b.size() - N.

An important question is what to call this method..."random_subset", "keep", or something else? Also, I just realized this could be accomplished with filter_bitset, but the performance benefit of a modifying method might be worth it.

It's hard to spec out and prioritise issues #29 #30. So run some profiles to spec them out.

You don't always want a target population when you are scheduling something for the future

I think it might be a good idea to store the set of valid states within CategoricalVariable. I have a use case where I'll need multiple categorical variables and having each object know its set of valid states makes rendering output easier. I'm considering writing a rendering prefab to take in an indeterminate number of categorical variables to output counts of each, which would be greatly eased by storing states within the variables (for outputting the marginals or intersections).

We currently don't have a good way to render output from double variable objects. For example, it would be nice if someone could do something like the following:

library(individual)
tsteps <- 50
x <- DoubleVariable$new(initial_values = rlnorm(n = 100))

process <- function(timestep) {
  oldx <- x$get_values()
  newx <- oldx + oldx*0.01
  x$queue_update(values = newx)
}

renderer <- Render$new(timesteps = tsteps)

render_process <- function(timestep) {
  renderer$render(name = "my_double_variable",value = x$get_values(),timestep = timestep)
}

simulation_loop(
  variables = list(x=x),
  processes = list(process, render_process),
  timesteps = tsteps
)

The design of Render$.vectors is for variables with discrete finite states. @giovannic have you worked out any good way to render output for DoubleVariable objects in the malaria sim package? I remember that has quite a few float variables, such as the various immunity types.

@giovannic I'm posting this here as a reminder for myself to fix this when I have time; noticed it in safir debugging.

We need to add a check to make sure that the max_size of the bitset is the same as the IntegerVariable (presumably this problem also occurs for DoubleVariable). Otherwise you can get a returned vector of nonsense values because it reads memory past the end of the vector.

library(individual)
x <- IntegerVariable$new(initial_values = 1:100)
b <- Bitset$new(1000)$insert(90:110)
x$get_values(b)

Clients can currently schedule an event for several individuals at a single time in the future.

Clients want to schedule a slightly different time for each individual in the future.

Their choice of model. (likely exponential)

Current version of individual has this behavior:

> variable <- DoubleVariable$new(1:10)
> variable$get_values()
 [1]  1  2  3  4  5  6  7  8  9 10
> variable$queue_update(values = 1:50)
> variable$.update()
> variable$get_values()
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
[30] 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

This does not agree with how it should work (this should be an error, according to the docs/design principles). I'm working on this and some enhancements for variable update speed ("shadow variable" in my fork https://github.com/slwu89/individual/tree/feat/variable-enhancements) and will correct this and include tests to check for this bug as part of that work.

• Investigate get_state and get_variable usages
• Figure out ways to perform c++ operations on them and reduce copies

This is related to openjournals/joss-reviews#3539

The plain text file for the LICENSE should contain this text:
https://opensource.org/licenses/MIT

In api$get_state we use unlist(...) which means that when debugging and passing in state names, things go badly. But it's not totally clear why unlist is needed here, as this level of list index is probably always wanted? Can it be removed?

Common processes include:

• fixed probability state transition
• exponential decay