diff --git a/Project.toml b/Project.toml index 9e4ea76..367e6f4 100644 --- a/Project.toml +++ b/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 , Avik Pal and contributors"] [deps] diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md index 5ab7f37..2c5fd53 100644 --- a/docs/src/_changelog.md +++ b/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 diff --git a/src/common_interface.jl b/src/common_interface.jl index 095a910..25b0c4c 100644 --- a/src/common_interface.jl +++ b/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`. diff --git a/src/expsumfit.jl b/src/expsumfit.jl index 6689304..a00dda0 100644 --- a/src/expsumfit.jl +++ b/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( diff --git a/src/king.jl b/src/king.jl index aaac21f..ea63d0d 100644 --- a/src/king.jl +++ b/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) diff --git a/src/linfit.jl b/src/linfit.jl index 68c7ea8..36701c3 100644 --- a/src/linfit.jl +++ b/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 diff --git a/src/nonlinfit.jl b/src/nonlinfit.jl index 0b132ea..6408092 100644 --- a/src/nonlinfit.jl +++ b/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 diff --git a/src/rationalfit.jl b/src/rationalfit.jl index 64c3679..adfc6a2 100644 --- a/src/rationalfit.jl +++ b/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) diff --git a/src/stats.jl b/src/stats.jl index 5ec0c39..5cbad83 100644 --- a/src/stats.jl +++ b/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) diff --git a/test/linfit_polynomial.jl b/test/linfit_polynomial.jl index 31cf99b..118080a 100644 --- a/test/linfit_polynomial.jl +++ b/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) diff --git a/test/nonlinfit_reinit.jl b/test/nonlinfit_reinit.jl index d4a24ee..e9283d3 100644 --- a/test/nonlinfit_reinit.jl +++ b/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 diff --git a/test/nonlinfit_weighted.jl b/test/nonlinfit_weighted.jl index c49fdc8..9665959 100644 --- a/test/nonlinfit_weighted.jl +++ b/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 diff --git a/test/stats.jl b/test/stats.jl index 6bd4658..031da6a 100644 --- a/test/stats.jl +++ b/test/stats.jl @@ -193,21 +193,18 @@ using Distributions: TDist, quantile @test all(stderror(sol_nc)[2:end] .> 0) @test confint(sol_nc)[1] == (0.0, 0.0) - # Modified King Fit (E^2 = A + B * U^n) - # x corresponds to E (Voltage) in Jacobian logic, but input data order is (U, E^2)? - # User creates CurveFitProblem(x, y). - # King assumes (Voltage, Velocity). - # My implementation: x=E, y=U in Jacobian. - # But in verify script I passed (U, E^2). - # So I will replicate verify script logic here. + # Modified King Fit: E² = A + B·Uⁿ, with x = E (voltage), y = U (velocity) A_true = 1.0; B_true = 2.0; n_true = 0.5 - x_k = collect(1.0:0.5:5.0) - y_k = @. A_true + B_true * x_k^n_true - prob_mod_king = CurveFitProblem(x_k, y_k) + U = collect(1.0:0.5:5.0) + E = sqrt.(A_true .+ B_true .* U .^ n_true) + prob_mod_king = CurveFitProblem(E, U) sol_mod_king = solve(prob_mod_king, ModifiedKingCurveFitAlgorithm()) + @test sol_mod_king.u ≈ [A_true, B_true, n_true] @test size(vcov(sol_mod_king)) == (3, 3) @test all(stderror(sol_mod_king) .> 0) + # Residuals are reported in velocity space, consistent with fitted() + @test residuals(sol_mod_king) ≈ U .- fitted(sol_mod_king) end @@ -288,6 +285,28 @@ using Distributions: TDist, quantile @test vcov(sol) ≈ vcov(sol; absolute_sigma = true) .* mse(sol) end + @testset "Residual-style problems (y = nothing)" begin + x = collect(1.0:10.0) + y = @. 2.0 + 3.0 * x + 0.1 * sin(x) + + f_oop(u, x) = @. u[1] + u[2] * x - $y + f_iip(resid, u, x) = @. resid = u[1] + u[2] * x - $y + fn_iip = NonlinearFunction{true}(f_iip; resid_prototype = similar(x)) + + # Smoke test that the stats functions handle problems without a `y` + for f in (f_oop, fn_iip) + sol = solve(NonlinearCurveFitProblem(f, [1.0, 1.0], x)) + + @test nobs(sol) == length(x) + @test dof_residual(sol) == length(x) - 2 + @test residuals(sol) == sol.resid + @test predict(sol) ≈ residuals(sol) + @test mse(sol) ≈ rss(sol) / dof_residual(sol) + @test length(stderror(sol)) == 2 + @test all(lo < c < hi for (c, (lo, hi)) in zip(coef(sol), confint(sol))) + end + end + @testset "Extended Algorithm Test Coverage" begin x_data = collect(1.0:0.5:5.0) y_data = 2.0 .* x_data .^ 0.5 .+ 0.1