Fix residuals space of modified King fits and add reinit!() checks

Project.toml

@@ -1,6 +1,6 @@
 name = "CurveFit"
 uuid = "5a033b19-8c74-5913-a970-47c3779ef25c"
-version = "1.10.0"
+version = "1.10.1"
 authors = ["Paulo José Saiz Jabardo <pjabardo@gmail.com>, Avik Pal <avikpal@mit.edu> and contributors"]
 
 [deps]

docs/src/_changelog.md

@@ -7,6 +7,25 @@ CurrentModule = CurveFit
 This documents notable changes in CurveFit.jl. The format is based on [Keep a
 Changelog](https://keepachangelog.com).
 
+## [v1.10.1] - 2026-07-01
+
+### Fixed
+- Fixed the residuals space of [`ModifiedKingCurveFitAlgorithm`](@ref),
+  previously they were in E² space but now they're calculated in velocity space
+  such that they're consistent with the standard definition of `residuals = y -
+  ŷ` ([#116]).
+- The `reinit!()` method for [`NonlinearCurveFitProblem`](@ref)'s will now
+  explicitly throw an error when invalid inputs are given (e.g. arrays of the
+  wrong length) to prevent internal corruption ([#116]).
+- Fixed the residuals and statistics calculations for polynomial fits,
+  previously they were computed from corrupted data ([#116]).
+- Solving with [`ExpSumFitAlgorithm`](@ref) will no longer modify the input x/y
+  arrays, and will return a vector of residuals instead of a `(n, 1)` matrix
+  ([#116]).
+- Fixed various statistics methods when using residual-only functions with
+  [`NonlinearCurveFitProblem`](@ref), as well as adding support for `sigma` with
+  such such functions ([#116])
+
 ## [v1.10.0] - 2026-06-28
 
 ### Added

src/common_interface.jl

@@ -99,8 +99,10 @@ optimization problem is solved:
 \argmin_u ~ \left\| f(u, x) - y \right\|_2
 ```
 
-If `y` is `nothing`, then it is treated as a zero vector. `f` is a generic Julia function or
-ideally a `NonlinearFunction` from [`SciMLBase.jl`](https://github.com/SciML/SciMLBase.jl).
+If `y` is `nothing`, then it is treated as a zero vector and `f` is expected to
+return residuals directly. `f` is a generic Julia function or ideally a
+`NonlinearFunction` from
+[`SciMLBase.jl`](https://github.com/SciML/SciMLBase.jl).
 
 Lower and upper bounds on the parameters can be passed through `lb` and `ub` keyword
 arguments. These should be arrays with the same length as `u0`.

src/expsumfit.jl

@@ -137,6 +137,9 @@ end
     prob <: CurveFitProblem
     alg
     kwargs
+    # Scaled working copies of prob.x/prob.y, the inputs are never modified
+    x <: AbstractVector
+    y <: AbstractVector
     Y <: AbstractMatrix
     S <: AbstractMatrix
     A <: AbstractVector
@@ -162,6 +165,8 @@ function CommonSolve.init(
         prob,
         alg,
         kwargs,
+        similar(prob.x),
+        similar(prob.y),
         similar(prob.x, nY, mY),
         similar(prob.x, nY, alg.n - 1),
         similar(prob.x, mY),
@@ -174,28 +179,25 @@ end
 
 # TODO: allocations in this function aren't optimized
 function CommonSolve.solve!(cache::ExpSumFitCache)
-    sc = __expsum_scale!(cache.prob.x, cache.prob.y)
+    copyto!(cache.x, cache.prob.x)
+    copyto!(cache.y, cache.prob.y)
+    sc = __expsum_scale!(cache.x, cache.y)
 
     __expsum_fill_Y!(
-        cache.Y, cache.S, cache.coeff, cache.prob.x, cache.prob.y, cache.alg.n, cache.alg.m
+        cache.Y, cache.S, cache.coeff, cache.x, cache.y, cache.alg.n, cache.alg.m
     )
 
     λ = __expsum_solve_λ(
-        cache.Y, cache.A, cache.Ā, cache.prob.y, cache.alg.n, cache.alg.m
+        cache.Y, cache.A, cache.Ā, cache.y, cache.alg.n, cache.alg.m
     )
 
     X = isreal(λ) ? cache.Xr : cache.Xc
-    __expsum_fill_X!(cache.prob.x, λ, X, cache.alg.n)
+    __expsum_fill_X!(cache.x, λ, X, cache.alg.n)
 
     qrX = qr!(X)
-    p = qrX \ cache.prob.y
+    p = qrX \ cache.y
     isreal(p) && (p = real(p))
 
-    for i in eachindex(cache.prob.x, cache.prob.y)
-        cache.prob.x[i] *= sc.x
-        cache.prob.y[i] *= sc.y
-    end
-
     for i in eachindex(p)
         p[i] *= sc.y
     end
@@ -214,7 +216,7 @@ function CommonSolve.solve!(cache::ExpSumFitCache)
     end
     backing = (; k, p, λ)
 
-    y_pred = k .+ sum(exp.(cache.prob.x * λ') .* p'; dims = 2)
+    y_pred = k .+ vec(sum(exp.(cache.prob.x * λ') .* p'; dims = 2))
     resid = cache.prob.y .- y_pred
 
     return CurveFitSolution(

src/king.jl

@@ -39,7 +39,6 @@ end
 # Common Solve Interface for ModifiedKingCurveFitAlgorithm
 @concrete struct ModifiedKingFitCache <: AbstractCurveFitCache
     initial_guess_cache <: Union{Nothing, KingFitCache}
-    nonlinear_cache
     prob <: CurveFitProblem
     alg <: ModifiedKingCurveFitAlgorithm
     kwargs
@@ -66,21 +65,7 @@ function CommonSolve.init(
         init(nobounds_prob, KingCurveFitAlgorithm(); kwargs...)
     end
 
-    nonlinear_cache = init(
-        NonlinearCurveFitProblem(
-            NonlinearFunction{true}(
-                __king_fun!;
-                resid_prototype = similar(prob.x)
-            ),
-            similar(prob.x, 3),
-            stack((prob.x, prob.y); dims = 1),
-            nothing;
-            lb = prob.lb, ub = prob.ub
-        ),
-        __FallbackNonlinearFitAlgorithm(alg.alg);
-        kwargs...
-    )
-    return ModifiedKingFitCache(initial_guess_cache, nonlinear_cache, prob, alg, kwargs)
+    return ModifiedKingFitCache(initial_guess_cache, prob, alg, kwargs)
 end
 
 function CommonSolve.solve!(cache::ModifiedKingFitCache)
@@ -105,7 +90,15 @@ function CommonSolve.solve!(cache::ModifiedKingFitCache)
     )
 
     sol = solve(nonlinear_prob, __FallbackNonlinearFitAlgorithm(cache.alg.alg); cache.kwargs...)
-    return CurveFitSolution(cache.alg, sol.u, sol.resid, cache.prob, sol.retcode, sol.original)
+
+    # The inner solve fits E² = A + B·Uⁿ, so its residual `A + B·Uⁿ − E²` is a
+    # difference of E² values. We instead report `residuals` as a difference of
+    # velocities, `U − Û`, so they match `fitted`/`sol(x)` (i.e.
+    # `residuals == prob.y .- fitted`). `vcov` rebuilds the E² residual it needs
+    # from `sol.u`.
+    A, B, n = sol.u
+    resid = @. cache.prob.y - ((cache.prob.x^2 - A) / B)^(1 / n)
+    return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode, sol.original)
 end
 
 function (sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm})(x)

src/linfit.jl

@@ -106,8 +106,9 @@ function CommonSolve.solve!(cache::PolynomialFitCache)
     __vandermondepoly!(cache.vandermondepoly_cache, cache.prob.x, cache.alg.degree)
     cache.linsolve_cache.A = cache.vandermondepoly_cache
     sol = solve!(cache.linsolve_cache)
-    # Always compute residuals manually as LinearSolve may return incorrect residuals
-    resid = cache.prob.y .- cache.vandermondepoly_cache * sol.u
+    # Compute residuals from the fitted polynomial. We can't reuse
+    # vandermondepoly_cache here as the in-place factorization overwrites it.
+    resid = cache.prob.y .- evalpoly.(cache.prob.x, Ref(sol.u))
     return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode)
 end
 

src/nonlinfit.jl

@@ -11,17 +11,22 @@ _unwrap_nonlinear_function(f::NonlinearSolveBase.AutoSpecializeCallable) = _unwr
 _unwrap_nonlinear_function(f::NonlinearSolveBase.BoundedWrapper) = _unwrap_nonlinear_function(f.f.f)
 _unwrap_nonlinear_function(f) = f
 
-# If `target` is nothing then we can completely ignore sigma
-__wrap_nonlinear_function(f::NonlinearFunction, ::Nothing, _) = f
+__wrap_nonlinear_function(f::NonlinearFunction, ::Nothing, ::Nothing) = f
 function __wrap_nonlinear_function(f::NonlinearFunction, target, sigma)
     internal_f = NonlinearFunctionWrapper{SciMLBase.isinplace(f)}(target, sigma, f.f)
     @set! f.f = internal_f
-    @set! f.resid_prototype = similar(target)
+    @set! f.resid_prototype = similar(isnothing(target) ? sigma : target)
     return f
 end
 
 # Out-of-place
 function (nlf::NonlinearFunctionWrapper{false})(p, X)
+    # `target` is only nothing for residual-style problems (y = nothing), in
+    # which case the wrapper exists purely to apply `sigma`.
+    if isnothing(nlf.target)
+        return nlf.f(p, X) ./ nlf.sigma
+    end
+
     resid = nlf.target .- nlf.f(p, X)
 
     if !isnothing(nlf.sigma)
@@ -34,8 +39,10 @@ end
 # In-place
 function (nlf::NonlinearFunctionWrapper{true})(resid, p, X)
     nlf.f(resid, p, X)
-    resid .= nlf.target .- resid
 
+    if !isnothing(nlf.target)
+        resid .= nlf.target .- resid
+    end
     if !isnothing(nlf.sigma)
         resid ./= nlf.sigma
     end
@@ -54,23 +61,29 @@ end
 
 function SciMLBase.reinit!(cache::GenericNonlinearCurveFitCache; u0 = nothing, x = nothing, y = nothing, sigma = nothing, kwargs...)
     if !isnothing(u0)
+        @assert size(u0) == size(cache.u0) "reiniting `u0` must keep the same size"
         kwargs = (; kwargs..., u0)
         copyto!(cache.u0, u0)
     end
 
     # x becomes `p` (parameter) in the NonlinearLeastSquaresProblem
     if !isnothing(x)
+        @assert size(x) == size(_get_cache(cache).p) "reiniting `x` must keep the same size"
         kwargs = (; kwargs..., p = x)
     end
 
     # Update `y` inplace
     wrapper = _unwrap_nonlinear_function(cache.cache.prob.f.f)
     if !isnothing(y)
+        @assert wrapper isa NonlinearFunctionWrapper && !isnothing(wrapper.target) "cannot reinit `y` for a problem created without a `y`"
+        @assert size(y) == size(wrapper.target) "reiniting `y` must keep the same size"
         copyto!(wrapper.target, y)
     end
 
     # Update `sigma` inplace
     if !isnothing(sigma)
+        @assert wrapper isa NonlinearFunctionWrapper && !isnothing(wrapper.sigma) "cannot reinit `sigma` for a problem created without a `sigma`"
+        @assert size(sigma) == size(wrapper.sigma) "reiniting `sigma` must keep the same size"
         copyto!(wrapper.sigma, sigma)
     end
 
@@ -133,6 +146,11 @@ function CommonSolve.solve!(cache::GenericNonlinearCurveFitCache)
 end
 
 function (sol::CurveFitSolution{<:__FallbackNonlinearFitAlgorithm})(x)
+    if SciMLBase.isinplace(sol.prob) && x isa AbstractArray
+        out = similar(sol.resid, length(x))
+        sol.prob.nlfunc(out, sol.u, x)
+        return out
+    end
     return sol.prob.nlfunc(sol.u, x)
 end
 

src/rationalfit.jl

@@ -48,7 +48,6 @@ end
 
 @concrete struct NonlinearRationalFitCache <: AbstractCurveFitCache
     initial_guess_cache <: Union{Nothing, LinearRationalFitCache}
-    nonlinear_cache
     prob <: CurveFitProblem
     alg <: RationalPolynomialFitAlgorithm
     kwargs
@@ -82,23 +81,7 @@ function CommonSolve.init(
             A, init(LinearProblem(A, prob.y), alg.alg; kwargs...), prob, alg, kwargs
         )
     end
-    nonlinear_cache = init(
-        NonlinearCurveFitProblem(
-            NonlinearFunction{true}(
-                __rational_fit_residual!(alg.num_degree, alg.den_degree);
-                resid_prototype = similar(prob.x)
-            ),
-            similar(prob.x, coeffs_length),
-            prob.x,
-            prob.y;
-            lb = prob.lb, ub = prob.ub
-        ),
-        __FallbackNonlinearFitAlgorithm(alg.alg);
-        kwargs...
-    )
-    return NonlinearRationalFitCache(
-        initial_guess_cache, nonlinear_cache, prob, alg, kwargs
-    )
+    return NonlinearRationalFitCache(initial_guess_cache, prob, alg, kwargs)
 end
 
 function CommonSolve.solve!(cache::LinearRationalFitCache)

src/stats.jl

@@ -53,7 +53,7 @@ end
 Return the number of observations used in the fit.
 """
 function StatsAPI.nobs(sol::CurveFitSolution)
-    return length(sol.prob.y)
+    return length(sol.resid)
 end
 
 # Return the positions in `sol.u` that are held fixed. These coefficients
@@ -226,19 +226,6 @@ function jacobian(sol::CurveFitSolution{<:ExpSumFitAlgorithm})
     return DifferentiationInterface.jacobian(f_pred, AutoForwardDiff(), u)
 end
 
-function jacobian(sol::CurveFitSolution{<:KingCurveFitAlgorithm})
-    # King's law: U = t² where t = (E²−A)/B, so ∂U/∂A = −2t/B and ∂U/∂B = −2t²/B
-    A, B = sol.u
-    x = sol.prob.x
-    J = Matrix{eltype(x)}(undef, length(x), 2)
-    @inbounds for i in eachindex(x)
-        t = (x[i]^2 - A) / B
-        J[i, 1] = -2t / B    # ∂U/∂A
-        J[i, 2] = -2t^2 / B  # ∂U/∂B
-    end
-    return J
-end
-
 function jacobian(sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm})
     # Modified King: E^2 = A + B * U^n
     # x corresponds to E (Voltage) - Model predicts x^2
@@ -443,6 +430,24 @@ function StatsAPI.vcov(sol::CurveFitSolution{<:KingCurveFitAlgorithm}; absolute_
     return _transformed_ols_vcov(t, Y, b, a, param_map, absolute_sigma)
 end
 
+# Modified King fits E² = A + B·Uⁿ by nonlinear least squares, so its covariance
+# is the Gauss–Newton form `σ²·(JᵀJ)⁻¹` built from the E² residual it minimized.
+# `residuals`/`mse` are reported as velocities, so we rebuild that E² residual
+# here instead of reusing `mse(sol)`.
+function StatsAPI.vcov(sol::CurveFitSolution{<:ModifiedKingCurveFitAlgorithm}; absolute_sigma::Bool = false)
+    covar = _vcov_from_jacobian(jacobian(sol))
+
+    if !absolute_sigma
+        A, B, n = sol.u
+        E, U = sol.prob.x, sol.prob.y
+        resid_fitspace = @. A + B * U^n - E^2
+        σ2 = sum(abs2, resid_fitspace) / dof_residual(sol)
+        covar .*= σ2
+    end
+
+    return covar
+end
+
 """
     stderror(sol::CurveFitSolution; absolute_sigma = false, rtol = NaN, atol = 0)
 

test/linfit_polynomial.jl

@@ -21,6 +21,8 @@ using LinearSolve
         @test sol(val) ≈ fn(val)
     end
 
+    @test sol.resid ≈ y .- sol.(x)
+
     @testset "ill-conditioned" begin
         true_coeffs = [80.0, -5.0e-18, -7.0e-20, -1.0e-36]
         x1 = 1.0e10 .* (0:0.1:5)

test/nonlinfit_reinit.jl

@@ -49,4 +49,22 @@ using NonlinearSolveBase: NonlinearSolveBase
     @test cache.cache.prob.f.f isa NonlinearSolveBase.AutoSpecializeCallable
     CurveFit.reinit!(cache; u0 = [1.0, 1.0, 1.0], x, y)
     @test solve!(cache).u ≈ a0 atol = 1.0e-7
+
+    # Test reinit!() input validation with an out-of-place function so the
+    # `y=nothing` cache can be constructed.
+    g(a, x) = @. a[1] + a[2] * x^a[3]
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x, y, sigma))
+
+    # Size changes are rejected
+    @test_throws AssertionError CurveFit.reinit!(cache; u0 = [1.0, 1.0])
+    @test_throws AssertionError CurveFit.reinit!(cache; x = 1.0:5.0)
+    @test_throws AssertionError CurveFit.reinit!(cache; y = y[1:5])
+    @test_throws AssertionError CurveFit.reinit!(cache; sigma = sigma[1:5])
+
+    # Can't reinit a variable the problem was created without
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x, y))
+    @test_throws AssertionError CurveFit.reinit!(cache; sigma)
+
+    cache = CurveFit.init(NonlinearCurveFitProblem(g, [0.5, 0.5, 0.5], x))
+    @test_throws AssertionError CurveFit.reinit!(cache; y)
 end

test/nonlinfit_weighted.jl

@@ -28,3 +28,25 @@ using SciMLBase
     @test err_weighted < err_no_weight
     @test SciMLBase.successful_retcode(sol_weighted)
 end
+
+@testset "Weighted residual-style problem (y = nothing)" begin
+    x = collect(1.0:10.0)
+    y = @. 2.0 + 3.0 * x + 0.1 * sin(x)
+    sigma = collect(range(0.5, 2.0; length = length(x)))
+    u0 = [1.0, 1.0]
+
+    pred(u, x) = @. u[1] + u[2] * x
+    sol_pred = solve(NonlinearCurveFitProblem(pred, u0, x, y, sigma))
+
+    # Residual-style problems (y = nothing) must apply sigma the same way as
+    # the equivalent prediction-style problem.
+    resid_oop(u, x) = y .- pred(u, x)
+    sol_oop = solve(NonlinearCurveFitProblem(resid_oop, u0, x, nothing, sigma))
+    @test sol_oop.u ≈ sol_pred.u
+    @test sol_oop.resid ≈ sol_pred.resid
+
+    resid_iip(resid, u, x) = resid .= y .- pred(u, x)
+    sol_iip = solve(NonlinearCurveFitProblem(resid_iip, u0, x, nothing, sigma))
+    @test sol_iip.u ≈ sol_pred.u
+    @test sol_iip.resid ≈ sol_pred.resid
+end