fix: call u0_constructor on resid_prototype

Author: AayushSabharwal

src/systems/nonlinear/nonlinearsystem.jl

@@ -345,16 +345,6 @@ function hessian_sparsity(sys::NonlinearSystem)
          unknowns(sys)) for eq in equations(sys)]
 end
 
-function calculate_resid_prototype(N, u0, p)
-    u0ElType = u0 === nothing ? Float64 : eltype(u0)
-    if SciMLStructures.isscimlstructure(p)
-        u0ElType = promote_type(
-            eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
-            u0ElType)
-    end
-    return zeros(u0ElType, N)
-end
-
 """
 ```julia
 SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(sys),
@@ -381,6 +371,7 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
         eval_module = @__MODULE__,
         sparse = false, simplify = false,
         initialization_data = nothing, cse = true,
+        resid_prototype = nothing,
         kwargs...) where {iip}
     if !iscomplete(sys)
         error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearFunction`")
@@ -402,12 +393,6 @@ function SciMLBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = unknowns(s
     observedfun = ObservedFunctionCache(
         sys; eval_expression, eval_module, checkbounds = get(kwargs, :checkbounds, false), cse)
 
-    if length(dvs) == length(equations(sys))
-        resid_prototype = nothing
-    else
-        resid_prototype = calculate_resid_prototype(length(equations(sys)), u0, p)
-    end
-
     NonlinearFunction{iip}(f;
         sys = sys,
         jac = _jac === nothing ? nothing : _jac,

src/systems/problem_utils.jl

@@ -1081,6 +1081,22 @@ function float_type_from_varmap(varmap, floatT = Bool)
     return float(floatT)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Calculate the `resid_prototype` for a `NonlinearFunction` with `N` equations and the
+provided `u0` and `p`.
+"""
+function calculate_resid_prototype(N::Int, u0, p)
+    u0ElType = u0 === nothing ? Float64 : eltype(u0)
+    if SciMLStructures.isscimlstructure(p)
+        u0ElType = promote_type(
+            eltype(SciMLStructures.canonicalize(SciMLStructures.Tunable(), p)[1]),
+            u0ElType)
+    end
+    return zeros(u0ElType, N)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -1292,7 +1308,14 @@ function process_SciMLProblem(
         end
         initialization_data = SciMLBase.remake_initialization_data(
             kwargs.initialization_data, kwargs, u0, t0, p, u0, p)
-        kwargs = merge(kwargs,)
+        kwargs = merge(kwargs, (; initialization_data))
+    end
+
+    if constructor <: NonlinearFunction && length(dvs) != length(eqs)
+        kwargs = merge(kwargs,
+            (;
+                resid_prototype = u0_constructor(calculate_resid_prototype(
+                    length(eqs), u0, p))))
     end
 
     f = constructor(sys, dvs, ps, u0; p = p,

test/nonlinearsystem.jl

@@ -442,3 +442,19 @@ end
         @test !in(D(y), vs)
     end
 end
+
+@testset "oop `NonlinearLeastSquaresProblem` with `u0 === nothing`" begin
+    @variables x y
+    @named sys = NonlinearSystem([0 ~ x - y], [], []; observed = [x ~ 1.0, y ~ 1.0])
+    prob = NonlinearLeastSquaresProblem{false}(complete(sys), nothing)
+    sol = solve(prob)
+    resid = sol.resid
+    @test resid == [0.0]
+    @test resid isa Vector
+    prob = NonlinearLeastSquaresProblem{false}(
+        complete(sys), nothing; u0_constructor = splat(SVector))
+    sol = solve(prob)
+    resid = sol.resid
+    @test resid == [0.0]
+    @test resid isa SVector
+end