diff --git a/Project.toml b/Project.toml
index 40ab38e07..617c99e80 100644
--- a/Project.toml
+++ b/Project.toml
@@ -60,6 +60,7 @@ DataFrames = "1.6"
 Distributed = "1.10"
 DocStringExtensions = "0.9"
 EnumX = "1"
+ForwardDiff = "0.10.36"
 FunctionWrappersWrappers = "0.1.3"
 IteratorInterfaceExtensions = "^1"
 LinearAlgebra = "1.10"
@@ -93,6 +94,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b"
@@ -109,4 +111,4 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Pkg", "Plots", "UnicodePlots", "PyCall", "PythonCall", "SafeTestsets", "Test", "StaticArrays", "StochasticDiffEq", "Aqua", "Zygote", "PartialFunctions", "DataFrames", "ModelingToolkit", "OrdinaryDiffEq"]
+test = ["Pkg", "Plots", "UnicodePlots", "PyCall", "PythonCall", "SafeTestsets", "Test", "StaticArrays", "StochasticDiffEq", "Aqua", "Zygote", "PartialFunctions", "DataFrames", "ModelingToolkit", "OrdinaryDiffEq", "ForwardDiff"]
diff --git a/test/downstream/modelingtoolkit_remake.jl b/test/downstream/modelingtoolkit_remake.jl
index fdc883f71..57cb0ae79 100644
--- a/test/downstream/modelingtoolkit_remake.jl
+++ b/test/downstream/modelingtoolkit_remake.jl
@@ -1,18 +1,22 @@
 using ModelingToolkit, SymbolicIndexingInterface
 using JumpProcesses
 using ModelingToolkit: t_nounits as t, D_nounits as D
+using Optimization
+using OptimizationOptimJL
 
-@parameters σ ρ β
+probs = []
+syss = []
+
+@parameters σ ρ β q
 @variables x(t) y(t) z(t)
 
-eqs = [D(D(x)) ~ σ * (y - x),
+eqs = [D(x) ~ σ * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z]
 
-@named sys = ODESystem(eqs, t)
-sys = structural_simplify(sys)
-u0 = [D(x) => 2.0,
-    x => 1.0,
+@named sys = ODESystem(eqs, t; parameter_dependencies = [q => 3β])
+sys = complete(sys)
+u0 = [x => 1.0,
     y => 0.0,
     z => 0.0]
 
@@ -21,145 +25,138 @@ p = [σ => 28.0,
     β => 8 / 3]
 
 tspan = (0.0, 100.0)
-oprob = ODEProblem(sys, u0, tspan, p, jac = true)
-
-@inferred typeof(oprob) remake(oprob; u0 = [x => 2.0], p = [σ => 29.0])
-oprob2 = remake(
-    oprob;
-    u0 = [x => 2.0, sys.y => 1.2, :z => 1.0],
-    p = [σ => 29.0, sys.ρ => 11.0, :β => 3.0]
-)
-@test oprob2.u0 isa Vector{<:Number}
-@test oprob2.p isa ModelingToolkit.MTKParameters
-@test oprob2[x] == oprob2[sys.x] == oprob2[:x] == 2.0
-@test oprob2[y] == oprob2[sys.y] == oprob2[:y] == 1.2
-@test oprob2[z] == oprob2[sys.z] == oprob2[:z] == 1.0
-@test getp(sys, σ)(oprob2) == 29.0
-@test getp(sys, sys.ρ)(oprob2) == 11.0
-@test getp(sys, :β)(oprob2) == 3.0
-
-oprob3 = remake(oprob; u0 = [x => 3.0], p = [σ => 30.0]) # partial update
-@test oprob3[x] == 3.0
-@test getp(sys, σ)(oprob3) == 30.0
-
-# SDEProblem.
-noiseeqs = [0.1 * x,
-    0.1 * y,
-    0.1 * z]
-@named noise_sys = SDESystem(sys, noiseeqs)
-noise_sys = complete(noise_sys)
-sprob = SDEProblem(noise_sys, u0, (0.0, 100.0), p)
-
-@inferred typeof(sprob) remake(sprob; u0 = [x => 2.0], p = [σ => 29.0])
-sprob2 = remake(
-    sprob;
-    u0 = [x => 2.0, sys.y => 1.2, :z => 1.0],
-    p = [σ => 29.0, sys.ρ => 11.0, :β => 3.0]
-)
-@test sprob2.u0 isa Vector{<:Number}
-@test sprob2.p isa ModelingToolkit.MTKParameters
-@test sprob2[x] == sprob2[sys.x] == sprob2[:x] == 2.0
-@test sprob2[y] == sprob2[sys.y] == sprob2[:y] == 1.2
-@test sprob2[z] == sprob2[sys.z] == sprob2[:z] == 1.0
-@test getp(sys, σ)(sprob2) == 29.0
-@test getp(sys, sys.ρ)(sprob2) == 11.0
-@test getp(sys, :β)(sprob2) == 3.0
-
-sprob3 = remake(sprob; u0 = [x => 3.0], p = [σ => 30.0]) # partial update
-@test sprob3[x] == 3.0
-@test getp(sys, σ)(sprob3) == 30.0
-
-# DiscreteProblem
-# @named de = DiscreteSystem(
-#     [D(x) ~ σ*(y-x),
-#     D(y) ~ x*(ρ-z)-y,
-#     D(z) ~ x*y - β*z],
-#     t,
-#     [x, y, z],
-#     [σ, ρ, β],
-# )
-# dprob = DiscreteProblem(de, u0, tspan, p)
-
-# @inferred typeof(dprob) remake(dprob; u0 = [x => 2.0], p = [σ => 29.0])
-# dprob2 = remake(
-#     dprob;
-#     u0 = [x => 2.0, sys.y => 1.2, :z => 1.0],
-#     p = [σ => 29.0, sys.ρ => 11.0, :β => 3.0]
-# )
-# @test dprob2.u0 isa Vector{<:Number}
-# @test dprob2.p isa ModelingToolkit.MTKParameters
-# @test dprob2[x] == dprob2[sys.x] == dprob2[:x] == 2.0
-# @test dprob2[y] == dprob2[sys.y] == dprob2[:y] == 1.2
-# @test dprob2[z] == dprob2[sys.z] == dprob2[:z] == 1.0
-# @test getp(de, σ)(dprob2) == 29.0
-# @test getp(de, sys.ρ)(dprob2) == 11.0
-# @test getp(de, :β)(dprob2) == 3.0
-
-# dprob3 = remake(dprob; u0 = [x => 3.0], p = [σ => 30.0]) # partial update
-# @test dprob3[x] == 3.0
-# @test getp(de, σ)(dprob3) == 30.0
-
-# NonlinearProblem
-@named ns = NonlinearSystem(
-    [0 ~ σ * (y - x),
-        0 ~ x * (ρ - z) - y,
-        0 ~ x * y - β * z],
-    [x, y, z],
-    [σ, ρ, β]
-)
-ns = complete(ns)
-nlprob = NonlinearProblem(ns, u0, p)
-
-@inferred typeof(nlprob) remake(nlprob; u0 = [x => 2.0], p = [σ => 29.0])
-nlprob2 = remake(
-    nlprob;
-    u0 = [x => 2.0, sys.y => 1.2, :z => 1.0],
-    p = [σ => 29.0, sys.ρ => 11.0, :β => 3.0]
-)
-@test nlprob2.u0 isa Vector{<:Number}
-@test nlprob2.p isa ModelingToolkit.MTKParameters
-@test nlprob2[x] == nlprob2[sys.x] == nlprob2[:x] == 2.0
-@test nlprob2[y] == nlprob2[sys.y] == nlprob2[:y] == 1.2
-@test nlprob2[z] == nlprob2[sys.z] == nlprob2[:z] == 1.0
-@test getp(ns, σ)(nlprob2) == 29.0
-@test getp(ns, sys.ρ)(nlprob2) == 11.0
-@test getp(ns, :β)(nlprob2) == 3.0
-
-nlprob3 = remake(nlprob; u0 = [x => 3.0], p = [σ => 30.0]) # partial update
-@test nlprob3[x] == 3.0
-@test getp(ns, σ)(nlprob3) == 30.0
-
-@parameters β γ
-@variables S(t) I(t) R(t)
-rate₁ = β * S * I
-affect₁ = [S ~ S - 1, I ~ I + 1]
-rate₂ = γ * I
-affect₂ = [I ~ I - 1, R ~ R + 1]
+
+push!(syss, sys)
+push!(probs, ODEProblem(sys, u0, tspan, p, jac = true))
+
+noise_eqs = [0.1x, 0.1y, 0.1z]
+@named sdesys = SDESystem(sys, noise_eqs)
+sdesys = complete(sdesys)
+
+push!(syss, sdesys)
+push!(probs, SDEProblem(sdesys, u0, tspan, p, jac = true))
+
+@named nsys = NonlinearSystem([0 ~ eq.rhs for eq in eqs], [x, y, z], [σ, β, ρ])
+nsys = complete(nsys)
+
+push!(syss, nsys)
+push!(probs, NonlinearProblem(nsys, u0, p, jac = true))
+
+rate₁ = β * x * y
+affect₁ = [x ~ x - σ, y ~ y + σ]
+rate₂ = ρ * y
+affect₂ = [y ~ y - 1, z ~ z + 1]
 j₁ = ConstantRateJump(rate₁, affect₁)
 j₂ = ConstantRateJump(rate₂, affect₂)
-j₃ = MassActionJump(2 * β + γ, [R => 1], [S => 1, R => -1])
-@named js = JumpSystem([j₁, j₂, j₃], t, [S, I, R], [β, γ])
+j₃ = MassActionJump(2 * β + ρ, [z => 1], [x => 1, z => -1])
+@named js = JumpSystem([j₁, j₂, j₃], t, [x, y, z], [σ, β, ρ])
 js = complete(js)
-u₀map = [S => 999, I => 1, R => 0.0]
-parammap = [β => 0.1 / 1000, γ => 0.01]
-tspan = (0.0, 250.0)
-jump_dprob = DiscreteProblem(js, u₀map, tspan, parammap)
-jprob = JumpProblem(js, jump_dprob, Direct())
-
-@inferred typeof(jprob) remake(jprob; u0 = [S => 900], p = [β => 0.2e-3])
-jprob2 = remake(
-    jprob;
-    u0 = [S => 900, js.I => 2, :R => 0.1],
-    p = [β => 0.2 / 1000, js.γ => 11.0]
-)
-@test jprob2.prob.u0 isa Vector{<:Number}
-@test jprob2.prob.p isa ModelingToolkit.MTKParameters
-@test jprob2[S] == jprob2[js.S] == jprob2[:S] == 900.0
-@test jprob2[I] == jprob2[js.I] == jprob2[:I] == 2.0
-@test jprob2[R] == jprob2[js.R] == jprob2[:R] == 0.1
-@test getp(js, β)(jprob2) == 0.2 / 1000
-@test getp(js, js.γ)(jprob2) == 11.0
-
-jprob3 = remake(jprob; u0 = [S => 901], p = [:β => 0.3 / 1000]) # partial update
-@test jprob3[S] == 901
-@test getp(js, β)(jprob3) == 0.3 / 1000
+jump_dprob = DiscreteProblem(js, u0, tspan, p)
+
+push!(syss, js)
+push!(probs, JumpProblem(js, jump_dprob, Direct()))
+
+@named optsys = OptimizationSystem(sum(eq.lhs for eq in eqs), [x, y, z], [σ, ρ, β])
+optsys = complete(optsys)
+push!(syss, optsys)
+push!(probs, OptimizationProblem(optsys, u0, p))
+
+k = ShiftIndex(t)
+@mtkbuild discsys = DiscreteSystem(
+    [x ~ x(k - 1) * ρ + y(k - 2), y ~ y(k - 1) * σ - z(k - 2), z ~ z(k - 1) * β + x(k - 2)],
+    t)
+push!(syss, discsys)
+push!(probs, DiscreteProblem(discsys, merge(Dict(unknowns(sys) .=> 0.0), Dict(u0)), p))
+
+for (sys, prob) in zip(syss, probs)
+    @test parameter_values(prob) isa ModelingToolkit.MTKParameters
+
+    @inferred typeof(prob) remake(prob; u0 = [x => 2.0], p = [σ => 29.0])
+
+    ugetter = getu(prob, [x, y, z])
+    prob2 = remake(prob; u0 = [x => 2.0, y => 3.0, z => 4.0])
+    @test ugetter(prob2) == [2.0, 3.0, 4.0]
+    prob2 = remake(prob; u0 = [sys.x => 2.0, sys.y => 3.0, sys.z => 4.0])
+    @test ugetter(prob2) == [2.0, 3.0, 4.0]
+    prob2 = remake(prob; u0 = [:x => 2.0, :y => 3.0, :z => 4.0])
+    @test ugetter(prob2) == [2.0, 3.0, 4.0]
+    prob2 = remake(prob; u0 = [x => 2.0, sys.y => 3.0, :z => 4.0])
+    @test ugetter(prob2) == [2.0, 3.0, 4.0]
+
+    prob2 = remake(prob; u0 = [x => 12.0])
+    @test ugetter(prob2) == [12.0, 0.0, 0.0]
+    prob2 = remake(prob; u0 = [sys.x => 12.0])
+    @test ugetter(prob2) == [12.0, 0.0, 0.0]
+    prob2 = remake(prob; u0 = [:x => 12.0])
+    @test ugetter(prob2) == [12.0, 0.0, 0.0]
+
+    pgetter = getp(prob, [σ, β, ρ, q])
+    prob2 = remake(prob; p = [σ => 0.1, β => 0.2, ρ => 0.3])
+    @test pgetter(prob2) == [0.1, 0.2, 0.3]
+    prob2 = remake(prob; p = [sys.σ => 0.1, sys.β => 0.2, sys.ρ => 0.3])
+    @test pgetter(prob2) == [0.1, 0.2, 0.3]
+    prob2 = remake(prob; p = [:σ => 0.1, :β => 0.2, :ρ => 0.3])
+    @test pgetter(prob2) == [0.1, 0.2, 0.3]
+    prob2 = remake(prob; p = [σ => 0.1, sys.β => 0.2, :ρ => 0.3])
+    @test pgetter(prob2) == [0.1, 0.2, 0.3]
+
+    prob2 = remake(prob; p = [σ => 0.5])
+    @test pgetter(prob2) == [0.5, 8 / 3, 10.0]
+    prob2 = remake(prob; p = [sys.σ => 0.5])
+    @test pgetter(prob2) == [0.5, 8 / 3, 10.0]
+    prob2 = remake(prob; p = [:σ => 0.5])
+    @test pgetter(prob2) == [0.5, 8 / 3, 10.0]
+end
+
+@variables ud(t) xd(t) yd(t)
+@parameters p1 p2::Int [tunable = false] p3
+dt = 0.1
+c = Clock(t, dt)
+k = ShiftIndex(c)
+
+eqs = [D(x) ~ Hold(ud)
+       ud ~ ud(k - 1) * p1 + yd
+       yd ~ p2 * yd(k - 1) + xd(k - 2) * p3
+       xd ~ Sample(t, dt)(x)]
+@mtkbuild sys = ODESystem(eqs, t; parameter_dependencies = [p3 => 2p1])
+prob = ODEProblem(sys, [x => 1.0], (0.0, 5.0),
+    [p1 => 1.0, p2 => 2, ud(k - 1) => 3.0, xd(k - 1) => 4.0, xd(k - 2) => 5.0])
+
+# parameter dependencies
+prob2 = remake(prob; p = [p1 => 2.0])
+@test prob2.ps[p1] == 2.0
+@test prob2.ps[p3] == 4.0
+@test prob2.ps[p2] isa Int # types preserved
+
+# ignore dependent parameter, preserve type
+prob2 = remake(prob; p = [p1 => 2.0, p3 => 1.0, p2 => 3])
+@test prob2.ps[p3] == 4.0
+@test prob2.ps[p2] isa Int
+
+# discrete portion
+prob2 = remake(prob; p = [ud => 2.5, xd => 3.5, xd(k - 1) => 4.5])
+@test prob2.ps[ud] == 2.5
+@test prob2.ps[xd] == 3.5
+@test prob2.ps[xd(k - 1)] == 4.5
+
+# Optimization
+@parameters p
+@mtkbuild sys = ODESystem([D(x) ~ -p * x], t)
+odeprob = ODEProblem(sys, [x => 1.0], (0.0, 10.0), [p => 0.5])
+
+ts = 0.0:0.5:10.0
+data = exp.(-2.5 .* ts)
+
+function loss(x, p)
+    prob = p[1]
+
+    prob = remake(prob; p = [prob.f.sys.p => x[1]], u0 = typeof(x)(prob.u0))
+    sol = solve(prob, Tsit5())
+    vals = sol(ts; idxs = prob.f.sys.x).u
+    return sum((data .- vals) .^ 2) / length(ts)
+end
+
+f = OptimizationFunction(loss, Optimization.AutoForwardDiff())
+prob = OptimizationProblem(f, [0.5], [odeprob])
+sol = solve(prob; BFGS())
+@test sol.u[1] ≈ 2.5
diff --git a/test/remake_tests.jl b/test/remake_tests.jl
index e5d6f6600..59aa2bba7 100644
--- a/test/remake_tests.jl
+++ b/test/remake_tests.jl
@@ -1,130 +1,206 @@
 using SciMLBase
 using SymbolicIndexingInterface
+using StaticArrays
+using ForwardDiff
 
+probs = []
+containerTypes = [Vector, Tuple, SVector{3}, MVector{3}, SizedVector{3}]
 # ODE
 function lorenz!(du, u, p, t)
     du[1] = p[1] * (u[2] - u[1])
     du[2] = u[1] * (p[2] - u[3]) - u[2]
     du[3] = u[1] * u[2] - p[3] * u[3]
 end
-u0 = [1.0; 0.0; 0.0]
+u0 = [1.0; 2.0; 3.0]
 tspan = (0.0, 100.0)
-p = [10.0, 28.0, 8 / 3]
-fn = ODEFunction(lorenz!; sys = SymbolCache([:x, :y, :z], [:a, :b, :c], :t))
-prob = ODEProblem(fn, u0, tspan, p)
-
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0, 4.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; u0 = [:x => 2.0, :z => 4.0, :y => 3.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; p = [11.0, 12.0, 13.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = [:a => 11.0, :c => 13.0, :b => 12.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = (11.0, 12.0, 13)).p == (11.0, 12.0, 13)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, 0.0, 0.0]
-@test remake(prob; p = [:b => 11.0]).p == [10.0, 11.0, 8 / 3]
-
-# BVP
-g = 9.81
-L = 1.0
-tspan = (0.0, pi / 2)
-function simplependulum!(du, u, p, t)
-    θ = u[1]
-    dθ = u[2]
-    du[1] = dθ
-    du[2] = -(p[1] / p[2]) * sin(θ)
+p = [10.0, 20.0, 30.0]
+sys = SymbolCache([:x, :y, :z], [:a, :b, :c], :t)
+fn = ODEFunction(lorenz!; sys)
+for T in containerTypes
+    push!(probs, ODEProblem(fn, u0, tspan, T(p)))
 end
-function bc1!(residual, u, p, t)
-    residual[1] = u[end ÷ 2][1] + pi / 2 # the solution at the middle of the time span should be -pi/2
-    residual[2] = u[end][1] - pi / 2 # the solution at the end of the time span should be pi/2
+
+function residual!(resid, u, p, t)
+    resid[1] = u[1] - 0.5
+    resid[2] = u[2] - 0.5
+    resid[3] = u[3] - 0.5
 end
-u0 = [pi / 2, pi / 2]
-p = [g, L]
-fn = BVPFunction(simplependulum!, bc1!; sys = SymbolCache([:x, :y], [:a, :b], :t))
-prob = BVProblem(fn, u0, tspan, p)
-
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0]).u0 == [2.0, 3.0]
-@test remake(prob; u0 = [:x => 2.0, :y => 3.0]).u0 == [2.0, 3.0]
-@test remake(prob; p = [11.0, 12.0]).p == [11.0, 12.0]
-@test remake(prob; p = [:a => 11.0, :b => 12.0]).p == [11.0, 12.0]
-@test remake(prob; p = (11.0, 12.0)).p == (11.0, 12.0)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, pi / 2]
-@test remake(prob; p = [:b => 11.0]).p == [g, 11.0]
-
-# SDE
-function sdef(du, u, p, t)
-    du[1] = p[1] * (u[2] - u[1])
-    du[2] = u[1] * (p[2] - u[3]) - u[2]
-    du[3] = u[1] * u[2] - p[3] * u[3]
+fn = BVPFunction(lorenz!, residual!; sys)
+for T in containerTypes
+    push!(probs, BVProblem(fn, u0, tspan, T(p)))
 end
-function sdeg(du, u, p, t)
+
+function noise!(du, u, p, t)
     du .= 0.1u
 end
+fn = SDEFunction(lorenz!, noise!; sys)
+for T in containerTypes
+    push!(probs, SDEProblem(fn, u0, tspan, T(p)))
+end
 
-u0 = [1.0, 0.0, 0.0]
-p = [10, 2.33, 26]
-tspan = (0, 100)
-fn = SDEFunction(sdef, sdeg; sys = SymbolCache([:x, :y, :z], [:a, :b, :c], :t))
-prob = SDEProblem(fn, u0, tspan, p)
-
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0, 4.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; u0 = [:x => 2.0, :z => 4.0, :y => 3.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; p = [11.0, 12.0, 13.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = [:a => 11.0, :c => 13.0, :b => 12.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = (11.0, 12.0, 13)).p == (11.0, 12.0, 13)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, 0.0, 0.0]
-@test remake(prob; p = [:b => 11.0]).p == [10.0, 11.0, 26.0]
-
-# OptimizationProblem
-function loss(u, p)
-    return (p[1] - u[1])^2 + p[2] * (u[2] - u[1]^2)^2
+function loss(x, p)
+    du = similar(x)
+    lorenz!(du, u, p, 0.0)
+    return sum(du)
 end
-u0 = [1.0, 2.0]
-p = [1.0, 100.0]
-fn = OptimizationFunction(loss; sys = SymbolCache([:x, :y], [:a, :b], :t))
-prob = OptimizationProblem(fn, u0, p)
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0]).u0 == [2.0, 3.0]
-@test remake(prob; u0 = [:x => 2.0, :y => 3.0]).u0 == [2.0, 3.0]
-@test remake(prob; p = [11.0, 12.0]).p == [11.0, 12.0]
-@test remake(prob; p = [:a => 11.0, :b => 12.0]).p == [11.0, 12.0]
-@test remake(prob; p = (11.0, 12.0)).p == (11.0, 12.0)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, 2.0]
-@test remake(prob; p = [:b => 11.0]).p == [1.0, 11.0]
-
-# NonlinearProblem
-function nlf(du, u, p)
-    du[1] = p[1] * (u[2] - u[1])
-    du[2] = u[1] * (p[2] - u[3]) - u[2]
-    du[3] = u[1] * u[2] - p[3] * u[3]
+
+fn = OptimizationFunction(loss; sys)
+for T in containerTypes
+    push!(probs, OptimizationProblem(fn, u0, T(p)))
+end
+
+function nllorenz!(du, u, p)
+    lorenz!(du, u, p, 0.0)
+end
+
+fn = NonlinearFunction(nllorenz!; sys)
+for T in containerTypes
+    push!(probs, NonlinearProblem(fn, u0, T(p)))
+end
+
+for T in containerTypes
+    push!(probs, NonlinearLeastSquaresProblem(fn, u0, T(p)))
+end
+
+for prob in deepcopy(probs)
+    prob2 = remake(prob)
+    @test prob2.u0 == u0
+    @test prob2.p == typeof(prob.p)(p)
+    for T in containerTypes
+        if T !== Tuple
+            local u0 = T([2.0, 3.0, 4.0])
+            prob2 = remake(prob; u0 = deepcopy(u0))
+            @test prob2.u0 == u0
+            @test prob2.u0 isa T
+        end
+        local p = T([11.0, 12.0, 13.0])
+        prob2 = remake(prob; p = deepcopy(p))
+        @test prob2.p == p
+        @test prob2.p isa T
+    end
+
+    for T in [Float32, Float64]
+        local u0 = [:x => T(2.0), :z => T(4.0), :y => T(3.0)]
+        prob2 = remake(prob; u0)
+        @test all(prob2.u0 .≈ T[2.0, 3.0, 4.0])
+        @test eltype(prob2.u0) == T
+
+        local u0 = [:x => T(2.0)]
+        prob2 = remake(prob; u0)
+        @test all(prob2.u0 .≈ [2.0, 2.0, 3.0])
+        @test eltype(prob2.u0) == Float64 # partial update promotes, since fallback is Float64
+
+        local p = [:a => T(11.0), :b => T(12.0), :c => T(13.0)]
+        prob2 = remake(prob; p)
+        @test all(prob2.p .≈ T[11.0, 12.0, 13.0])
+        @test eltype(prob2.p) == T
+
+        local p = [:a => T(11.0)]
+        prob2 = remake(prob; p)
+        @test all(prob2.p .≈ [11.0, 20.0, 30.0])
+        if prob.p isa Tuple
+            @test prob2.p isa Tuple{T, Float64, Float64}
+        else
+            @test eltype(prob2.p) == Float64
+        end
+    end
+
+    # constant defaults
+    begin
+        prob.f.sys.defaults[:a] = 0.1
+        prob.f.sys.defaults[:x] = 0.1
+        # remake with no updates should use existing values
+        prob2 = remake(prob)
+        @test prob2.u0 == u0
+        @test prob2.p == typeof(prob.p)(p)
+
+        # respect defaults (:x), fallback to existing value (:z)
+        prob2 = remake(prob; u0 = [:y => 0.2])
+        @test prob2.u0 ≈ [0.1, 0.2, 3.0]
+        @test prob2.p == typeof(prob.p)(p) # params unaffected
+
+        # override defaults
+        prob2 = remake(prob; u0 = [:x => 0.2])
+        @test prob2.u0 ≈ [0.2, 2.0, 3.0]
+        @test prob2.p == typeof(prob.p)(p)
+
+        prob2 = remake(prob; p = [:b => 0.2])
+        @test prob2.u0 == u0
+        @test all(prob2.p .≈ [0.1, 0.2, 30.0])
+
+        prob2 = remake(prob; p = [:a => 0.2])
+        @test prob2.u0 == u0
+        @test all(prob2.p .≈ [0.2, 20.0, 30.0])
+
+        empty!(prob.f.sys.defaults)
+    end
+
+    # dependent defaults
+    begin
+        prob.f.sys.defaults[:b] = :(3a)
+        prob.f.sys.defaults[:y] = :(3x)
+        # remake with no updates should use existing values
+        prob2 = remake(prob)
+        @test prob2.u0 == u0
+        @test prob2.p == typeof(prob.p)(p)
+
+        # respect defaults (:y), fallback to existing value (:z)
+        prob2 = remake(prob; u0 = [:x => 0.2])
+        @test prob2.u0 ≈ [0.2, 0.6, 3.0]
+        @test prob2.p == typeof(prob.p)(p) # params unaffected
+
+        # override defaults
+        prob2 = remake(prob; u0 = [:y => 0.2])
+        @test prob2.u0 ≈ [1.0, 0.2, 3.0]
+        @test prob2.p == typeof(prob.p)(p)
+
+        prob2 = remake(prob; p = [:a => 0.2])
+        @test prob2.u0 == u0
+        @test all(prob2.p .≈ [0.2, 0.6, 30.0])
+
+        prob2 = remake(prob; p = [:b => 0.2])
+        @test prob2.u0 == u0
+        @test all(prob2.p .≈ [10.0, 0.2, 30.0])
+
+        empty!(prob.f.sys.defaults)
+    end
+
+    # defaults dependent on each other (params <-> states)
+    begin
+        prob.f.sys.defaults[:b] = :(3x)
+        prob.f.sys.defaults[:y] = :(3a)
+        # remake with no updates should use existing values
+        prob2 = remake(prob)
+        @test prob2.u0 == u0
+        @test prob2.p == typeof(prob.p)(p)
+
+        # need to pass empty `Dict()` to prevent defaulting to existing values
+        prob2 = remake(prob; u0 = [:x => 0.2], p = Dict())
+        @test prob2.u0 ≈ [0.2, 30.0, 3.0]
+        @test all(prob2.p .≈ [10.0, 0.6, 30.0])
+
+        # override defaults
+        prob2 = remake(prob; u0 = [:y => 0.2], p = Dict())
+        @test prob2.u0 ≈ [1.0, 0.2, 3.0]
+        @test all(prob2.p .≈ [10.0, 3.0, 30.0])
+
+        prob2 = remake(prob; p = [:a => 0.2], u0 = Dict())
+        @test prob2.u0 ≈ [1.0, 0.6, 3.0]
+        @test all(prob2.p .≈ [0.2, 3.0, 30.0])
+
+        prob2 = remake(prob; p = [:b => 0.2], u0 = Dict())
+        @test prob2.u0 ≈ [1.0, 30.0, 3.0]
+        @test all(prob2.p .≈ [10.0, 0.2, 30.0])
+
+        empty!(prob.f.sys.defaults)
+    end
+
+    if !isa(prob.p, Tuple)
+        function fakeloss!(p)
+            prob2 = remake(prob; p = [:a => p])
+            @test eltype(prob2.p) <: ForwardDiff.Dual
+            return prob2.ps[:a]
+        end
+        ForwardDiff.derivative(fakeloss!, 1.0)
+    end
 end
-u0 = [1.0, 0.0, 0.0]
-p = [10, 2.33, 26]
-fn = NonlinearFunction(nlf; sys = SymbolCache([:x, :y, :z], [:a, :b, :c], :t))
-prob = NonlinearProblem(fn, u0, p)
-
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0, 4.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; u0 = [:x => 2.0, :z => 4.0, :y => 3.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; p = [11.0, 12.0, 13.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = [:a => 11.0, :c => 13.0, :b => 12.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = (11.0, 12.0, 13)).p == (11.0, 12.0, 13)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, 0.0, 0.0]
-@test remake(prob; p = [:b => 11.0]).p == [10.0, 11.0, 26.0]
-
-# NonlinearLeastSquaresProblem
-prob = NonlinearLeastSquaresProblem(fn, u0, p)
-@test remake(prob).u0 == u0
-@test remake(prob).p == p
-@test remake(prob; u0 = [2.0, 3.0, 4.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; u0 = [:x => 2.0, :z => 4.0, :y => 3.0]).u0 == [2.0, 3.0, 4.0]
-@test remake(prob; p = [11.0, 12.0, 13.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = [:a => 11.0, :c => 13.0, :b => 12.0]).p == [11.0, 12.0, 13.0]
-@test remake(prob; p = (11.0, 12.0, 13)).p == (11.0, 12.0, 13)
-@test remake(prob; u0 = [:x => 2.0]).u0 == [2.0, 0.0, 0.0]
-@test remake(prob; p = [:b => 11.0]).p == [10.0, 11.0, 26.0]
diff --git a/test/runtests.jl b/test/runtests.jl
index 3fcbc36b4..3cd39c34a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -64,6 +64,10 @@ end
         @time @safetestset "Problem building tests" begin
             include("problem_building_test.jl")
         end
+    end
+
+    if !is_APPVEYOR &&
+       (GROUP == "Core" || GROUP == "All" || GROUP == "SymbolicIndexingInterface")
         @time @safetestset "Remake" begin
             include("remake_tests.jl")
         end
@@ -106,15 +110,15 @@ end
         @time @safetestset "Partial Functions" begin
             include("downstream/partial_functions.jl")
         end
-        @time @safetestset "ModelingToolkit Remake" begin
-            include("downstream/modelingtoolkit_remake.jl")
-        end
     end
 
     if !is_APPVEYOR && (GROUP == "Downstream" || GROUP == "SymbolicIndexingInterface")
         if GROUP != "Downstream"
             activate_downstream_env()
         end
+        @time @safetestset "ModelingToolkit Remake" begin
+            include("downstream/modelingtoolkit_remake.jl")
+        end
         @time @safetestset "Symbol and integer based indexing of interpolated solutions" begin
             include("downstream/symbol_indexing.jl")
         end