Merge pull request #3605 from AayushSabharwal/as/copy-guesses

feat: use guesses as temporary values for variables solved by initialization

Author: AayushSabharwal

src/systems/problem_utils.jl

@@ -938,8 +938,7 @@ function maybe_build_initialization_problem(
         is_parameter_solvable(p, pmap, defs, guesses) || continue
         get(op, p, missing) === missing || continue
         p = unwrap(p)
-        stype = symtype(p)
-        op[p] = get_temporary_value(p, floatT)
+        op[p] = getu(initializeprob, p)(initializeprob)
         if iscall(p) && operation(p) === getindex
             arrp = arguments(p)[1]
             op[arrp] = collect(arrp)
@@ -948,7 +947,7 @@ function maybe_build_initialization_problem(
 
     if is_time_dependent(sys)
         for v in missing_unknowns
-            op[v] = get_temporary_value(v, floatT)
+            op[v] = getu(initializeprob, v)(initializeprob)
         end
         empty!(missing_unknowns)
     end

test/initial_values.jl

@@ -281,3 +281,31 @@ end
     @test prob.p isa Vector{Float64}
     @test length(prob.p) == 5
 end
+
+@testset "Temporary values for solved variables are guesses" begin
+    @parameters σ ρ β=missing [guess = 8 / 3]
+    @variables x(t) y(t) z(t) w(t) w2(t)
+
+    eqs = [D(D(x)) ~ σ * (y - x),
+        D(y) ~ x * (ρ - z) - y,
+        D(z) ~ x * y - β * z,
+        w ~ x + y + z + 2 * β,
+        0 ~ x^2 + y^2 - w2^2
+    ]
+
+    @mtkbuild sys = ODESystem(eqs, t)
+
+    u0 = [D(x) => 2.0,
+        x => 1.0,
+        y => 0.0,
+        z => 0.0]
+
+    p = [σ => 28.0,
+        ρ => 10.0]
+
+    tspan = (0.0, 100.0)
+    prob = ODEProblem(sys, u0, tspan, p, jac = true, guesses = [w2 => -1.0],
+        warn_initialize_determined = false)
+    @test prob[w2] ≈ -1.0
+    @test prob.ps[β] ≈ 8 / 3
+end

test/initializationsystem.jl

@@ -1186,10 +1186,10 @@ end
     @mtkbuild sys = ODESystem([D(x) ~ x * p + q, x^3 + y^3 ~ 3], t)
     prob = ODEProblem(
         sys, [], (0.0, 1.0), [p => 1.0]; guesses = [x => 1.0, y => 1.0, q => 1.0])
-    @test prob[x] == 0.0
-    @test prob[y] == 0.0
+    @test prob[x] == 1.0
+    @test prob[y] == 2.0
     @test prob.ps[p] == 1.0
-    @test prob.ps[q] == 0.0
+    @test prob.ps[q] == 3.0
     integ = init(prob)
     @test integ[x] ≈ 1 / cbrt(3)
     @test integ[y] ≈ 2 / cbrt(3)