Parameters to function in argument events to trajectory function

[munozav]

Dear Aaron, Thanks so much for your great package. At this moment I am passing an event function, that is dependable on the params list, into the argument events to trajectory.

The events function works perfectly if before running the trajectory function, params list is saved in R environment with a name different to "params", let's say "pardef" and it is called by this name in the event function.

However currently, I am trying to run several simulations that only differ on the params list, so I want to create an R function that globally runs all the related pomp functions in one step and whose input is the params list "pardef".

In previous issues, you explained that "additional arguments to those functions are simply passed down to deSolve::ode when the deterministic skeleton is a vector field" but unfortunately if I run the global function, params or pardef is not recognized inside the event function that I am passing.

It would appreciate whatever information that can help me on this issue.

[kingaa]

@munozav: I'm afraid pomp doesn't really support the full range of dynamical systems that deSolve does. In particular, pomp only allows for deterministic skeletons that are either vectorfields or maps. It is possible, using deSolve, to handle a broader class of systems, including those for which state variables undergo finite jumps at specified "events". Why not use deSolve directly, eliminating the pomp middle-man?

If you really must use pomp, perhaps you can do something cheeky, like define a global pardef matrix as a copy of the parameter matrix you wish to use, then have your events function read from this. You could pass the column number of the pardef matrix as a state variable (with zero derivative). It's ugly, but it should work.

[munozav]

Many thanks for your prompt reply, unfortunately, I require to use pomp, because of its great functionality :).
Let me see if I understood you well. The parameter set that I want to pass into the event function is the same as the one introduced into the argument params in the trajectory function.

pardef.v<-unlist(pardef)

trajectory( pomp.det,
                   params = pardef,
                    events = list(func = eventfunc,
                                           time =time.window)

eventfunc<- function(times, y, params) {
               N<- pardef.v[“Beta”]
               x <- y * pardef.v[“Epsilon”] 
               ….. 
               }

Do you know where I should introduce pardef.v in the trajectory function? Many thanks for any advise.

[kingaa]

I was imagining you would define a matrix of parameters pardef, such that each column is a distinct parameter vector. Then, you would expand the state vector to include a new state variable, say index, and the parameter vector to contain a new parameter, say Index. You would define the derivative of index to be identically zero when specifying the deterministic skeleton. In rinit, index would simply be set equal to Index. Then, eventfunc could look something like

eventfunc <- function(times, y, params) {
    N <- pardef[“Beta”,as.integer(y["index"])]
    x <- y * pardef[“Epsilon”,as.integer(y["index"])] 
      …
}

You'd have to be careful to make sure that as.integer(Index) was what you wanted: it can be dodgy switching back and forth between integers and floating point numbers.

As I say, it's not pretty, but it should get the job done.

But all this would be unnecessary, no?, if one could get the parameters one wanted into eventfunc. What is the value of params in the eventfunc that you use? Remind me: why does the event function not receive the parameters?

[kingaa]

@munozav: could you provide a very simple example? I'm now curious as to why deSolve should not be able to handle this situation.

[kingaa]

@munozav: Consider the following:

library(pomp)
library(tidyverse)
theme_set(theme_bw())

pardef <- parmat(c(b=0.6,dose=10,Index=1),nrep=20)
pardef["b",] <- seq(0.1,1,length=20)
pardef["Index",] <- 1:20

tibble(
  time=seq(0,20,by=0.1)
) %>%
  mutate(x=NA) %>%
  pomp(
    times="time",
    t0=0,
    skeleton=vectorfield(
      Csnippet("Dconc = -b*conc; Dindex = 0;")
    ),
    rinit=Csnippet("conc = 20; index = nearbyint(Index);"),
    statenames=c("conc","index"),
    paramnames=c("b","Index")
  ) -> po

po %>%
  trajectory(
    format="d",
    params=pardef
  ) %>%
  left_join(
    pardef %>% melt() %>% spread(variable,value),
    by=c(index="Index")
  ) %>%
  ggplot(aes(x=time,y=conc,group=index,color=index))+
  geom_line()

po %>%
  trajectory(
    format="d",
    params=pardef,
    events=list(
      func=function(t,y,parms,...){
        yy <- matrix(y,nrow=2)
        i <- as.integer(yy[2,])
        yy[1,] <- yy[1,] + pardef["dose",i]
        dim(yy) <- NULL
        yy
      },
      times=seq(1,10,by=2)
    )
  ) %>%
  left_join(
    pardef %>% melt() %>% spread(variable,value),
    by=c(index="Index")
  ) %>%
  ggplot(aes(x=time,y=conc,group=b,color=b))+
  geom_line()

[munozav]

Aaron you are awesome!!, Thanks a lot for this example!!!. Many many thanks for your great package and for your kind support!!!!
Have a wonderful weekend.

[gokhale616]

Hi Dr. @kingaa ,

I have been faced with a problem similar to the one posted here. I am working with a trajectory matching problem and thus have rely on the trajectory() workhorse for my analysis.

Problem: I wish to generate a pulse of individuals entering into a state variable at a time t > t0.

Looking around, it became clear that pomp doesn't play nicely with all of the features available in the package desolve. One such feature is the events feature. Among other things, users can generate state-variable specific impulses at discrete time-points using this feature. As a work around, I decided to define my system as being externally forced with an impulse that moves a certain number of individuals to one of the state variables at a pre-specified time point.

As a minimal reproducible example, let's assume a compartmental model with 3 state variables (TIV). I initialize my system with no individuals in the V compartment. I wish to have 1 individual to magically migrate to V at time t = 2. For that, similar to the solution proposed above, I define a dummy variable V_event. This variable should remain 0 for all values of t except for the ones closest to t = 2. I say closest, because time specified by me, in this case 2, might not be used by the integrator to evaluate the solution.

For finicky a solution, this method should still work! Or generate at least one impulse. To verify whether I was able to find the right time, I have added an rprint() statement in my Csnippet that defines the skeleton. On execution, I could see that the dummy variables remain 0 for all values of time t except for the ones closest to t = 2.

But upon inspecting the trajectory output no changes are observed in the state variables! Am I missing something obvious here? Is there a better way of executing this solution? I am adding my codes for reference. I would greatly appreciate some help here. Thank you!

library(tidyverse)
library(ggplot)
library(pomp)

det_skel <- "
  /* ====== pseudo condition for introduction of the virus: events forcing ======*/
  double V_event;

  if(fabs(t-t_intro) < t_tol) {
    V_event = 1;
  } else {
    V_event = 0; 
  }

  Rprintf(\" V_event = %lg, t_val = %lg, abs_t_diff = %lg\\n\", V_event, t, fabs(t-t_intro));

  /*====== Balancing System of Ordinary Differential Equations ======*/
  /* Host cells */
  /*Compartmental shifts: Target Cells*/
  DT =  -beta_1*V*T;

  /*Compartmental shifts: Infected Cells*/
  DI = beta_1*V*T - delta*I; 

  /* Virus cells */
  /*Compartmental shifts: Virus particles*/
  DV = V_event + p*I - c*V;  
"

rinit <- "
  T = nearbyint(T_0);
  I = nearbyint(I_0);
  V = nearbyint(V_0);

" 

# State names:
state_names <- c("T", "I", "V")

# Parameter names:
rp_names <- c("beta_1", "delta", "p", "c", "t_intro", "t_tol")

ip_names <- c("T_0", "I_0", "V_0")

p_names <- c(rp_names, ip_names)

rp_vals <- c(beta_1 = 2e-5, delta = 3, p = 0.35, c = 20, t_intro = 2, t_tol = 1e-2)

ip_vals <- c(T_0 = 7e7, I_0 = 0, V_0 = 0)

p_vals <- c(rp_vals, ip_vals)

po <- data.frame(Days = seq(0, 7, by = 2), 
                 Obs = NA) %>% 
  pomp(t0 = 0, 
       times = "Days",
       skeleton = vectorfield(Csnippet(det_skel)),
       rinit = Csnippet(rinit), 
       statenames = state_names, 
       paramnames = p_names,
       params = p_vals, 
       #cdir = getwd(),
       #cfile = "min_TIV"

  ) 

po %>% 
  trajectory(format="d", include.data=FALSE, verbose = TRUE) %>% 
  select(-`.id`) %>% 
  gather(key = "Compartment", value = "Value", -Days) %>% 
  ggplot(aes(x = Days, y = Value)) +
  geom_line()+
  facet_wrap(.~Compartment)+
  theme_bw(base_size = 13)

[kingaa]

@gokhale616: You are correct that pomp, though it uses deSolve, does not attempt to provide all the functionality of the latter. However, pomp does support models that are time-inhomogeneous (i.e., non-autonomous), so one can do the job. However, when one has a non-smooth time-dependence (as you do) and one uses the sophisticated integrators in deSolve (which in effect assume smooth time-dependence), then one can run into the sort of problem you have encountered. The essence of the matter is this: the sophisticated integrators use variable-length timesteps to achieve speed without sacrificing accuracy. Because your forcing function V_event is zero outside of a small interval, the solvers have no awareness that it is not zero everywhere. In particular, unless they happen to evaluate this function within the interval [t_intro, t_intro+tol), they will never encounter a situation where V_event != 0.

You can work around this in a number of ways.

1. Force the solver to take smaller steps. One can do this using the hmax parameter. If you furnish hmax=tol to the solver, and you use the default integrator (lsoda), then the latter will never take a step larger than hmax, whence it will always necessarily evaluate V_event in the interval where it is nonzero.
2. Use flow to solve the equation over the interval of your choice. In essence, you recast the model as a discrete-time model, where the map in question is evaluated using flow. In this case, you could implement the impulsive forcing outside of the call to flow.
3. Define a covariate which is the number of impulses that have occurred at each time. Then, define a state variable which tracks the number of impulses that have occurred. Apply the impulse whenever the two do not match, then update the state variable.

Looking at these three options, I think the most efficient would be the last one. It might introduce a slight non-determinism into the trajectories, since you would not have full control over when the impulse is applied, however.

[gokhale616]

This made a lot of sense @kingaa.

Both the 1st and 3rd method were quite straight forward to implement and as you mention both of them introduce a slight non-determinism into the system trajectories.

Thanks lot for a quick response!

[kingaa]

Am closing this now. Feel free to reopen if more discussion is warranted.