ciemss / pyciemss

Causal and probabilistic reasoning with continuous time dynamical systems
Other
17 stars 6 forks source link

Bug report: Optimize integration interface gives IndexError when given multiple interventions #392

Closed liunelson closed 9 months ago

liunelson commented 1 year ago

When I provide more than one intervention to the load_and_optimize_and_sample_petri_model(...) function, I get the following IndexError: Screenshot 2023-10-17 at 18 10 51

For reference, these are the inputs:

num_samples = 5
timepoints = np.arange(5.0, dtype = float)
start_time = timepoints[0] - 1e-5

# Interventions over the parameters
# interventions = [(2.1, "beta")]
# intervention_bounds = [[0.0], [3.0]]
interventions = [(2.5, "beta"), (3.2, "beta")] # <=========
intervention_bounds = [[0.0, 0.5], [3.0, 1.0]]

# Initial condition
start_state = {}
for var in model['semantics']['ode']['initials']:
    start_state[var['target']] = var['expression']
    for param in model['semantics']['ode']['parameters']:
        start_state[var['target']] = start_state[var['target']].replace(param['id'], str(param['value']))
    start_state[var['target']] = float(start_state[var['target']])

# Optimization objective function to be minimized
# e.g. L1 norm of intervention parameters
objective_function = lambda x: np.sum(np.abs(x))

# Quantity of interest that is targeted by the optimization
# tuple of the form (callable function, state variable name, function arguments
# "scenario2dec_nday_average" computes the average of the n last days
qoi = ("scenario2dec_nday_average", "Infected_sol", 2)

optimize_result = load_and_optimize_and_sample_petri_model(
    model_path, 
    num_samples, 
    timepoints = timepoints, 
    interventions = interventions,
    bounds = intervention_bounds,
    qoi = qoi, 
    risk_bound = 10.0,
    objfun = objective_function,
    initial_guess = 1.0, # initial guess of the parameter for the optimizer
    verbose = True,
    start_time = start_time,
    start_state = start_state,
    n_samples_ouu = int(100),
    method = "euler",
    alpha_qs = [0.25, 0.5, 0.75, 0.95]
)

The model is just BIOMD0000000955_askenet.json

anirban-chaudhuri commented 10 months ago

This should be fixed once the new optimize interface is merged into main.