From 5d6afa4f552d7c5bc5b8579a70f65cfa14f351a3 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Wed, 24 Jun 2026 15:22:58 +0200 Subject: [PATCH 1/4] Always apply `yfun_inverse` to linear fit intercepts Otherwise `yfun(b)` will be returned instead of `b`, which is confusing (e.g. `log(b)` instead of `b` for exp/power fits) --- docs/src/_changelog.md | 13 ++++++++++++- src/common_interface.jl | 4 ++++ src/linfit.jl | 5 +++-- test/linfit_exp.jl | 2 +- test/linfit_power.jl | 2 +- 5 files changed, 21 insertions(+), 5 deletions(-) diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md index 0118379a..c88919f9 100644 --- a/docs/src/_changelog.md +++ b/docs/src/_changelog.md @@ -7,7 +7,18 @@ CurrentModule = CurveFit This documents notable changes in CurveFit.jl. The format is based on [Keep a Changelog](https://keepachangelog.com). -## v[1.9.1] - 2026-04-25 +## Unreleased + +### Fixed +- Fixed returned parameters of [`ExpCurveFitAlgorithm`](@ref) and + [`PowerCurveFitAlgorithm`](@ref) so they aren't log-transformed ([#112]). + +## [v1.9.2] - 2026-06-24 + +### Changed +- Various improvements to CI. + +## [v1.9.1] - 2026-04-25 ### Fixed - Fixed support for `reinit!()`'ing nonlinear fit caches when using the new diff --git a/src/common_interface.jl b/src/common_interface.jl index f0d938ff..fa3866c2 100644 --- a/src/common_interface.jl +++ b/src/common_interface.jl @@ -223,6 +223,10 @@ be converted to a linear fit in a specific function space by choosing appropriat `xfun` and `yfun`. The `yfun_inverse` is used to convert the fitted values back to the original space (can be specified by defining `InverseFunctions.inverse`). +When `yfun` is not `identity`, `yfun_inverse` is applied to the fitted intercept before +storing it, so that `sol.u = (a, b)` directly satisfies the original-space formula +(e.g. `y = b * exp(a*x)` for [`ExpCurveFitAlgorithm`](@ref)). + This algorithm does not support bounds constraints (`lb`/`ub`). """ function LinearCurveFitAlgorithm(; diff --git a/src/linfit.jl b/src/linfit.jl index 46e365d1..68c7ea8a 100644 --- a/src/linfit.jl +++ b/src/linfit.jl @@ -65,14 +65,15 @@ function CommonSolve.solve!(cache::GenericLinearFitCache) resid ./= cache.prob.sigma end + b_stored = cache.alg.yfun_inverse(b) return CurveFitSolution( - cache.alg, (a, b), resid, cache.prob, ReturnCode.Success + cache.alg, (a, b_stored), resid, cache.prob, ReturnCode.Success ) end function (sol::CurveFitSolution{<:LinearCurveFitAlgorithm})(x) a, b = sol.u - return sol.alg.yfun_inverse.(b .+ a .* sol.alg.xfun.(x)) + return sol.alg.yfun_inverse.(sol.alg.yfun(b) .+ a .* sol.alg.xfun.(x)) end # Polynomial Fit diff --git a/test/linfit_exp.jl b/test/linfit_exp.jl index 5662df89..8ff0ea6b 100644 --- a/test/linfit_exp.jl +++ b/test/linfit_exp.jl @@ -11,7 +11,7 @@ using Test sol = solve(prob, ExpCurveFitAlgorithm()) @test sol.u[1] ≈ 0.8 - @test sol.u[2] ≈ log(2.0) + @test sol.u[2] ≈ 2.0 @testset for val in (0.0, 1.5, 4.5, 10.0) @test sol(val) ≈ fn(val) diff --git a/test/linfit_power.jl b/test/linfit_power.jl index 5134c81b..a8ac4fae 100644 --- a/test/linfit_power.jl +++ b/test/linfit_power.jl @@ -11,7 +11,7 @@ using Test sol = solve(prob, PowerCurveFitAlgorithm()) @test sol.u[1] ≈ 0.8 - @test sol.u[2] ≈ log(2.0) + @test sol.u[2] ≈ 2.0 @testset for val in (0.0, 1.5, 4.5, 10.0) @test sol(val) ≈ fn(val) From fe8210a6c64238c8f65557a5c055aff769028154 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Wed, 24 Jun 2026 15:51:23 +0200 Subject: [PATCH 2/4] Fix parameter handling of King curve fits Fitting King curves is done by linearizing the model, but recovering the parameters can't be handled by `LinearCurveFitAlgorithm` applying the inverse y-transform alone because the `A` and `B` parameters of the King model are related to each other. Previously this also gave inaccurate initial guesses for `ModifiedKingCurveFitAlgorithm`. The Jacobian was also incorrect since it was computed using the linear Jacobian formula. Now there's a `KingCurveFitAlgorithm` type to dispatch on and recover the parameters, and the correct analytical Jacobian is used. --- docs/src/_changelog.md | 5 +++-- src/common_interface.jl | 4 ++-- src/king.jl | 35 +++++++++++++++++++++++++++++++++-- src/stats.jl | 13 +++++++++++++ test/king_law.jl | 3 +++ test/stats.jl | 7 +------ 6 files changed, 55 insertions(+), 12 deletions(-) diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md index c88919f9..e86c274b 100644 --- a/docs/src/_changelog.md +++ b/docs/src/_changelog.md @@ -10,8 +10,9 @@ Changelog](https://keepachangelog.com). ## Unreleased ### Fixed -- Fixed returned parameters of [`ExpCurveFitAlgorithm`](@ref) and - [`PowerCurveFitAlgorithm`](@ref) so they aren't log-transformed ([#112]). +- Fixed parameter handling of [`ExpCurveFitAlgorithm`](@ref), + [`PowerCurveFitAlgorithm`](@ref), and [`KingCurveFitAlgorithm`](@ref) so they + don't return transformed parameters ([#112]). ## [v1.9.2] - 2026-06-24 diff --git a/src/common_interface.jl b/src/common_interface.jl index fa3866c2..b9ee82e6 100644 --- a/src/common_interface.jl +++ b/src/common_interface.jl @@ -295,7 +295,7 @@ ExpCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = identity, yfun = log) Represents a king curve fitting problem where `x` and `y` are the data points to fit. This algorithm does not support passing weights through `sigma` in [`CurveFitProblem`](@ref). This algorithm does not support bounds constraints (`lb`/`ub`). -We want to solve for `a` and `b` according to original King's law (1910) that represents +We want to solve for `A` and `B` according to original King's law (1910) that represents the relationship between voltage (E) and velocity (U) in a hotwire anemometer: ```math @@ -308,7 +308,7 @@ or x^2 = A + B y^{1/2} ``` """ -KingCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt) +struct KingCurveFitAlgorithm <: AbstractCurveFitAlgorithm end @doc doc""" ModifiedKingCurveFitAlgorithm(alg::Union{Nothing, AbstractNonlinearAlgorithm} = nothing) diff --git a/src/king.jl b/src/king.jl index 7e00e088..210d149f 100644 --- a/src/king.jl +++ b/src/king.jl @@ -5,9 +5,40 @@ function __king_fun!(resid, p, x) return nothing end -# Common Solve Interface +# Common Solve Interface for KingCurveFitAlgorithm +@concrete struct KingFitCache <: AbstractCurveFitCache + prob <: CurveFitProblem + alg <: KingCurveFitAlgorithm + kwargs +end + +function CommonSolve.init(prob::CurveFitProblem, alg::KingCurveFitAlgorithm; kwargs...) + @assert !is_nonlinear_problem(prob) "King's law fitting doesn't work with nlfunc specification." + @assert prob.u0 === nothing "King's law fit doesn't support initial guess (u0) specification" + sigma_not_supported(prob) + bounds_not_supported(prob) + return KingFitCache(prob, alg, kwargs) +end + +function CommonSolve.solve!(cache::KingFitCache) + # Fit sqrt(U) = b + a*E^2, then recover King's constants A and B: + # sqrt(U) = -A/B + (1/B)*E^2 => B = 1/a, A = -b/a + b, a = __linear_fit_internal(abs2, cache.prob.x, sqrt, cache.prob.y, nothing) + B = 1 / a + A = -b * B + y_pred = ((cache.prob.x .^ 2 .- A) ./ B) .^ 2 + resid = cache.prob.y .- y_pred + return CurveFitSolution(cache.alg, (A, B), resid, cache.prob, ReturnCode.Success) +end + +function (sol::CurveFitSolution{<:KingCurveFitAlgorithm})(x::Number) + A, B = sol.u + return ((x^2 - A) / B)^2 +end + +# Common Solve Interface for ModifiedKingCurveFitAlgorithm @concrete struct ModifiedKingFitCache <: AbstractCurveFitCache - initial_guess_cache <: Union{Nothing, GenericLinearFitCache} + initial_guess_cache <: Union{Nothing, KingFitCache} nonlinear_cache prob <: CurveFitProblem alg <: ModifiedKingCurveFitAlgorithm diff --git a/src/stats.jl b/src/stats.jl index c44cd8c6..fd0e222d 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -215,6 +215,19 @@ 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 diff --git a/test/king_law.jl b/test/king_law.jl index 2f7f24e8..e7476aa6 100644 --- a/test/king_law.jl +++ b/test/king_law.jl @@ -13,6 +13,9 @@ using Test prob = CurveFitProblem(E, U) sol = solve(prob, KingCurveFitAlgorithm()) + @test sol.u[1] ≈ A + @test sol.u[2] ≈ B + @testset for val in range(minimum(E), stop = maximum(E), length = 10) @test sol(val) ≈ fn(val) end diff --git a/test/stats.jl b/test/stats.jl index fc2d6e8d..207480c6 100644 --- a/test/stats.jl +++ b/test/stats.jl @@ -245,12 +245,7 @@ using LinearSolve x_data = collect(1.0:0.5:5.0) y_data = 2.0 .* x_data .^ 0.5 .+ 0.1 - # 1. King Fit (Linear alias) - # y = A + B*sqrt(x) -> Linear fit with yfun=sqrt? No, King is xfun=abs2, yfun=sqrt? - # Wait, King is LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt) - # Check source: KingCurveFitAlgorithm() = LinearCurveFitAlgorithm(; xfun = abs2, yfun = sqrt) - # This implies sqrt(y) = A + B*x^2. - # This is strictly linear in transformed space, so vcov should work. + # 1. King Fit prob_king = CurveFitProblem(x_data, y_data) sol_king = solve(prob_king, KingCurveFitAlgorithm()) @test size(vcov(sol_king)) == (2, 2) From 0fea6a43cc68ee0f55d9755e04682c5086278a08 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Wed, 24 Jun 2026 16:11:00 +0200 Subject: [PATCH 3/4] Always compute rational polynomial fits in the original space The old branch would have optionally taken the residuals from the liinearized model, which would not be correct. --- src/rationalfit.jl | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/rationalfit.jl b/src/rationalfit.jl index 18cd3a3a..cbff63b6 100644 --- a/src/rationalfit.jl +++ b/src/rationalfit.jl @@ -107,23 +107,19 @@ function CommonSolve.solve!(cache::LinearRationalFitCache) ) cache.linsolve_cache.A = cache.mat sol = solve!(cache.linsolve_cache) - resid = if sol.resid === nothing - # Linear problem: y ≈ p/q => y*q - p ≈ 0 (linearized residual) - # But StatsAPI expects y - p/q (nonlinear residual) or linearized? - # Standard definition is y - y_pred. - # So we should compute y - p(x)/q(x) using the fitted params. - - # We need to construct the RationalPolynomial to eval it. - # Helper function from rationalfit.jl isn't easily accessible inside solve! without alloc. - # But we can reuse the logic from call: - rpoly = RationalPolynomial( - view(sol.u, 1:(cache.alg.num_degree + 1)), - vcat(one(eltype(sol.u)), view(sol.u, (cache.alg.num_degree + 2):length(sol.u))) - ) - cache.prob.y .- rpoly.(cache.prob.x) - else - sol.resid - end + # Linear problem: y ≈ p/q => y*q - p ≈ 0 (linearized residual) + # But StatsAPI expects y - p/q (nonlinear residual) or linearized? + # Standard definition is y - y_pred. + # So we should compute y - p(x)/q(x) using the fitted params. + + # We need to construct the RationalPolynomial to eval it. + # Helper function from rationalfit.jl isn't easily accessible inside solve! without alloc. + # But we can reuse the logic from call: + rpoly = RationalPolynomial( + view(sol.u, 1:(cache.alg.num_degree + 1)), + vcat(one(eltype(sol.u)), view(sol.u, (cache.alg.num_degree + 2):length(sol.u))) + ) + resid = cache.prob.y .- rpoly.(cache.prob.x) return CurveFitSolution(cache.alg, sol.u, resid, cache.prob, sol.retcode) end From 110b70d220e6945f90b9760b8181a0711f725217 Mon Sep 17 00:00:00 2001 From: JamesWrigley Date: Wed, 24 Jun 2026 16:25:07 +0200 Subject: [PATCH 4/4] Fix the Jacobian calculations for linear fits with a y-transform The residuals of a linear fit are stored in original y-space, so if a transform was applied to y we would end up computing the Jacobian for a linear model on them. --- Project.toml | 2 +- docs/src/_changelog.md | 16 ++++++++++++---- src/stats.jl | 23 ++++++++++++++++++----- test/stats.jl | 26 ++++++++++++++++++++++++++ 4 files changed, 57 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 87e28cbd..68b80a37 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "CurveFit" uuid = "5a033b19-8c74-5913-a970-47c3779ef25c" -version = "1.9.2" +version = "1.9.3" authors = ["Paulo José Saiz Jabardo , Avik Pal and contributors"] [deps] diff --git a/docs/src/_changelog.md b/docs/src/_changelog.md index e86c274b..13d02187 100644 --- a/docs/src/_changelog.md +++ b/docs/src/_changelog.md @@ -7,12 +7,20 @@ CurrentModule = CurveFit This documents notable changes in CurveFit.jl. The format is based on [Keep a Changelog](https://keepachangelog.com). -## Unreleased +## [v1.9.3] ### Fixed -- Fixed parameter handling of [`ExpCurveFitAlgorithm`](@ref), - [`PowerCurveFitAlgorithm`](@ref), and [`KingCurveFitAlgorithm`](@ref) so they - don't return transformed parameters ([#112]). +- [`LinearCurveFitAlgorithm`](@ref) will now automatically invert the intercept + when `yfun` is given to ensure that the returned parameters match the values + being fitted. Also affects [`ExpCurveFitAlgorithm`](@ref) and + [`PowerCurveFitAlgorithm`](@ref). *This is considered a bugfix rather than a + breaking change.* +- Fixed the parameter handling and Jacobian of [`KingCurveFitAlgorithm`](@ref) + ([#112]). +- Fixed the statistics methods for [`LinearCurveFitAlgorithm`](@ref) when a + transform is applied with `yfun` ([#112]). Previously the Jacobian for a + linear function would be computed on the residuals stored in the original + y-space. ## [v1.9.2] - 2026-06-24 diff --git a/src/stats.jl b/src/stats.jl index fd0e222d..5c7cd4e4 100644 --- a/src/stats.jl +++ b/src/stats.jl @@ -102,11 +102,24 @@ end function jacobian(sol::CurveFitSolution{<:LinearCurveFitAlgorithm}) x = sol.prob.x - xfun = sol.alg.xfun - J = Matrix{eltype(x)}(undef, length(x), 2) - J[:, 1] .= xfun.(x) # Slope - J[:, 2] .= 1 # Intercept - return J + alg = sol.alg + if alg.yfun === identity + # Linear model, compute Jacobian analytically + J = Matrix{eltype(x)}(undef, length(x), 2) + J[:, 1] .= alg.xfun.(x) + J[:, 2] .= 1 + return J + else + # When there's a y-transform applied (e.g. power/exp fits) the residuals + # are in the original y-space so it would be incorrect to use the linear + # analytical Jacobian, instead we use DI. + a, b_stored = sol.u + u = [a, b_stored] + f_pred = u_curr -> begin + alg.yfun_inverse.(alg.yfun(u_curr[2]) .+ u_curr[1] .* alg.xfun.(x)) + end + return DifferentiationInterface.jacobian(f_pred, AutoForwardDiff(), u) + end end function jacobian(sol::CurveFitSolution{<:PolynomialFitAlgorithm}) diff --git a/test/stats.jl b/test/stats.jl index 207480c6..e7df97b9 100644 --- a/test/stats.jl +++ b/test/stats.jl @@ -276,4 +276,30 @@ using LinearSolve @test size(vcov(sol_rat_nl)) == (3, 3) @test all(stderror(sol_rat_nl) .> 0) end + + @testset "Jacobian in original y-space for transformed fits" begin + x = collect(1.0:10.0) + + # Power fit: y = b * x^a => ∂y/∂a = y*log(x), ∂y/∂b = y/b + a_true, b_true = 1.5, 3.0 + y_pow = @. b_true * x^a_true + sol_pow = solve(CurveFitProblem(x, y_pow), PowerCurveFitAlgorithm()) + a, b = sol_pow.u + @test [a, b] ≈ [a_true, b_true] + y_pred = sol_pow.(x) + J_pow = CurveFit.jacobian(sol_pow) + @test J_pow[:, 1] ≈ y_pred .* log.(x) + @test J_pow[:, 2] ≈ y_pred ./ b + + # Exp fit: y = b * exp(a*x) => ∂y/∂a = y*x, ∂y/∂b = y/b + a_true, b_true = 0.3, 2.0 + y_exp = @. b_true * exp(a_true * x) + sol_exp = solve(CurveFitProblem(x, y_exp), ExpCurveFitAlgorithm()) + a, b = sol_exp.u + @test [a, b] ≈ [a_true, b_true] + y_pred = sol_exp.(x) + J_exp = CurveFit.jacobian(sol_exp) + @test J_exp[:, 1] ≈ y_pred .* x + @test J_exp[:, 2] ≈ y_pred ./ b + end end