[v10] refactor: use ImplicitDiscreteSystem to implement affects in callback systems
#3452
Conversation
vyudu
force-pushed
the
callback
branch
2 times, most recently
from
4a6f2f6 to
b7997b0
Compare
vyudu
changed the title
ImplicitDiscreteSystem to implement affects in callback systems ImplicitDiscreteSystem to implement affects in callback systems
|
@ChrisRackauckas @AayushSabharwal This is very close, just needs a release for ImplicitDiscreteSolve and getting all the doc blocks working. One of the FMI tests is failing because the reinitialization isn't accurate enough (like 1.008e-8 but the tolerance is 1e-8), but I am unsure if that is related. |
vyudu
changed the title
ImplicitDiscreteSystem to implement affects in callback systems ImplicitDiscreteSystem to implement affects in callback systems
|
@AayushSabharwal I think this is basically ready to merge into the v10 branch. The FMI failure is the thing we discussed where the initialization fails with some error very close to the tolerance; the other test failures seem unrelated |
|
|
||
| @mtkbuild sys = ODESystem( | ||
| D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]]) | ||
| D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ Pre(c) + 1]]) |
There was a problem hiding this comment.
Is c updated in this case, even if it is not saved? That's actually pretty useful if it works this way.
| end | ||
| outvar = :p | ||
| function (u, t, integ) | ||
| if DiffEqBase.isinplace(integ.sol.prob) |
There was a problem hiding this comment.
| if DiffEqBase.isinplace(integ.sol.prob) | |
| if DiffEqBase.isinplace(SciMLBase.get_sol(integ).prob) |
This seems safer
| else | ||
| initialize = handle_optional_setup_fn( | ||
| map(fn -> fn.initialize, affect_functions), SciMLBase.INITIALIZE_DEFAULT) | ||
| return CallbackSet(compiled_callbacks...) |
There was a problem hiding this comment.
I see. Yeah, this is fine then. Thanks for the clarification.
| u_getters = [getsym(sys, aff_map[u]) for u in dvs_to_access] | ||
| p_getters = [getsym(sys, unPre(p)) for p in ps_to_access] | ||
| u_setters = [setsym(sys, u) for u in dvs_to_update] | ||
| p_setters = [setsym(sys, p) for p in ps_to_update] | ||
| affu_getters = [getsym(affsys, sys_map[u]) for u in dvs_to_update] | ||
| affp_getters = [getsym(affsys, sys_map[p]) for p in ps_to_update] |
There was a problem hiding this comment.
All of these can just call getsym or setsym with an array of symbolics. That way the getter just returns an array of values, instead of having to call an array of getters.
| pmap = Pair[p => getp(integ) for (p, getp) in zip(ps_to_access, p_getters)] | ||
| u0map = Pair[u => getu(integ) | ||
| for (u, getu) in zip(dvs_to_access, u_getters)] | ||
| affprob = remake(affprob, u0 = u0map, p = pmap, tspan = (integ.t, integ.t)) |
There was a problem hiding this comment.
Symbolic remake is still very slow. We know exactly what the u0 and p should be for this, so they should be passed directly to remake. This can leverage setsym_oop. So we'd have ugetter = getsym(sys, ...) and usetter = setsym_oop(affsys, ...) and do
new_us = ugetter(integ)
new_u0, _ = usetter(affprob, new_us)And similarly,
new_ps = pgetter(integ)
_, new_p = psetter(affprob, new_ps)Then affprob = remake(affprob, u0 = new_u0, p = new_p, tspan = (integ.t, integ.t))
The reason I suggest two different usetter and psetter is because getsym(sys, [dvs..., ps...]) is not type-stable, so we need two getters and consequently two setters.
|
|
||
| alg_eqs = filter(eq -> eq.lhs isa Union{Symbolic, Number} && !is_diff_equation(eq), | ||
| deqs) | ||
| cont_callbacks = to_cb_vector(continuous_events; CB_TYPE = SymbolicContinuousCallback, | ||
| iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) | ||
| disc_callbacks = to_cb_vector(discrete_events; CB_TYPE = SymbolicDiscreteCallback, | ||
| iv = iv, alg_eqs = alg_eqs, warn_no_algebraic = false) |
There was a problem hiding this comment.
Since we're doing the same thing in the SDEProblem constructor (and should do it in DDEProblem and SDDEProblem as well) can this code block be extracted into a function that is called in all 4 problem constructors? That way we don't have to be concerned about them falling out of sync.
Close #2639 #2612