Skip to content

Commit

Permalink
store_trace without show_trace; maxTime (#229)
Browse files Browse the repository at this point in the history
  • Loading branch information
@viktmar
viktmar authored Oct 12, 2022
1 parent e9b9e87 commit 863bbca
Show file tree
Hide file tree
Showing 2 changed files with 61 additions and 53 deletions.
56 changes: 28 additions & 28 deletions src/curve_fit.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
struct LsqFitResult{P,R,J,W<:AbstractArray,T}
struct LsqFitResult{P,R,J,W <: AbstractArray,T}
param::P
resid::R
jacobian::J
Expand Down Expand Up @@ -28,30 +28,30 @@ end
# provide a method for those who have their own Jacobian function
function lmfit(f, g, p0::AbstractArray, wt::AbstractArray; kwargs...)
r = f(p0)
R = OnceDifferentiable(f, g, p0, copy(r); inplace = false)
R = OnceDifferentiable(f, g, p0, copy(r); inplace=false)
lmfit(R, p0, wt; kwargs...)
end

#for inplace f and inplace g
# for inplace f and inplace g
function lmfit(f!, g!, p0::AbstractArray, wt::AbstractArray, r::AbstractArray; kwargs...)
R = OnceDifferentiable(f!, g!, p0, copy(r); inplace = true)
R = OnceDifferentiable(f!, g!, p0, copy(r); inplace=true)
lmfit(R, p0, wt; kwargs...)
end

#for inplace f only
# for inplace f only
function lmfit(
f,
p0::AbstractArray,
wt::AbstractArray,
r::AbstractArray;
autodiff = :finite,
autodiff=:finite,
kwargs...,
)
R = OnceDifferentiable(f, p0, copy(r); inplace = true, autodiff = autodiff)
R = OnceDifferentiable(f, p0, copy(r); inplace=true, autodiff=autodiff)
lmfit(R, p0, wt; kwargs...)
end

function lmfit(f, p0::AbstractArray, wt::AbstractArray; autodiff = :finite, kwargs...)
function lmfit(f, p0::AbstractArray, wt::AbstractArray; autodiff=:finite, kwargs...)
# this is a convenience function for the curve_fit() methods
# which assume f(p) is the cost functionj i.e. the residual of a
# model where
Expand All @@ -68,15 +68,15 @@ function lmfit(f, p0::AbstractArray, wt::AbstractArray; autodiff = :finite, kwar
# construct Jacobian function, which uses finite difference method
r = f(p0)
autodiff = autodiff == :forwarddiff ? :forward : autodiff
R = OnceDifferentiable(f, p0, copy(r); inplace = false, autodiff = autodiff)
R = OnceDifferentiable(f, p0, copy(r); inplace=false, autodiff=autodiff)
lmfit(R, p0, wt; kwargs...)
end

function lmfit(
R::OnceDifferentiable,
p0::AbstractArray,
wt::AbstractArray;
autodiff = :finite,
autodiff=:finite,
kwargs...,
)
results = levenberg_marquardt(R, p0; kwargs...)
Expand Down Expand Up @@ -125,7 +125,7 @@ function curve_fit(
xdata::AbstractArray,
ydata::AbstractArray,
p0::AbstractArray;
inplace = false,
inplace=false,
kwargs...,
)
check_data_health(xdata, ydata)
Expand All @@ -147,7 +147,7 @@ function curve_fit(
xdata::AbstractArray,
ydata::AbstractArray,
p0::AbstractArray;
inplace = false,
inplace=false,
kwargs...,
)
check_data_health(xdata, ydata)
Expand All @@ -171,7 +171,7 @@ function curve_fit(
ydata::AbstractArray,
wt::AbstractArray,
p0::AbstractArray;
inplace = false,
inplace=false,
kwargs...,
)
check_data_health(xdata, ydata)
Expand All @@ -195,7 +195,7 @@ function curve_fit(
ydata::AbstractArray,
wt::AbstractArray,
p0::AbstractArray;
inplace = false,
inplace=false,
kwargs...,
)
check_data_health(xdata, ydata)
Expand Down Expand Up @@ -271,7 +271,7 @@ function estimate_covar(fit::LsqFitResult)
return covar
end

function StatsBase.stderror(fit::LsqFitResult; rtol::Real = NaN, atol::Real = 0)
function StatsBase.stderror(fit::LsqFitResult; rtol::Real=NaN, atol::Real=0)
# computes standard error of estimates from
# fit : a LsqFitResult from a curve_fit()
# atol : absolute tolerance for approximate comparisson to 0.0 in negativity check
Expand All @@ -283,21 +283,21 @@ function StatsBase.stderror(fit::LsqFitResult; rtol::Real = NaN, atol::Real = 0)
if !isapprox(
vratio,
0.0,
atol = atol,
rtol = isnan(rtol) ? Base.rtoldefault(vratio, 0.0, 0) : rtol,
atol=atol,
rtol=isnan(rtol) ? Base.rtoldefault(vratio, 0.0, 0) : rtol,
) && vratio < 0.0
error("Covariance matrix is negative for atol=$atol and rtol=$rtol")
end
return sqrt.(abs.(vars))
end

function margin_error(fit::LsqFitResult, alpha = 0.05; rtol::Real = NaN, atol::Real = 0)
function margin_error(fit::LsqFitResult, alpha=0.05; rtol::Real=NaN, atol::Real=0)
# computes margin of error at alpha significance level from
# fit : a LsqFitResult from a curve_fit()
# alpha : significance level, e.g. alpha=0.05 for 95% confidence
# atol : absolute tolerance for approximate comparisson to 0.0 in negativity check
# rtol : relative tolerance for approximate comparisson to 0.0 in negativity check
std_errors = stderror(fit; rtol = rtol, atol = atol)
std_errors = stderror(fit; rtol=rtol, atol=atol)
dist = TDist(dof(fit))
critical_values = quantile(dist, 1 - alpha / 2)
# scale standard errors by quantile of the student-t distribution (critical values)
Expand All @@ -306,25 +306,25 @@ end

function confidence_interval(
fit::LsqFitResult,
alpha = 0.05;
rtol::Real = NaN,
atol::Real = 0,
alpha=0.05;
rtol::Real=NaN,
atol::Real=0,
)
# computes confidence intervals at alpha significance level from
# fit : a LsqFitResult from a curve_fit()
# alpha : significance level, e.g. alpha=0.05 for 95% confidence
# atol : absolute tolerance for approximate comparisson to 0.0 in negativity check
# rtol : relative tolerance for approximate comparisson to 0.0 in negativity check
std_errors = stderror(fit; rtol = rtol, atol = atol)
margin_of_errors = margin_error(fit, alpha; rtol = rtol, atol = atol)
std_errors = stderror(fit; rtol=rtol, atol=atol)
margin_of_errors = margin_error(fit, alpha; rtol=rtol, atol=atol)
confidence_intervals =
collect(zip(coef(fit) - margin_of_errors, coef(fit) + margin_of_errors))
end

@deprecate standard_errors(args...; kwargs...) stderror(args...; kwargs...)
@deprecate estimate_errors(
fit::LsqFitResult,
confidence = 0.95;
rtol::Real = NaN,
atol::Real = 0,
) margin_error(fit, 1 - confidence; rtol = rtol, atol = atol)
confidence=0.95;
rtol::Real=NaN,
atol::Real=0,
) margin_error(fit, 1 - confidence; rtol=rtol, atol=atol)
58 changes: 33 additions & 25 deletions src/levenberg_marquardt.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,19 +79,21 @@ Comp & Applied Math).
function levenberg_marquardt(
df::OnceDifferentiable,
initial_x::AbstractVector{T};
x_tol::Real = 1e-8,
g_tol::Real = 1e-12,
maxIter::Integer = 1000,
lambda = T(10),
tau = T(Inf),
lambda_increase::Real = 10.0,
lambda_decrease::Real = 0.1,
min_step_quality::Real = 1e-3,
good_step_quality::Real = 0.75,
show_trace::Bool = false,
lower::AbstractVector{T} = Array{T}(undef, 0),
upper::AbstractVector{T} = Array{T}(undef, 0),
avv!::Union{Function,Nothing,Avv} = nothing,
x_tol::Real=1e-8,
g_tol::Real=1e-12,
maxIter::Integer=1000,
maxTime::Float64=Inf,
lambda=T(10),
tau=T(Inf),
lambda_increase::Real=10.0,
lambda_decrease::Real=0.1,
min_step_quality::Real=1e-3,
good_step_quality::Real=0.75,
show_trace::Bool=false,
store_trace::Bool=false,
lower::AbstractVector{T}=Array{T}(undef, 0),
upper::AbstractVector{T}=Array{T}(undef, 0),
avv!::Union{Function,Nothing,Avv}=nothing,
) where {T}

# First evaluation
Expand Down Expand Up @@ -155,14 +157,18 @@ function levenberg_marquardt(

# Maintain a trace of the system.
tr = LMTrace{LevenbergMarquardt}()
if show_trace
if show_trace || store_trace
d = Dict("lambda" => lambda)
os = LMState{LevenbergMarquardt}(iterCt, sum(abs2, value(df)), NaN, d)
push!(tr, os)
println(os)
if show_trace
println(os)
end
end

while (~converged && iterCt < maxIter)
startTime = time()

while (~converged && iterCt < maxIter && maxTime > time() - startTime)
# jacobian! will check if x is new or not, so it is only actually
# evaluated if x was updated last iteration.
jacobian!(df, x) # has alias J
Expand All @@ -175,7 +181,7 @@ function levenberg_marquardt(
# It is additionally useful to bound the elements of DtD below to help
# prevent "parameter evaporation".

DtD = vec(sum(abs2, J, dims = 1))
DtD = vec(sum(abs2, J, dims=1))
for i = 1:length(DtD)
if DtD[i] <= MIN_DIAGONAL
DtD[i] = MIN_DIAGONAL
Expand All @@ -187,23 +193,23 @@ function levenberg_marquardt(
@simd for i = 1:n
@inbounds JJ[i, i] += lambda * DtD[i]
end
#n_buffer is delta C, JJ is g compared to Mark's code
# n_buffer is delta C, JJ is g compared to Mark's code
mul!(n_buffer, transpose(J), value(df))
rmul!(n_buffer, -1)

v .= JJ \ n_buffer


if avv! != nothing
#GEODESIC ACCELERATION PART
# GEODESIC ACCELERATION PART
avv!(dir_deriv, x, v)
mul!(a, transpose(J), dir_deriv)
rmul!(a, -1) #we multiply by -1 before the decomposition/division
LAPACK.potrf!('U', JJ) #in place cholesky decomposition
LAPACK.potrs!('U', JJ, a) #divides a by JJ, taking into account the fact that JJ is now the `U` cholesky decoposition of what it was before
rmul!(a, -1) # we multiply by -1 before the decomposition/division
LAPACK.potrf!('U', JJ) # in place cholesky decomposition
LAPACK.potrs!('U', JJ, a) # divides a by JJ, taking into account the fact that JJ is now the `U` cholesky decoposition of what it was before
rmul!(a, 0.5)
delta_x .= v .+ a
#end of the GEODESIC ACCELERATION PART
# end of the GEODESIC ACCELERATION PART
else
delta_x = v
end
Expand Down Expand Up @@ -263,12 +269,14 @@ function levenberg_marquardt(
iterCt += 1

# show state
if show_trace
if show_trace || store_trace
g_norm = norm(J' * value(df), Inf)
d = Dict("g(x)" => g_norm, "dx" => copy(delta_x), "lambda" => lambda)
os = LMState{LevenbergMarquardt}(iterCt, sum(abs2, value(df)), g_norm, d)
push!(tr, os)
println(os)
if show_trace
println(os)
end
end

# check convergence criteria:
Expand Down

0 comments on commit 863bbca

Please sign in to comment.