Merge pull request #136 from SciML/myb/more

Add optimal_parameter_intervention_for_reach

Author: ChrisRackauckas

src/EasyModelAnalysis.jl

@@ -25,6 +25,7 @@ export get_sensitivity, create_sensitivity_plot, get_sensitivity_of_maximum
 export stop_at_threshold, get_threshold
 export model_forecast_score
 export optimal_threshold_intervention, prob_violating_threshold,
-       optimal_parameter_intervention_for_threshold, optimal_parameter_threshold
+       optimal_parameter_intervention_for_threshold, optimal_parameter_threshold,
+       optimal_parameter_intervention_for_reach
 
 end

src/intervention.jl

@@ -165,3 +165,107 @@ function optimal_parameter_intervention_for_threshold(prob, obs, threshold,
     ss = duration_constraint(optx, [], Val(true))
     Dict(ps .=> optx), ss, ret
 end
+
+"""
+    optimal_parameter_intervention_for_reach(prob, obs, reach, cost, ps,
+        lb, ub, intervention_tspan, duration; ineq_cons = nothing, maxtime=60)
+
+## Arguments
+
+  - `prob`: An ODEProblem.
+  - `obs`: The observation symbolic expression.
+  - `reach`: The reach for the observation, i.e., the constraint enforces that `obs` reaches `reach`.
+  - `cost`: the cost function for minimization, e.g. `α + 20 * β`.
+  - `ps`: the parameters that appear in the cost, e.g. `[α, β]`.
+  - `lb`: the lower bounds of the parameters e.g. `[-10, -5]`.
+  - `ub`: the uppwer bounds of the parameters e.g. `[5, 10]`.
+  - `intervention_tspan`: intervention time span, e.g. `(20.0, 30.0)`. Defaults to `prob.tspan`.
+  - `duration`: Duration for the evaluation of intervention. Defaults to `prob.tspan[2] - prob.tspan[1]`.
+
+
+## Keyword Arguments
+
+  - `maxtime`: Maximum optimzation time. Defaults to `60`.
+  - `ineq_cons`: a vector of symbolic expressions in terms of symbolic
+      parameters. The optimizer will enforce `ineq_cons .< 0`.
+
+# Returns
+
+  - `opt_p`: Optimal intervention parameters.
+  - `(s1, s2, s3)`: Pre-intervention, intervention, post-intervention solutions.
+  - `ret`: Return code from the optimization.
+"""
+function optimal_parameter_intervention_for_reach(prob, obs, reach,
+                                                  symbolic_cost, ps, lb, ub,
+                                                  intervention_tspan = prob.tspan,
+                                                  duration = abs(-(prob.tspan...));
+                                                  maxtime = 60, ineq_cons = nothing,
+                                                  kw...)
+    t0 = prob.tspan[1]
+    ti_start, ti_end = intervention_tspan
+    symbolic_cost = Symbolics.unwrap(symbolic_cost)
+    #ps = collect(ModelingToolkit.vars(symbolic_cost))
+    _cost = Symbolics.build_function(symbolic_cost, ps, expression = Val{false})
+    _cost(prob.p) # just throw when something is wrong during the setup.
+
+    cost = let _cost = _cost
+        (x, grad) -> _cost(x)
+    end
+
+    function duration_constraint(x::Vector, grad::Vector, ::Val{p} = Val(false)) where {p}
+        tf = t0 + duration
+        prob_preintervention = remake(prob, tspan = (t0, ti_start))
+        if p
+            sol_preintervention = solve(prob_preintervention; kw...)
+        else
+            sol_preintervention = stop_at_threshold(prob_preintervention, obs, reach; kw...)
+            reach_time = ti_start - sol_preintervention.t[end]
+            reach_time > 0 && return sol_preintervention.t[end] - tf
+        end
+
+        prob_intervention = remake(prob, u0 = sol_preintervention.u[end], p = ps .=> x,
+                                   tspan = (ti_start, ti_end))
+        if p
+            sol_intervention = solve(prob_intervention; kw...)
+        else
+            sol_intervention = stop_at_threshold(prob_intervention, obs, reach; kw...)
+            reach_time = ti_end - sol_intervention.t[end]
+            reach_time > 0 && return sol_intervention.t[end] - tf
+        end
+
+        prob_postintervention = remake(prob, u0 = sol_intervention.u[end],
+                                       tspan = (ti_end, t0 + duration))
+        if p
+            sol_postintervention = solve(prob_postintervention; kw...)
+            sol_preintervention, sol_intervention, sol_postintervention
+        else
+            sol_postintervention = stop_at_threshold(prob_postintervention, obs, reach;
+                                                     kw...)
+            reach_time = tf - sol_postintervention.t[end]
+            10.0
+        end
+    end
+
+    opt = Opt(:GN_ISRES, length(ps))
+    opt.lower_bounds = lb
+    opt.upper_bounds = ub
+    opt.xtol_rel = 1e-4
+
+    opt.min_objective = cost
+    init_x = @. (lb + ub) / 2
+    duration_constraint(init_x, [])
+    inequality_constraint!(opt, duration_constraint, 1e-16)
+    if ineq_cons !== nothing
+        for con in ineq_cons
+            _con = Symbolics.build_function(Symbolics.unwrap(con), ps,
+                                            expression = Val{false})
+            _con(init_x)
+            ineq_con = (x, _) -> _con(x)
+            inequality_constraint!(opt, ineq_con, 1e-16)
+        end
+    end
+    opt.maxtime = maxtime
+    (optf, optx, ret) = NLopt.optimize(opt, init_x)
+    ss = duration_constraint(optx, [], Val(true))
+    Dict(ps .=> optx), ss, ret
+end

test/threshold.jl

@@ -59,6 +59,13 @@ opt_ps, s2, ret = optimal_parameter_threshold(prob, x, 3,
 opt_ps, s2, ret = optimal_parameter_threshold(prob, x, 3,
                                               -p, [p],
                                               [-1.0], [1.0]);
+opt_ps, (s1, s2, s3), ret = optimal_parameter_intervention_for_reach(remake(prob,
+                                                                            tspan = (0,
+                                                                                     1.0)),
+                                                                     x, 300,
+                                                                     p, [p],
+                                                                     [0.0], [30.0]);
+@test 10 < opt_ps[p] < 11
 @variables t x(t) y(t)
 @parameters p1 p2
 D = Differential(t)