Unexpected results returned by Optimize

[liunelson]

I'm trying to do some test runs of Optimize following the interface notebook.

Here, I'm applying an intervention on a basic SIR model by changing the value of p_cbeta at t = 15.0 days, optimizing this value to minimize the risk that the 3-day average of I(t) at t = 50.0 days is > 50.0.

This function call runs without error.

However, the results are unexpected:

• a plot of the intervened parameter (persistent_p_cbeta) over time shows no step change in its value at t = 15.0 days or any other time; I expected a step-like curve.
• a plot of the state variables (S, I, R) over time show that an intervention is applied at t = 15.0; however, the intervention is on a rate constant and thus should only manifest as a change in the slope of the curves; what the plot shows seems to involve a change in the value of the state variables.

MODELS_PATH = "https://raw.githubusercontent.com/DARPA-ASKEM/simulation-integration/main/data/models/"
model3 = os.path.join(MODELS_PATH, "SIR_stockflow.json")

# Specify time
start_time = 0.0
end_time = 50.0
logging_step_size = 1.0

intervention_time = torch.tensor(15.0)

intervened_params = "p_cbeta"
p_cbeta_current = model3_tm.parameters[intervened_params].value
initial_guess_interventions = p_cbeta_current
bounds_interventions = [
    [model3_tm.parameters["p_cbeta"].distribution.parameters["minimum"]], 
    [model3_tm.parameters["p_cbeta"].distribution.parameters["maximum"]]
]

# Define QoI
observed_params = ["I_state"]
qoi = lambda x: obs_nday_average_qoi(x, observed_params, 3)
risk_bound = 50.0

objfun = lambda x: np.abs(p_cbeta_current - x)
static_parameter_interventions = {intervention_time: intervened_params}

opt_result = pyciemss.optimize(
    model3, 
    end_time, 
    logging_step_size, 
    qoi, 
    risk_bound, 
    static_parameter_interventions, 
    objfun, 
    initial_guess_interventions = initial_guess_interventions, 
    bounds_interventions = bounds_interventions, 
    start_time = 0.0, 
    n_samples_ouu = int(1e2), 
    maxiter = 1, 
    maxfeval = 20, 
    solver_method = "euler"
)

[liunelson]

@anirban-chaudhuri Does this make sense to you?

[anirban-chaudhuri]

I am not sure. Can you send me the entire file with the plotting code? I can take a look to see if something is going wrong.

[liunelson]

You can take a look here: https://github.com/liunelson/pyciemss/blob/nl-test/notebook/integration_demo/nliu/review_interfaces.py

[anirban-chaudhuri]

I believe the parameter value is not getting saved properly in the results dictionary when there is an intervention with the sample interface, That is why you are seeing the strange result for the p_cbeta plots. I have created #499 to address the issue with the saved parameter value. I am still looking through the optimize interface to see if it is giving strange results. I think it might be doing what it is supposed to but there is no policy that satisfies the constraint in the provided bounds. Let me check this though. When I ran it, it shows that the constraints are not satisfied. I believe for this SIR problem it is might not be possible to get the (risk associated with) infections <50.

[anirban-chaudhuri]

I think that optimize interface is working. I believe the issue was that there was no optimal policy to be found in the given range for the parameters. I modified the bound and the initial guess (so that it starts at a feasible location) and added your code along with the plots in https://github.com/ciemss/pyciemss/blob/ac-review-opt-TA4interface/docs/source/interfaces.ipynb (I added your code bits to the end of the notebook) @liunelson Can you take a look to see if this resolves the issue? We might have to give more guidance on how to interpret the optimize interface results (e.g., at least the OptResults key needs to be checked to see whether it gave a feasible solution).

[liunelson]

Your explanation from yesterday makes sense; switching to dopri5 and looking at the message in OptResults help me understand what happened.

[liunelson]

I'll get Terarium to surface the message in OptResult