Recycling trajectory issues when using 'now' function

[chrisreading01]

Any help/advice much appreciated!

I have a simmer model for a patient waiting list, and depending how long a patient has been waiting they are assigned to a different capacity pot. E.g. a patient waiting for treatment at 44 days is added to the 40-80 capacity queue, a patient at 85 days is added to the 80-120 day capacity queue, and so on. This is all done using the now() function to compare the current simulation time with the patient's arrival time. It's to help quantify the impact of not always treating patients strictly from longest wait to shortest.

The trajectory the patient follows is dynamically created using some trajectory blocks stored as strings, and it takes 50-60 seconds to parse together. I would like to only do this once, and then recycle that trajectory across however many iterations I simulate. But re-using the trajectory poses difficulties when trying to functionalise the process, I presume because of the different environments in play.

See code examples below, which are a simplified version just printing the current simulation time in different ways:

## This way works as expected - creating a trajectory and a simmer in global env
t1 <- trajectory() %>% 
  timeout(1) %>% 
  log_(function(){as.character(simmer::now(sim))})

sim <- simmer() 
sim %>% 
  add_generator("test",t1,at(0,1,2,3,4)) %>% 
  run()


## This way also works - creating a trajectory and a simmer within function environment
runASim1 <- function(){
  t2 <- trajectory() %>% 
    timeout(1) %>% 
    log_(function(){as.character(simmer::now(sim2))})
  
  sim2 <- simmer() 
  sim2 %>% 
    add_generator("test",t2,at(0,1,2,3,4)) %>% 
    run()
}

runASim1()


## This way falls over -- error says that sim3 cannot be found when executing simmer::now()
# ideally this is the way which would let me create the trajectory once, and then recycle as many times as needed
t3 <- trajectory() %>% 
  timeout(1) %>% 
  log_(function(){as.character(simmer::now(sim3))})

runASim2 <- function(traj){
  sim3 <- simmer() 
  sim3 %>% 
    add_generator("test",traj,at(0,1,2,3,4)) %>% 
    run()
} 
  
runASim2(t3)

In the last example, even when the function environment contains an object called 'sim3', the trajectory doesn't think it exists when it tries to run now(sim3). The middle example works fine, but in my actual model this would mean I'd have to recreate the trajectory inside every function call, at a cost of 50-60 seconds a time.

Any suggestions on how I could share the trajectory across these environments, or if there's a better way to achieve this would be gratefully received :)

[Enchufa2]

Yes, functions in trajectories are tied to the environment where the trajectory was defined, and there's also some cleaning involved to avoid reference loops that could make your memory grow indefinitely.

A simple solution to your third example would be something like this:

sim3 <- NULL # <---
t3 <- trajectory() %>% 
  timeout(1) %>% 
  log_(function(){as.character(simmer::now(sim3))})

runASim2 <- function(traj){
  sim3 <<- simmer() # <---
  sim3 %>% 
    add_generator("test",traj,at(0,1,2,3,4)) %>% 
    run()
} 

runASim2(t3)

Note that the first line is not strictly needed, but it helps if you define the trajectory in another place different from the global environment. The best solution also depends on whether you want different simulations, or multiple replications of the same simulation (in this case, I would just define a single simulation and then call reset and run repeteadly), and whether you plan to parallelize (in this case, you can only share the trajectories on Unix systems with multicore).

BTW, I'm curious why you have the trajectories as strings. What's the application?

[chrisreading01]

Thanks so much for the rapid response! The global assignment I'm guessing will only work if the iterations are run sequentially? My goal is to parallelize the simulations, my current cruder version of the model runs 50 iterations in mclapply but has to reparse the trajectories each time. It's running on an ubuntu server. Is there a way to achieve this in parallel, without re-building the trajectories each time?

The trajectories as strings was a crude solution to dynamically building the trajectories. The application is a patient waiting list simulation, and as we don't have document patient treatment pathways which could be 'hard coded' as trajectories, I have a set of SQL/R scripts which extract real patient journeys (e.g. Clock Start -> Outpatient New appt -> Inpatient admission -> Clock Stop) and then convert these into text strings of predefined simmer trajectory bricks. The text bricks are stored in a SQL table, and then retrieved and parsed/evaluated when needed for use in a simulation.

For each of our patient waiting lists there is a surprising diversity of patient journeys, so there would have been too many trajectories to set up manually, so this was how I ended up 'mining' our historic data to generate trajectories dynamically. I've attached an example diagram of how messy the variety of trajectories can be!

[Enchufa2]

The approach above should work for mclapply without issue, with the only addition that you should add wrap() after run() (see Replication). It won't work for any other type of parallelization though.

As for the trajectories, when I have to build dynamic trajectories, I parametrize them in a function, and then call this function n times with a grid of parameters (see examples here). It should be much faster than parsing and evaluating text. And you could also try |> instead of %>%, but the latter should be pretty fast already these days.

But I don't think it's worth it to generate a distinct trajectory per patient (that's my understanding at least). I would implement the graph that you shared as your only trajectory for this problem: I just see a bunch of visit blocks in a chain of 7 branches. You could even implement them with a single branch and a rollback, keeping track of the number of jumps with an attribute. Something like this:

probs <- array()

t <- trajectory() %>%
  set_attribute("current", 1) %>% # 1 is "clock start"
  set_attribute("iteration", 1) %>%
  branch(
    tag = "branch",
    function() {
      idx <- get_attribute(sim, c("current", "iteration"))
      sample(1:n, 1, prob=probs[idx[1], , idx[2]])
    }, continue = rep(TRUE, n),
    trajectory() %>%
      set_attribute("current", 1) %>%
      visit("op 1", function() ...), # sample from proper distribution
    trajectory() %>%
      set_attribute("current", 2) %>%
      visit("op 2", function() ...), # sample from proper distribution
    trajectory() %>%
      set_attribute("current", 3) %>%
      visit("op 3", function() ...), # sample from proper distribution
    ...,
    trajectory() %>% # n is "stop"
      leave(1)
  ) %>%
  set_attribute("iteration", 1, mod="+") %>%
  rollback("branch")

where probs should be a multidimensional array, i.e. layers of matrices that represent the probability of transitioning from i to j at iteration k. You have the numbers in your diagram, so you just need to divide them by the total number of patients.

[chrisreading01]

Thank you, I'll give the wrap/run approach a go (I've seen the wrap function in the documentation but struggled a bit to understand what it's for).

I'll definitely have a look at your example for dynamic trajectories. The text/parsing approach was a working solution given my limited R knowledge when I set up the first version of this model. It's not producing a trajectory per patient, but it does have a trajectory for every unique start to finish journey through the above diagram, a lot of which is repetitive. I previously used a single trajectory with LOTS of branching within the model, but the volume of random number generation in each branching point was causing slowdown in big simulations.

That style of branching trajectory for N activities is a really elegant solution which I would never have thought of. Once my current model version is working, I'll come back to this and hopefully it can increase the efficiency and bypass my parsing time problem altogether 👍 Getting quite excited!

And thank you for maintaining this package and providing support, it's already transforming how my NHS trust does waiting list forecasting :) It was great to see you on the NHS-R Community conference last year.

[Enchufa2]

Thanks, nice to see a growing community of users! :)

About wrap(), the technical explanation is that external pointers (i.e. objects created in the C++ side) cannot be exported/imported between R sessions). This method collects all the results from a simulation (so that they can be extracted later with get_mon_*()) getting rid of the simulator itself, and therefore this can be exported (and shared between workers when you are parallelizing).