remake failures

[TorkelE]

Just trying to enumerate cases where remake fails and where there are some dependency between variables/parameters. Not 100% sure everything is declared correctly should work (and if not, would be useful to know so I can adjust my workflows). However, I have tried various versions, and I do think you always run into some errors one way or another.

If someone tells me what should/shout not work I can rework, and then we can incorporate into the tests to ensure stuff works.

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

# Defines the model.
@parameters k1 k2 V0
@variables X1(t) X2(t) Y1(t) Y2(t) V(t) = V0 W(t)
@parameters v = V w = W Γ[1:2] = missing [guess = [1.0, 1.0]]
eqs = [
    D(X1) ~ k1 * (Γ[1] - X1) - k2 * X1,
    D(Y1) ~ k1 * (Γ[2] - Y1) - k2 * Y1,
    D(V) ~ 1 - V,
    D(W) ~ 1 - W,
]
observed = [X2 ~ Γ[1] - X1, Y2 ~ Γ[2] - Y1]
@mtkbuild osys = ODESystem(eqs, t, [X1, X2, Y1, Y2, V, W], [k1, k2, V0, v, w, Γ]; observed)

# Create the`ODEProblem`s.
u0 = [X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0]
ps = [k1 => 0.1, k2 => 0.2, V0 => 3.0]
oprob1 = ODEProblem(osys, u0, 1.0, ps)
oprob2 = remake(oprob1, u0 = [X1 => 10.0])
oprob3 = remake(oprob2, u0 = [X2 => 20.0])
oprob4 = remake(oprob1, u0 = [X2 => 20.0, Y1 => 30.0])
oprob5 = remake(oprob1, u0 = [X1 => 10.0, X2 => 20.0])
oprob6 = remake(oprob1, u0 = [Y2 => 40.0], p = [k1 => 0.4])
oprob7 = remake(oprob1, u0 = [X1 => 10.0, X2 => 20.0], p = [V0 => 50.0])
oprob8 = remake(oprob1, u0 = [W => 60.0])

# Creates a testing function.
function test_vals(prob, us_correct::Dict, ps_correct::Dict)
    integ = init(prob)
    for u in keys(us_correct)
        @test prob[u] == us_correct[u]
        @test integ[u] == us_correct[u]
    end
    for p in keys(ps_correct)
        @test prob.ps[p] == ps_correct[p]
        @test integ.ps[p] == ps_correct[p]
    end
end

# Tests oprob1 without using test function.
@test oprob1[X1] == 1.0
@test oprob1[X2] == 2.0
@test oprob1[Y1] == 3.0
@test oprob1[Y2] == 4.0
@test oprob1[V] == 5.0
@test oprob1[W] == 6.0
@test oprob1.ps[k1] == 0.1
@test oprob1.ps[k2] == 0.2
@test oprob1.ps[V0] == 3.0
@test oprob1.ps[v] == 5.0 
@test oprob1.ps[w] == 6.0
@test oprob1.ps[Γ[1]] == 3.0
@test oprob1.ps[Γ[2]] == 7.0
integ1 = init(oprob1)
@test integ1[X1] == 1.0
@test integ1[X2] == 2.0
@test integ1[Y1] == 3.0
@test integ1[Y2] == 4.0
@test integ1[V] == 5.0
@test integ1[W] == 6.0
@test integ1.ps[k1] == 0.1
@test integ1.ps[k2] == 0.2
@test integ1.ps[V0] == 3.0
@test integ1.ps[v] == 5.0 
@test integ1.ps[w] == 6.0
@test integ1.ps[Γ[1]] == 3.0
@test integ1.ps[Γ[2]] == 7.0

# Test remaining `ODEProblems`.
test_vals(oprob1, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 7.0))
test_vals(oprob2, 
    Dict(X1 => 10.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 12.0, Γ[2] => 7.0))
test_vals(oprob3, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))
test_vals(oprob4, 
    Dict(X1 => 1.0, X2 => 20.0, Y1 => 30.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 21.0, Γ[2] => 34.0))
test_vals(oprob5, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))
test_vals(oprob6, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 30.0, Y2 => 40.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.4, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 43.0))
test_vals(oprob7, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 40.0, V => 50.0, W => 6.0), 
    Dict(k1 => 0.4, k2 => 0.2, V0 => 50.0, v => 55.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))
test_vals(oprob8, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 60.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 50.0, v => 50.0, w => 60.0, Γ[1] => 30.0, Γ[2] => 7.0))

[AayushSabharwal]

@test oprob1[X1] == 1.0
@test oprob1[X2] == 2.0
@test oprob1[Y1] == 3.0
@test oprob1[Y2] == 4.0
@test oprob1[V] == 5.0
@test oprob1[W] == 6.0
# ...

This requires #3290 to be fixed. Trivial initialization isn't run in the constructor yet. As such, for the explanations below I've updated the test_vals function to:

function test_vals(prob, us_correct::Dict, ps_correct::Dict)
    integ = init(prob)
    for u in keys(us_correct)
        # @test prob[u] == us_correct[u]
        @test integ[u] == us_correct[u]
    end
    for p in keys(ps_correct)
        # @test prob.ps[p] == ps_correct[p]
        @test integ.ps[p] == ps_correct[p]
    end
end

---

test_vals(oprob1, 
    Dict(X1 => 10.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 12.0, Γ[2] => 7.0))

Isn't this test incorrect? oprob1 was created with the constraint X1 => 1.0, so why would X1 be 10.0? And Gamma[1] would consequently be 3.0 (X1 + X2). V0 was also given as V0 => 3.0, so testing it is 5.0 doesn't make sense?

julia> test_vals(oprob1,
           Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 7.0))

This is the correct test call.

---

test_vals(oprob2, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

X1 => 10.0 is correct here, since that's the condition provided to remake. But why would X2 be 20.0? It retains the initial condition X2 => 2.0 provided to oprob1. Gamma[1] would also be 12.0 and V0 should be 3.0.

julia> test_vals(oprob2,
           Dict(X1 => 10.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 12.0, Γ[2] => 7.0))

---

test_vals(oprob3, 
    Dict(X1 => 1.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 21.0, Γ[2] => 7.0))

This is also incorrect. oprob2 had the condition X1 => 10.0, so oprob3 will retain this condition as well. Corresponding changes to Gamma[1] and the above error with V0 also apply.

julia> test_vals(oprob3,
           Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

---

test_vals(oprob4, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

oprob4 was made from oprob1, so it retains X1 => 1.0. It was given Y1 => 30.0 so that should be the initial condition here? Corresponding changes to Gamma[1], Gamma[2] and V0 also apply.

julia> test_vals(oprob4,
           Dict(X1 => 1.0, X2 => 20.0, Y1 => 30.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 21.0, Γ[2] => 34.0))

---

test_vals(oprob5, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 30.0, Y2 => 4.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 34.0))

Should be

julia> test_vals(oprob5,
           Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

For the same reasons. Retain the conditions of the original problem, add the constraints provided to remake, update solved parameters and V0.

---

Correcting the remake call for oprob6 from

oprob6 = remake(oprob1, u0 = [Y2 => 40.0], ps = [k1 => 0.4])

to

julia> oprob6 = remake(oprob1, u0 = [Y2 => 40.0], p = [k1 => 0.4])

test_vals(oprob6, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 40.0, V => 5.0, W => 6.0), 
    Dict(k1 => 0.4, k2 => 0.2, V0 => 5.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 43.0))

This is a complicated case. remake retains symbolic defaults. Because remake wasn't given a value for v, it uses the symbolic default V. Because remake wasn't given a value for V it uses the symbolic default V0. Thus, v gets a value of 3.0 (that of V0). Now, initialization is trivial and hence it runs and updates V to 5.0 since that is a hard constraint. But v isn't a solvable parameter, so it retains its value.

julia> test_vals(oprob6,
           Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 40.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.4, k2 => 0.2, V0 => 3.0, v => 3.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 43.0))

---

Similarly fixing the remake call of oprob7

test_vals(oprob7, 
    Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 50.0, W => 6.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 50.0, v => 50.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

V won't be 50 here, since the default V => V0 was overridden during construction of oprob1 by V => 5.0. So that is a hard constraint during initialization and V retains its value of 5.0.

julia> test_vals(oprob7,
           Dict(X1 => 10.0, X2 => 20.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 6.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 50.0, v => 50.0, w => 6.0, Γ[1] => 30.0, Γ[2] => 7.0))

---

test_vals(oprob8, 
    Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 60.0), 
    Dict(k1 => 0.1, k2 => 0.2, V0 => 5.0, v => 5.0, w => 60.0, Γ[1] => 3.0, Γ[2] => 7.0))

Since remake wasn't given the keyword argument p, it retains the original parameter object without modifications. The default w => W isn't respected. To perform this update, the call should be oprob8 = remake(oprob1, u0 = [W => 60.0], p = Dict()) or w should be made a solvable parameter. The fixed test call is

julia> test_vals(oprob8,
           Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 5.0, W => 60.0),
           Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 5.0, w => 6.0, Γ[1] => 3.0, Γ[2] => 7.0))

---

I hope this clarifies the issue.

[TorkelE]

Thanks a lot! Yeah, there was definitely some stuff that was messed up in the code. I have updated now.

Seems we will have to wait for #3290 and #3395 for a starter. Hopefully we can get this all working eventually. Either way, it might make sense to incorporate something like this into the test suite as well?

[AayushSabharwal]

#3395 is resolved in #3390. I'm working on fixing other CI issues in that PR, one of which is particularly elusive. #3290 is actually really easy to do - we just need to call prob = remake(prob) at the end of the constructors. I'll do that in a separate PR, since #3390 has been held up for a while now.