From fa25f04f8b200ddf563379885cf4b94a77a0a741 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 5 Dec 2024 03:12:26 +0530
Subject: [PATCH] fix: fix `remake_initialization_data` on problems with no
 initprob

---
 src/systems/nonlinear/initializesystem.jl | 26 +++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl
index 726d171bd0..6e51a198c6 100644
--- a/src/systems/nonlinear/initializesystem.jl
+++ b/src/systems/nonlinear/initializesystem.jl
@@ -303,6 +303,18 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p,
                 diff_idxs[i] || continue
                 u0map[dvs[i]] = newu0[i]
             end
+        else
+            # The user did pass `u0` to `remake`, so they want to use
+            # the specified values. If we fill the rest of `u0map` with
+            # old values, the initialization can end up overdetermined
+            # (e.g. if the user specified a value for an algebraic variable)
+            # and wants to solve for a differential variable
+
+            # instead, fill the guesses with old values. That way, the user's
+            # u0map are the only constraints (along with defaults) and
+            # if the system is underdetermined, variables will likely retain
+            # their old values.
+            merge!(guesses, Dict(dvs .=> newu0))
         end
         if p === missing
             # the user didn't pass `p` to `remake`, so they want to retain
@@ -311,6 +323,20 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p,
             for p in ps
                 pmap[p] = getp(sys, p)(newp)
             end
+        else
+            # we can't just fill guesses for parameters like we did for unknowns
+            # because that can unintentionally make parameters solvable. If the
+            # parameter is solvable with the existing specification, fill its
+            # guess and fill `pmap` otherwise.
+            pmap_fallbacks = Dict()
+            for p in ps
+                if is_parameter_solvable(p, pmap, defs, guesses)
+                    guesses[p] = getp(sys, p)(newp)
+                else
+                    pmap_fallbacks[p] = getp(sys, p)(newp)
+                end
+            end
+            add_fallbacks!(pmap, ps, pmap_fallbacks)
         end
         # all non-solvable parameters need values regardless
         for p in ps