ukhsa-collaboration / pygom

ODE modelling in Python
GNU General Public License v2.0
27 stars 20 forks source link

The passage of time is not a stochastic process: periodic functions may present some dificulty when using .simulate_jump() #60

Open m-d-grunnill opened 2 years ago

m-d-grunnill commented 2 years ago

I am trying to think of how you would stochastically simulate a model with a periodic term within pygom. The common_models module intialisis a determiistic SIS model with periodic transmision in the following way.

state = ['S', 'E', 'I', 'tau']
param_list = ['mu', 'alpha', 'gamma', 'beta_0', 'beta_1']
derived_param = [('beta_S', 'beta_0 * (1 + beta_1*cos(2*pi*tau))')]
ode = [
    Transition(origin='S', equation='mu - beta_S*S*I - mu*S',
               transition_type=TransitionType.ODE),
    Transition(origin='E', equation='beta_S*S*I - (mu + alpha)*E',
               transition_type=TransitionType.ODE),
    Transition(origin='I', equation='alpha*E - (mu + gamma)*I',
               transition_type=TransitionType.ODE),
    Transition(origin='tau', equation='1',
               transition_type=TransitionType.ODE)
    ]
# initialize the model
ode_obj = DeterministicOde(state, param_list,
                           derived_param=derived_param,
                           ode=ode)

From this example I would intialise a stochastic version of this model in the following way:

actual_states=['S', 'E', 'I', 'R']
state = actual_states + ['time']
param_list = ['mu', 'alpha', 'gamma', 'beta_0', 'beta_1']
derived_param = [('beta_S', 'beta_0 * (1 + beta_1*cos(2*pi*time))')]
transitions = [
    Transition(origin='S', destination='E', equation='beta_S*S*I',
               transition_type=TransitionType.T),
    Transition(origin='E', destination='I', equation='alpha*E',
               transition_type=TransitionType.T),
    Transition(origin='I', estination='R', equation='gamma*I',
               transition_type=TransitionType.T),
    ]
bd_list = [
   Transition(origin=states['S'][i], equation='mu*(S + E + I + R)', transition_type=TransitionType.B)]
   Transition(origin='time', equation='1', transition_type=TransitionType.B)
   ]

bd_list += [Transition(origin=item, equation='mu*'+item, transition_type=TransitionType.D) for item in actual_states]
    # initialize the model
    model = SimulateOde(state, param_list,
                               derived_param=derived_param,
                               transition=transition,
                               birth_death=bd_list)

The issue is that this would stochisticise time if the model.simulate_jump() method was called.

twomagpi commented 3 months ago

Yes need to seriously revisit how time dependent processes are done. It is too complex for deterministic models and (probably) impossible for stochastic ones. @j-c-gibson I think you had some ideas on this?

j-c-gibson commented 3 months ago

Hopefully there is a solution which doesn't require us to define the dummy state for time. In the worst case, we can at least do this behind the scenes so that the end user doesn't have to.

The issue of time evolving stochastically is tough. In fact, I don't think there is any way for PyGOM to currently mix stochastically and deterministically evolving variables.

Also, we might wish to make a distinction between states (which are reserved for populations) and variables. Otherwise, defining time as a state is a bit tricky conceptually.