[WIP] forward sensitivity transform #398
Conversation
|
Do we need this transform? It feels like we can achieve the same thing by providing |
|
Personally, in my work I need to do a lot of nested differentiation towards different subsets of parameters, up to 4th order sometimes. And going from |
I mean, that is ForwardSensitivity, but is that optimal? What are the equations like? How do you go higher order? This will be pretty much as fast as it gets, so it'll be a good thing to have for multiple reasons, benchmarking and all of that. It's also a good teaching tool since it generates the equations in a tangible way. |
| iv = sys.iv() | ||
| f = du_dt = [eq.rhs for eq ∈ sys.eqs] | ||
| u = [state(iv) for state ∈ states(sys)] | ||
| sen_symbols = [Symbol(:d,state.op.name,:_d,par.op.name) for state ∈ u, par in p] |
There was a problem hiding this comment.
The canonical naming we're using for this is https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/first_order_transform.jl#L3 which is what would match Differential(t)(x) in the lowered form
There was a problem hiding this comment.
Ideally what I would want is to use the transform to generate a system like this:
using ModelingToolkit
using DifferentialEquations
@parameters t
@parameters u₀
@parameters k
@variables u(t)
@derivatives Dₜ'~t
@derivatives Dᵤ₀'~u₀
a = Dₜ(u)
b = Dᵤ₀(u)
eqs = [Dₜ(u) ~ k*u,
Dₜ(Dᵤ₀(u)) ~ 1.0]
sys = ODESystem(eqs)
u0 = [u => 2.0,
u₀=> 1.0]
p = [k => 1.0]
tspan = (0.0,1.0)
prob = ODEProblem(sys,u0,tspan,p)
sol = solve(prob,Tsit5())
using Plots; plot(sol)I think this would would make it easier to detect identical states (witch switched order of partial derivatives) in higher order derivatives.
There was a problem hiding this comment.
Yes, that's what using the convention allows to work. Breaking the convention would break this example.
| u = [state(iv) for state ∈ states(sys)] | ||
| sen_symbols = [Symbol(:d,state.op.name,:_d,par.op.name) for state ∈ u, par in p] | ||
| du_dp = [Variable(sen_symbol)(iv) for sen_symbol ∈ sen_symbols] | ||
| @derivatives D'~iv |
There was a problem hiding this comment.
| @derivatives D'~iv | |
| D = Differential(iv.name) |
so the naming is kept.
|
One advantage of this approach compared to A disadvantage is that I can not think of a straightforward way to do |
|
I could use some advice on this: is it possible/useful to use this transform combined with |
No because adaptive time stepping and all of that wouldn't be handled well. That's better handled by AD tooling and such. |
| u = [state(iv) for state ∈ states(sys)] | ||
| sen_symbols = [Symbol(:d,state.op.name,:_d,par.op.name) for state ∈ u, par in p] | ||
| du_dp = [Variable(sen_symbol)(iv) for sen_symbol ∈ sen_symbols] | ||
| @derivatives D'~iv |
There was a problem hiding this comment.
| @derivatives D'~iv | |
| D = Differential(iv.op.name) |
| sen_symbols = [Symbol(:d,state.op.name,:_d,par.op.name) for state ∈ u, par in p] | ||
| du_dp = [Variable(sen_symbol)(iv) for sen_symbol ∈ sen_symbols] |
There was a problem hiding this comment.
| sen_symbols = [Symbol(:d,state.op.name,:_d,par.op.name) for state ∈ u, par in p] | |
| du_dp = [Variable(sen_symbol)(iv) for sen_symbol ∈ sen_symbols] | |
| names = [Symbol(state.op.name,:_,par.op.name) for state ∈ u, par ∈ p] | |
| du_dp = [Variable(name)(iv) for name ∈ names] |
This does not specify the parametric type of Variable{Type}. Is basing this on vartype(state.op) good enough? Does the type of par not matter?
Is it really needed to use this naming scheme at this level already? Can we not keep the sensitivity states as Differential types?
D_p = [Differential(par.op.name) for par ∈ p]
[D_par(state) for state ∈ u, D_par ∈ D_p]
I think this would be beneficial for higher order derivatives.
There was a problem hiding this comment.
I think you can just keep them as differentials. That should work.
Should work for forward sensitivity towards (a subset) of the parameters present in the dynamic function.
Sensitivity towards initial conditions also works out of the box by adding them to
p(thendf/dp = 0).If a parameter influencing the dynamics is also thought to affect the initial condition of the system, this can be taken into account by appropriately modelling
u_0 = h(p).