Dynamic intervention behaviour should account for the slope at trigger threshold crossings

[liunelson]

Relevant to @anirban-chaudhuri @sabinala

Dynamic interventions are triggered by a model variable crossing some threshold value. In pyciemss, this involves a chirho ZeroEvent.

Pascale and Dan on our team raised this point recently:

• chirho ZeroEvent doesn't care whether the approach/crossing happens above or below the threshold https://github.com/ciemss/chirho/blob/master/chirho/dynamical/handlers/interruption.py
• If a model variable is non-monotonic, then a dynamic intervention would be triggered multiple times
• This also means I cannot actually use dynamic parameter interventions to represent the hospitalization-dependent mitigation policy in Monthly Demo Scenario (https://github.com/DARPA-ASKEM/program-milestones/tree/main/6-month-milestone/hackathon)
• More precisely, the "start-mitigation" dynamic parameter intervention would be triggered twice, when H (1) exceeds 0.8 * 3000 and (2) decreases below 0.8 * 3000; however, if I use a "end-mitigation" dynamic parameter intervention to end the mitigation policy, then it'd also get triggered twice and at the same time; both interventions would clash.
• I think we should be able to set the trigger variable, the trigger threshold, and a trigger slope, akin to how triggering works in oscilloscopes
  

[anirban-chaudhuri]

@liunelson @djinnome Is this issue resolved now?

[liunelson]

I don't think this issue is resolved.

@djinnome shared this chirho doc page about modeling state-dependent interventions
https://basisresearch.github.io/chirho/dynamical_intro.html#Modeling-a-State-Dependent-Intervention

@sabinala mentioned the interface notebooks https://github.com/ciemss/pyciemss/blob/main/docs/source/optimize_interface.ipynb

In both cases, I don't see how I could set up an intervention policy that runs while a target variable is greater than a threshold value (e.g. set beta to 0.05 while I > 100). To do so, I'd need be able to specify two interventions, (1) beta = 0.05 triggered when I increases from below to above 100 and (2) beta = 0.5 when I decreases from above to below 100.

The chirho example relies on R(t) to trigger an end to the policy, which is not how the problem is formulated it.

From the notebooks, it seems that, currently, I can only define dynamic interventions that:

• are triggered once, whenever var - thresh = 0 for the first time
• cannot be triggered multiple times even when var oscillates around thresh

Taking the example from interfaces.ipynb, how can I set up an intervention policy where (1) beta = 0.20 when I(t) increases above 150 and (2) beta = 0.25 when I(t) decreases below 180?

This screenshot shows the wrong behaviour:

[liunelson]

The intention is to create an intervention policy like this:

[sabinala]

@liunelson let me try coding up an example of this, it should be doable!

[sabinala]

@liunelson I've tried messing around with this quite a bit and I'm not finding a way around this issue. Let's take it up with Sam when he returns next week.

[SamWitty]

@liunelson , thanks for the detailed discussion and requests. This scenario is on the border of what ChiRho (and thus PyCIEMSS) can handle, but I think we can get a workable solution without major upstream changes.

chirho ZeroEvent doesn't care whether the approach/crossing happens above or below the threshold https://github.com/ciemss/chirho/blob/master/chirho/dynamical/handlers/interruption.py

This is true. In our examples thusfar we've been able to work around this by knowing something about the domain. For example, in general we know that infections increase over time, so if we have an event function defined by infections crossing some threshold then we know that that is approaching that threshold from below. If this ends up being a real blocker for this scenario then we could refine this.

If a model variable is non-monotonic, then a dynamic intervention would be triggered multiple times

This is not true. ChiRho's interruptions inherit from ShallowMessengers behind the scenes (see https://github.com/BasisResearch/chirho/blob/master/chirho/dynamical/internals/_utils.py#L85), which remove themselves after a single application. In other words, events only happen once. My suspicion is that this will resolve much of (if not all of) your issues with this scenario.

I think we should be able to set the trigger variable, the trigger threshold, and a trigger slope, akin to how triggering works in oscilloscopes

You absolutely can do this. You just need to have the dynamical system have a state variable for the derivative of the state variable, as you would anyway if you wanted to encode a higher order differential equation into its first order form. Then you would just have the event function refer to that state variable representing the derivative of another state variable. This does require modification of the model, which could in theory be automated to include terms for derivatives.

This screenshot shows the wrong behaviour:

Can you share the code that generated this behavior? It's hard to say what's going wrong as is.

[sabinala]

@SamWitty How would this work exactly?

"'I think we should be able to set the trigger variable, the trigger threshold, and a trigger slope, akin to how triggering works in oscilloscopes'

You absolutely can do this. You just need to have the dynamical system have a state variable for the derivative of the state variable, as you would anyway if you wanted to encode a higher order differential equation into its first order form. Then you would just have the event function refer to that state variable representing the derivative of another state variable. This does require modification of the model, which could in theory be automated to include terms for derivatives."

Let's say you want an intervention to trigger when a certain state variable crosses the threshold from above (with negative slope). Would you write the intervention like this? (var2 = name of derivative of var1 in the model)

# Define the threshold for when the intervention should be applied
def make_var_threshold(var: str, var2: str, threshold: torch.Tensor):
    def var_threshold(time, state):
        if state[var2] <0:
              return state[var] - threshold
       else
              return (something that is not zero)
    return var_threshold
    
infection_threshold = make_var_threshold("I", "rate_of_change_I", torch.tensor(150.0))
dynamic_parameter_interventions1 = {infection_threshold: {"p_cbeta": torch.tensor(0.3)}}

result = pyciemss.sample(model3, end_time, logging_step_size, num_samples, start_time=start_time, 
                         dynamic_parameter_interventions=dynamic_parameter_interventions1, 
                         solver_method="dopri5")

[SamWitty]

@sabinala Possibly, but it's generally better to avoid conditional branching when you want things to be autodiffed (as we do here).

Here's another way of doing this (that is very untested):

def make_var_threshold(var: str, slope_var: str, threshold: torch.Tensor):
    def var_threshold(time, state):
        # This line might not work, but would be easy to fix
        return state[var] - threshold + (state[slope_var] > 0)
    return var_threshold

[sabinala]

@SamWitty I would like to try your suggestion, but I haven't been able to make a state variable for the slope (state[slope_var]) with the templates available for an AMR. Is there something I'm missing?

[SamWitty]

@sabinala , I'm realizing now that this is a bit harder than I thought. It's easy to add the integral of an existing state variable to the state, but not the derivative.

[sabinala]

True, I made a few attempts at getting it to work, but no luck. Leaving this for now.