From 55d8df58d7154e432df5d3fabb3fa45bef6691d8 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 13:24:25 -0600 Subject: [PATCH 01/27] refactor fit!(::LinearMixedModel) to support a more flexible optimization backend --- src/MixedModels.jl | 1 + src/linearmixedmodel.jl | 72 +------------------------ src/mixedmodel.jl | 27 ++++++++-- src/nlopt.jl | 98 +++++++++++++++++++++++++++++++++ src/optsummary.jl | 116 ++++++++++++++++++---------------------- src/serialization.jl | 7 +++ 6 files changed, 184 insertions(+), 137 deletions(-) create mode 100644 src/nlopt.jl diff --git a/src/MixedModels.jl b/src/MixedModels.jl index f3a9f7d4b..266f9e925 100644 --- a/src/MixedModels.jl +++ b/src/MixedModels.jl @@ -205,6 +205,7 @@ include("grouping.jl") include("mimeshow.jl") include("serialization.jl") include("profile/profile.jl") +include("nlopt.jl") include("prima.jl") # COV_EXCL_START diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 14f297835..d205aa753 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -177,7 +177,6 @@ function LinearMixedModel( θ = foldl(vcat, getθ(c) for c in reterms) optsum = OptSummary(θ, lbd) optsum.sigma = isnothing(σ) ? nothing : T(σ) - fill!(optsum.xtol_abs, 1.0e-10) return LinearMixedModel( form, reterms, @@ -459,55 +458,10 @@ function StatsAPI.fit!( ArgumentError("The response is constant and thus model fitting has failed") ) end - opt = Opt(optsum) optsum.REML = REML optsum.sigma = σ - prog = ProgressUnknown(; desc="Minimizing", showspeed=true) - # start from zero for the initial call to obj before optimization - iter = 0 + xmin, fmin = optimize!(m; progress, thin) fitlog = optsum.fitlog - function obj(x, g) - isempty(g) || throw(ArgumentError("g should be empty for this objective")) - iter += 1 - val = if isone(iter) && x == optsum.initial - optsum.finitial - else - try - objective(updateL!(setθ!(m, x))) - catch ex - # This can happen when the optimizer drifts into an area where - # there isn't enough shrinkage. Why finitial? Generally, it will - # be the (near) worst case scenario value, so the optimizer won't - # view it as an optimum. Using Inf messes up the quadratic - # approximation in BOBYQA. - ex isa PosDefException || rethrow() - optsum.finitial - end - end - progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - !isone(iter) && iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) - return val - end - NLopt.min_objective!(opt, obj) - try - # use explicit evaluation w/o calling opt to avoid confusing iteration count - optsum.finitial = objective(updateL!(setθ!(m, optsum.initial))) - catch ex - ex isa PosDefException || rethrow() - # give it one more try with a massive change in scaling - @info "Initial objective evaluation failed, rescaling initial guess and trying again." - @warn """Failure of the initial evaluation is often indicative of a model specification - that is not well supported by the data and/or a poorly scaled model. - """ - optsum.initial ./= - (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * - maximum(response(m)) - optsum.finitial = objective(updateL!(setθ!(m, optsum.initial))) - end - empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) - fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial)) - ProgressMeter.finish!(prog) ## check if small non-negative parameter values can be set to zero xmin_ = copy(xmin) lb = optsum.lowerbd @@ -518,7 +472,7 @@ function StatsAPI.fit!( end loglength = length(fitlog) if xmin ≠ xmin_ - if (zeroobj = obj(xmin_, T[])) ≤ (fmin + optsum.ftol_zero_abs) + if (zeroobj = objective!(m, xmin_)) ≤ (fmin + optsum.ftol_zero_abs) fmin = zeroobj copyto!(xmin, xmin_) elseif length(fitlog) > loglength @@ -529,11 +483,8 @@ function StatsAPI.fit!( ## ensure that the parameter values saved in m are xmin updateL!(setθ!(m, xmin)) - optsum.feval = opt.numevals optsum.final = xmin optsum.fmin = fmin - optsum.returnvalue = ret - _check_nlopt_return(ret) return m end @@ -828,25 +779,6 @@ function objective(m::LinearMixedModel{T}) where {T} return isempty(wts) ? val : val - T(2.0) * sum(log, wts) end -""" - objective!(m::LinearMixedModel, θ) - objective!(m::LinearMixedModel) - -Equivalent to `objective(updateL!(setθ!(m, θ)))`. - -When `m` has a single, scalar random-effects term, `θ` can be a scalar. - -The one-argument method curries and returns a single-argument function of `θ`. - -Note that these methods modify `m`. -The calling function is responsible for restoring the optimal `θ`. -""" -function objective! end - -function objective!(m::LinearMixedModel) - return Base.Fix1(objective!, m) -end - function objective!(m::LinearMixedModel{T}, θ) where {T} return objective(updateL!(setθ!(m, θ))) end diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl index 36e11fc81..16618be7c 100644 --- a/src/mixedmodel.jl +++ b/src/mixedmodel.jl @@ -117,6 +117,25 @@ StatsAPI.modelmatrix(m::MixedModel) = m.feterm.x StatsAPI.nobs(m::MixedModel) = length(m.y) +""" + objective!(m::MixedModel, θ) + objective!(m::MixedModel) + +Equivalent to `objective(updateL!(setθ!(m, θ)))`. + +When `m` has a single, scalar random-effects term, `θ` can be a scalar. + +The one-argument method curries and returns a single-argument function of `θ`. + +Note that these methods modify `m`. +The calling function is responsible for restoring the optimal `θ`. +""" +function objective! end + +function objective!(m::MixedModel) + return Base.Fix1(objective!, m) +end + StatsAPI.predict(m::MixedModel) = fitted(m) function retbl(mat, trm) @@ -131,13 +150,13 @@ StatsAPI.adjr2(m::MixedModel) = r2(m) function StatsAPI.r2(m::MixedModel) @error ( """There is no uniquely defined coefficient of determination for mixed models - that has all the properties of the corresponding value for classical + that has all the properties of the corresponding value for classical linear models. The GLMM FAQ provides more detail: - + https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#how-do-i-compute-a-coefficient-of-determination-r2-or-an-analogue-for-glmms - Alternatively, MixedModelsExtras provides a naive implementation, but + Alternatively, MixedModelsExtras provides a naive implementation, but the warnings there and in the FAQ should be taken seriously! """ ) @@ -148,7 +167,7 @@ end raneftables(m::MixedModel; uscale = false) Return the conditional means of the random effects as a `NamedTuple` of Tables.jl-compliant tables. - + !!! note The API guarantee is only that the NamedTuple contains Tables.jl tables and not on the particular concrete type of each table. """ diff --git a/src/nlopt.jl b/src/nlopt.jl new file mode 100644 index 000000000..43d258ac9 --- /dev/null +++ b/src/nlopt.jl @@ -0,0 +1,98 @@ +push!(OPTIMIZATION_BACKENDS, :nlopt) + +const NLoptBackend = Val{:nlopt} + +function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thin::Int=tyepmax(Int)) + optsum = m.optsum + opt = Opt(optsum) + + function obj(x, g) + isempty(g) || throw(ArgumentError("g should be empty for this objective")) + iter += 1 + val = if isone(iter) && x == optsum.initial + optsum.finitial + else + try + objective!(m, x) + catch ex + # This can happen when the optimizer drifts into an area where + # there isn't enough shrinkage. Why finitial? Generally, it will + # be the (near) worst case scenario value, so the optimizer won't + # view it as an optimum. Using Inf messes up the quadratic + # approximation in BOBYQA. + ex isa PosDefException || rethrow() + optsum.finitial + end + end + progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) + !isone(iter) && iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + return val + end + NLopt.min_objective!(opt, obj) + prog = ProgressUnknown(; desc="Minimizing", showspeed=true) + # start from zero for the initial call to obj before optimization + iter = 0 + fitlog = optsum.fitlog + + try + # use explicit evaluation w/o calling opt to avoid confusing iteration count + optsum.finitial = objective!(m, optsum.initial) + catch ex + ex isa PosDefException || rethrow() + # give it one more try with a massive change in scaling + @info "Initial objective evaluation failed, rescaling initial guess and trying again." + @warn """Failure of the initial evaluation is often indicative of a model specification + that is not well supported by the data and/or a poorly scaled model. + """ + optsum.initial ./= + (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * + maximum(response(m)) + optsum.finitial = objective!(m, optsum.initial) + end + empty!(fitlog) + push!(fitlog, (copy(optsum.initial), optsum.finitial)) + fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial)) + ProgressMeter.finish!(prog) + optsum.feval = opt.numevals + optsum.returnvalue = ret + _check_nlopt_return(ret) + return xmin, fmin +end + +function NLopt.Opt(optsum::OptSummary) + lb = optsum.lowerbd + + opt = NLopt.Opt(optsum.optimizer, length(lb)) + NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective + NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective + NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values + if length(optsum.xtol_abs) == length(lb) # not true for fast=false optimization in GLMM + NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values + end + NLopt.lower_bounds!(opt, lb) + NLopt.maxeval!(opt, optsum.maxfeval) + NLopt.maxtime!(opt, optsum.maxtime) + if isempty(optsum.initial_step) + optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb)) + else + NLopt.initial_step!(opt, optsum.initial_step) + end + return opt +end + + +const _NLOPT_FAILURE_MODES = [ + :FAILURE, + :INVALID_ARGS, + :OUT_OF_MEMORY, + :FORCED_STOP, + :MAXEVAL_REACHED, + :MAXTIME_REACHED, +] + +function _check_nlopt_return(ret, failure_modes=_NLOPT_FAILURE_MODES) + ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited") + if ret ∈ failure_modes + @warn("NLopt optimization failure: $ret") + end +end diff --git a/src/optsummary.jl b/src/optsummary.jl index 864f09cd6..574cdd830 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -1,61 +1,83 @@ """ OptSummary -Summary of an `NLopt` optimization +Summary of an optimization # Fields + +## Tolerances, initial and final values * `initial`: a copy of the initial parameter values in the optimization * `finitial`: the initial value of the objective * `lowerbd`: lower bounds on the parameter values -* `ftol_rel`: as in NLopt -* `ftol_abs`: as in NLopt -* `xtol_rel`: as in NLopt -* `xtol_abs`: as in NLopt -* `initial_step`: as in NLopt -* `maxfeval`: as in NLopt (`maxeval`) -* `maxtime`: as in NLopt * `final`: a copy of the final parameter values from the optimization * `fmin`: the final value of the objective * `feval`: the number of function evaluations -* `optimizer`: the name of the optimizer used, as a `Symbol` -* `returnvalue`: the return value, as a `Symbol` + Available backends can be examined via [`OPTIMIZATION_BACKENDS`](@ref). +* `returnvalue`: the return value, as a `Symbol`. The available return values will differ between backends. * `xtol_zero_abs`: the tolerance for a near zero parameter to be considered practically zero * `ftol_zero_abs`: the tolerance for change in the objective for setting a near zero parameter to zero +* `maxfeval`: as in NLopt (`maxeval`) and PRIMA (`maxfun`) + +## Choice of optimizer and backend +* `optimizer`: the name of the optimizer used, as a `Symbol` +* `backend`: the optimization library providing the optimizer, default is `NLoptBackend`. + +## Backend-specific fields +* `ftol_rel`: as in NLopt, not used in PRIMA +* `ftol_abs`: as in NLopt, not used in PRIMA +* `xtol_rel`: as in NLopt, not used in PRIMA +* `xtol_abs`: as in NLopt, not used in PRIMA +* `initial_step`: as in NLopt, not used in PRIMA +* `maxtime`: as in NLopt, not used in PRIMA +* `rhobeg`: as in PRIMA, not used in NLopt +* `rhoend`: as in PRIMA, not used in NLopt + +## MixedModels-spcific fields, unrelated to the optimizer * `fitlog`: A vector of tuples of parameter and objectives values from steps in the optimization * `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs * `REML`: use the REML criterion for LMM fits * `sigma`: a priori value for the residual standard deviation for LMM -The last three fields are MixedModels functionality and not related directly to the `NLopt` package or algorithms. - !!! note The internal storage of the parameter values within `fitlog` may change in the future to use a different subtype of `AbstractVector` (e.g., `StaticArrays.SVector`) for each snapshot without being considered a breaking change. + +!!! note + The exact order and number of fields may change as support for additional backends and features + thereof are added. In other words: use the keyword constructor and do **not** use the positional + constructor. """ Base.@kwdef mutable struct OptSummary{T<:AbstractFloat} initial::Vector{T} + finitial::T = Inf * one(eltype(initial)) lowerbd::Vector{T} + final::Vector{T} = copy(initial) + fmin::T = Inf * one(eltype(initial)) + feval::Int = -1 + returnvalue::Symbol = :FAILURE + xtol_zero_abs::T = eltype(initial)(0.001) + ftol_zero_abs::T = eltype(initial)(1.e-5) + maxfeval::Int = -1 + + + optimizer::Symbol = :LN_BOBYQA + backend::Symbol = :nlopt + # the @kwdef macro isn't quite smart enough for us to use the type parameter # for the default values, but we can fake it - finitial::T = Inf * one(eltype(initial)) ftol_rel::T = eltype(initial)(1.0e-12) ftol_abs::T = eltype(initial)(1.0e-8) xtol_rel::T = zero(eltype(initial)) xtol_abs::Vector{T} = zero(initial) .+ 1e-10 initial_step::Vector{T} = empty(initial) - maxfeval::Int = -1 maxtime::T = -one(eltype(initial)) - feval::Int = -1 - final::Vector{T} = copy(initial) - fmin::T = Inf * one(eltype(initial)) - optimizer::Symbol = :LN_BOBYQA - returnvalue::Symbol = :FAILURE - xtol_zero_abs::T = eltype(initial)(0.001) - ftol_zero_abs::T = eltype(initial)(1.e-5) + + rhobeg::T = one(T) + rhoend::T = rhobeg / 1_000_000 + # not SVector because we would need to parameterize on size (which breaks GLMM) fitlog::Vector{Tuple{Vector{T},T}} = [(initial, fmin)] - # don't really belong here but I needed a place to store them nAGQ::Int = 1 REML::Bool = false sigma::Union{T,Nothing} = nothing @@ -121,50 +143,18 @@ function Base.show(io::IO, ::MIME"text/plain", s::OptSummary) return println(io, "Return code: ", s.returnvalue) end -Base.show(io::IO, s::OptSummary) = Base.show(io, MIME"text/plain"(), s) - -function NLopt.Opt(optsum::OptSummary) - lb = optsum.lowerbd - - opt = NLopt.Opt(optsum.optimizer, length(lb)) - NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective - NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective - NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameter values - if length(optsum.xtol_abs) == length(lb) # not true for fast=false optimization in GLMM - NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameter values - end - NLopt.lower_bounds!(opt, lb) - NLopt.maxeval!(opt, optsum.maxfeval) - NLopt.maxtime!(opt, optsum.maxtime) - if isempty(optsum.initial_step) - optsum.initial_step = NLopt.initial_step(opt, optsum.initial, similar(lb)) - else - NLopt.initial_step!(opt, optsum.initial_step) - end - return opt -end - -StructTypes.StructType(::Type{<:OptSummary}) = StructTypes.Mutable() -StructTypes.excludes(::Type{<:OptSummary}) = (:lowerbd,) - -const _NLOPT_FAILURE_MODES = [ - :FAILURE, - :INVALID_ARGS, - :OUT_OF_MEMORY, - :FORCED_STOP, - :MAXEVAL_REACHED, - :MAXTIME_REACHED, -] - -function _check_nlopt_return(ret, failure_modes=_NLOPT_FAILURE_MODES) - ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited") - if ret ∈ failure_modes - @warn("NLopt optimization failure: $ret") - end -end +Base.show(io::IO, s::OptSummary) = Base.show(io, MIME("text/plain"), s) function Base.:(==)(o1::OptSummary{T}, o2::OptSummary{T}) where {T} return all(fieldnames(OptSummary)) do fn return getfield(o1, fn) == getfield(o2, fn) end end + +""" + OPTIMIZATION_BACKENDS + +A list of currently available optimization backends. +""" +const OPTIMIZATION_BACKENDS = Symbol[] +optimize!(m::MixedModel; kwargs...) = optimize!(m, Val(m.optsum.backend); kwargs...) diff --git a/src/serialization.jl b/src/serialization.jl index 708504c18..e0e44d9f0 100644 --- a/src/serialization.jl +++ b/src/serialization.jl @@ -79,6 +79,9 @@ function restoreoptsum!(ops::OptSummary{T}, dict::AbstractDict) where {T} :ftol_zero_abs, # added in v4.25.0 :sigma, # added in v4.1.0 :fitlog, # added in v4.1.0 + :backend, # added in v4.30.0 + :rhobeg, # added in v4.30.0 + :rhoend, # added in v4.30.0 ) nmdiff = setdiff( propertynames(ops), # names in freshly created optsum @@ -120,6 +123,10 @@ function restoreoptsum!(ops::OptSummary{T}, dict::AbstractDict) where {T} return ops end + +StructTypes.StructType(::Type{<:OptSummary}) = StructTypes.Mutable() +StructTypes.excludes(::Type{<:OptSummary}) = (:lowerbd,) + """ saveoptsum(io::IO, m::MixedModel) saveoptsum(filename, m::MixedModel) From 55f291c3b6e9f4f9be9682a44635e798826ad6bb Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 19:48:44 +0000 Subject: [PATCH 02/27] Update src/optsummary.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/optsummary.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/optsummary.jl b/src/optsummary.jl index 574cdd830..96417da92 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -32,7 +32,7 @@ Summary of an optimization * `rhobeg`: as in PRIMA, not used in NLopt * `rhoend`: as in PRIMA, not used in NLopt -## MixedModels-spcific fields, unrelated to the optimizer +## MixedModels-specific fields, unrelated to the optimizer * `fitlog`: A vector of tuples of parameter and objectives values from steps in the optimization * `nAGQ`: number of adaptive Gauss-Hermite quadrature points in deviance evaluation for GLMMs * `REML`: use the REML criterion for LMM fits From eba1859b0d557acf63b3f55da5c54dba5977bf79 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 14:49:40 -0600 Subject: [PATCH 03/27] fix crossref --- src/optsummary.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/optsummary.jl b/src/optsummary.jl index 96417da92..5eae6375e 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -12,7 +12,7 @@ Summary of an optimization * `final`: a copy of the final parameter values from the optimization * `fmin`: the final value of the objective * `feval`: the number of function evaluations - Available backends can be examined via [`OPTIMIZATION_BACKENDS`](@ref). + Available backends can be examined via `OPTIMIZATION_BACKENDS`. * `returnvalue`: the return value, as a `Symbol`. The available return values will differ between backends. * `xtol_zero_abs`: the tolerance for a near zero parameter to be considered practically zero * `ftol_zero_abs`: the tolerance for change in the objective for setting a near zero parameter to zero From 179dd7d4ba5d546684cbe4fe016723b8b27611cf Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 14:52:23 -0600 Subject: [PATCH 04/27] NEWS and version bump --- NEWS.md | 5 +++++ Project.toml | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index c34ee21d1..2b0e3dc7e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +MixedModels v4.30.0 Release Notes +============================== +- Refactor calls to backend optimizer to make it easier to add and use different optimization backends. [#802] + MixedModels v4.29.1 Release Notes ============================== - Populate `optsum` in `prfit!` call. [#801] @@ -595,3 +599,4 @@ Package dependencies [#795]: https://github.com/JuliaStats/MixedModels.jl/issues/795 [#799]: https://github.com/JuliaStats/MixedModels.jl/issues/799 [#801]: https://github.com/JuliaStats/MixedModels.jl/issues/801 +[#802]: https://github.com/JuliaStats/MixedModels.jl/issues/802 diff --git a/Project.toml b/Project.toml index 6fc061d03..ef14b3f15 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MixedModels" uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" author = ["Phillip Alday ", "Douglas Bates "] -version = "4.29.1" +version = "4.30.0" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" From 06bdf1af48240bcbd0520a555174143a7bb6ba23 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 15:34:21 -0600 Subject: [PATCH 05/27] OptSummary show specialization --- src/nlopt.jl | 2 ++ src/optsummary.jl | 51 +++++++++++++++++++++++++++++++++++++---------- 2 files changed, 43 insertions(+), 10 deletions(-) diff --git a/src/nlopt.jl b/src/nlopt.jl index 43d258ac9..8b74ec7d1 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -96,3 +96,5 @@ function _check_nlopt_return(ret, failure_modes=_NLOPT_FAILURE_MODES) @warn("NLopt optimization failure: $ret") end end + +opt_params(::NLoptBackend) = (:ftol_rel, :ftol_abs, :xtol_rel, :xtol_abs, :initial_step, :maxfeval, :maxtime) diff --git a/src/optsummary.jl b/src/optsummary.jl index 5eae6375e..f693d140f 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -60,7 +60,6 @@ Base.@kwdef mutable struct OptSummary{T<:AbstractFloat} ftol_zero_abs::T = eltype(initial)(1.e-5) maxfeval::Int = -1 - optimizer::Symbol = :LN_BOBYQA backend::Symbol = :nlopt @@ -127,20 +126,22 @@ function Base.show(io::IO, ::MIME"text/plain", s::OptSummary) println(io, "Initial parameter vector: ", s.initial) println(io, "Initial objective value: ", s.finitial) println(io) - println(io, "Optimizer (from NLopt): ", s.optimizer) + println(io, "Backend: ", s.backend) + println(io, "Optimizer: ", s.optimizer) println(io, "Lower bounds: ", s.lowerbd) - println(io, "ftol_rel: ", s.ftol_rel) - println(io, "ftol_abs: ", s.ftol_abs) - println(io, "xtol_rel: ", s.xtol_rel) - println(io, "xtol_abs: ", s.xtol_abs) - println(io, "initial_step: ", s.initial_step) - println(io, "maxfeval: ", s.maxfeval) - println(io, "maxtime: ", s.maxtime) + + for param in opt_params(Val(s.backend)) + println(io, rpad(string(param, ":"), length("Initial parameter vector: ")), + getfield(s, param)) + end println(io) println(io, "Function evaluations: ", s.feval) + println(io, "xtol_zero_abs: ", s.xtol_zero_abs) + println(io, "ftol_zero_abs: ", s.ftol_zero_abs) println(io, "Final parameter vector: ", s.final) println(io, "Final objective value: ", s.fmin) - return println(io, "Return code: ", s.returnvalue) + println(io, "Return code: ", s.returnvalue) + return nothing end Base.show(io::IO, s::OptSummary) = Base.show(io, MIME("text/plain"), s) @@ -158,3 +159,33 @@ A list of currently available optimization backends. """ const OPTIMIZATION_BACKENDS = Symbol[] optimize!(m::MixedModel; kwargs...) = optimize!(m, Val(m.optsum.backend); kwargs...) + +""" + optimize!(::LinearMixedModel, ::Val{backend}; kwargs...) + optimize!(::GeneralizedLinearMixedModel, ::Val{backend}; kwargs...) + +Perform optimization on a mixed model, minimizing the objective. + +`optimize!` set ups the call to the backend optimizer using the options contained in the +model's `optsum` field. It then calls that optimizer and returns `xmin, fmin`. +Providing support for a new backend involves defining appropriate `optimize!` methods +with the second argument of type `::Val{:backend_name}` and adding `:backend_name` to +`OPTIMIZATION_BACKENDS`. Similarly, a method `opt_params(::Val{:backend_name})` should +be defined, which returns the optimization paramers (e.g. `xtol_abs` or `rho_end`) used +by the backend. + +Common keyword arguments are `progress` to show a progress meter and `thin` to control thinning +of the fitlog. +""" +function optimize! end + + +""" + opt_params(::Val{backend}) + +Return a collection of the fields of [`OptSummary`](@ref) used by backend. + +`:xtol_zero_abs`, `:ftol_zero_abs` do not need to be specified because +they are used _after_ optimization and are thus shared across backends. +""" +function opt_params end From 2c0580be2d2abceb925c7012c4af624382478c7b Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 16:08:55 -0600 Subject: [PATCH 06/27] update PRIMA backend to use new infrastructure --- ext/MixedModelsPRIMAExt.jl | 69 +++++++++++++++++++++++++++++--------- src/nlopt.jl | 2 +- src/prima.jl | 21 +----------- test/prima.jl | 2 ++ 4 files changed, 57 insertions(+), 37 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 6e80d8995..e08318a66 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -5,21 +5,32 @@ using MixedModels: ProgressMeter, ProgressUnknown using PRIMA: PRIMA function MixedModels.prfit!(m::LinearMixedModel; - progress::Bool=true, - REML::Bool=m.optsum.REML, - σ::Union{Real,Nothing}=m.optsum.sigma, - thin::Int=1) - optsum = m.optsum - copyto!(optsum.final, optsum.initial) - optsum.REML = REML - optsum.sigma = σ - optsum.finitial = objective!(m, optsum.initial) + kwargs...) + + unfit!(m) + m.optsum.optimizer = :bobyqa + m.optsum.backend = :prima + + return fit!(m; kwargs...) +end + +prima_optimizer!(::Val{:bobyqa}, args...; kwargs...) = PRIMA.bobyqa!(args...; kwargs...) +prima_optimizer!(::Val{:cobyla}, args...; kwargs...) = PRIMA.cobyla(args...; kwargs...) +prima_optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa(args...; kwargs...) +push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) + +const PRIMABackend = Val{:prima} + +function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bool=true, thin::Int=tyepmax(Int)) + optsum = m.optsum prog = ProgressUnknown(; desc="Minimizing", showspeed=true) # start from zero for the initial call to obj before optimization iter = 0 - fitlog = empty!(optsum.fitlog) + fitlog = optsum.fitlog + function obj(x) + isempty(g) || throw(ArgumentError("g should be empty for this objective")) iter += 1 val = if isone(iter) && x == optsum.initial optsum.finitial @@ -37,19 +48,45 @@ function MixedModels.prfit!(m::LinearMixedModel; end end progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - if isone(iter) || iszero(rem(iter, thin)) - push!(fitlog, (copy(x), val)) - end + !isone(iter) && iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) return val end + try + # use explicit evaluation w/o calling opt to avoid confusing iteration count + optsum.finitial = objective!(m, optsum.initial) + catch ex + ex isa PosDefException || rethrow() + # give it one more try with a massive change in scaling + @info "Initial objective evaluation failed, rescaling initial guess and trying again." + @warn """Failure of the initial evaluation is often indicative of a model specification + that is not well supported by the data and/or a poorly scaled model. + """ + optsum.initial ./= + (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * + maximum(response(m)) + optsum.finitial = objective!(m, optsum.initial) + end + empty!(fitlog) + push!(fitlog, (copy(optsum.initial), optsum.finitial)) + info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; xl=m.optsum.lowerbd) ProgressMeter.finish!(prog) - info = PRIMA.bobyqa!(obj, optsum.final; xl=m.optsum.lowerbd) optsum.feval = info.nf optsum.fmin = info.fx optsum.returnvalue = Symbol(info.status) - optsum.optimizer = :PRIMA_BOBYQA - return m + + return optsum.final, optsum.fmin end + +function _check_prima_return(info::PRIMA.Info) + if !PRIMA.issuccess(info) + @warn "PRIMA optimization failure: $(ret)\n$(PRIMA.reason(info))" + end + + return nothing +end + +MixedModels.opt_params(::PRIMABackend) = (:rhobeg, :rhoend) + end # module diff --git a/src/nlopt.jl b/src/nlopt.jl index 8b74ec7d1..eaa90b134 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -6,6 +6,7 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thi optsum = m.optsum opt = Opt(optsum) + prog = ProgressUnknown(; desc="Minimizing", showspeed=true) function obj(x, g) isempty(g) || throw(ArgumentError("g should be empty for this objective")) iter += 1 @@ -29,7 +30,6 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thi return val end NLopt.min_objective!(opt, obj) - prog = ProgressUnknown(; desc="Minimizing", showspeed=true) # start from zero for the initial call to obj before optimization iter = 0 fitlog = optsum.fitlog diff --git a/src/prima.jl b/src/prima.jl index 7dd89d11d..0d9cbe823 100644 --- a/src/prima.jl +++ b/src/prima.jl @@ -1,6 +1,6 @@ """ - prfit!(m::LinearMixedModel) + prfit!(m::LinearMixedModel; kwargs...) Fit a mixed model using the [PRIMA](https://github.com/libprima/PRIMA.jl) implementation of the BOBYQA optimizer. @@ -13,25 +13,6 @@ of the BOBYQA optimizer. and potentially breaking changes it would take to fully replace the current NLopt backend with a PRIMA backend or a backend supporting a range of optimizers. -!!! note "OptSummary" - As part of this experimental foray, the structure of [`OptSummary`](@ref) is - not changed. This means that some fields of `OptSummary` are inaccurate when - examining a model fit with PRIMA. The following fields are unused with PRIMA - fits: - - - all tolerances (`ftol_rel`, `ftol_abs`, `xtol_rel`, `xtol_abs`) - - optimization timeouts (`maxfeval`, `maxtime`) - - `initial_step` - - The following fields have a different meaning when used with PRIMA: - - - `returnvalue` is populated with a symbol representing the PRIMA - return value, which PRIMA represents as an enum. - - `optimizer` is populated with a dummy value, indicating that a model was - fit with PRIMA. If you wish to refit the model with the NLOpt backend, - you will need to update the field with an appropriate NLOpt optimizer, - e.g. `:LN_BOBYQA` - !!! note "Package extension" In order to reduce the dependency burden, all methods of this function are implemented in a package extension and are only defined when PRIMA.jl is loaded diff --git a/test/prima.jl b/test/prima.jl index ff32cf03d..c1c9b1718 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -10,4 +10,6 @@ include("modelcache.jl") prmodel = prfit!(LinearMixedModel(formula(model), dataset(:sleepstudy)); progress=false) @test isapprox(loglikelihood(model), loglikelihood(prmodel)) + @test prmodel.optsum.optimizer == :bobyqa + @test prmodel.optsum.backend == :prima end From 847334ea1e311a9863253f29523f00d535b9fc2e Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 22:11:15 +0000 Subject: [PATCH 07/27] Update src/optsummary.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/optsummary.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/optsummary.jl b/src/optsummary.jl index f693d140f..f5c1dd930 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -171,7 +171,7 @@ model's `optsum` field. It then calls that optimizer and returns `xmin, fmin`. Providing support for a new backend involves defining appropriate `optimize!` methods with the second argument of type `::Val{:backend_name}` and adding `:backend_name` to `OPTIMIZATION_BACKENDS`. Similarly, a method `opt_params(::Val{:backend_name})` should -be defined, which returns the optimization paramers (e.g. `xtol_abs` or `rho_end`) used +be defined, which returns the optimization parameters (e.g. `xtol_abs` or `rho_end`) used by the backend. Common keyword arguments are `progress` to show a progress meter and `thin` to control thinning From cbad1c45f1ff341baf0ddc46e62e0e90f9b50f72 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 16:19:24 -0600 Subject: [PATCH 08/27] cleanup --- NEWS.md | 3 ++- ext/MixedModelsPRIMAExt.jl | 12 +++++++----- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index 2b0e3dc7e..5d8c8d614 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,7 @@ MixedModels v4.30.0 Release Notes ============================== -- Refactor calls to backend optimizer to make it easier to add and use different optimization backends. [#802] +- Refactor calls to backend optimizer to make it easier to add and use different optimization backends. + The structure of `OptSummary` has been expanded and `prfit!` has been updated to use this new structure. [#802] MixedModels v4.29.1 Release Notes ============================== diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index e08318a66..16a3510f9 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -1,13 +1,13 @@ module MixedModelsPRIMAExt -using MixedModels: MixedModels, LinearMixedModel, objective! -using MixedModels: ProgressMeter, ProgressUnknown +using MixedModels +using MixedModels: ProgressMeter, ProgressUnknown, objective! using PRIMA: PRIMA function MixedModels.prfit!(m::LinearMixedModel; kwargs...) - unfit!(m) + MixedModels.unfit!(m) m.optsum.optimizer = :bobyqa m.optsum.backend = :prima @@ -69,7 +69,9 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo end empty!(fitlog) push!(fitlog, (copy(optsum.initial), optsum.finitial)) - info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; xl=m.optsum.lowerbd) + info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; + xl=optsum.lowerbd, maxfun=optsum.maxfeval, + optsum.rhoend, optsumrhobeg) ProgressMeter.finish!(prog) optsum.feval = info.nf optsum.fmin = info.fx @@ -87,6 +89,6 @@ function _check_prima_return(info::PRIMA.Info) return nothing end -MixedModels.opt_params(::PRIMABackend) = (:rhobeg, :rhoend) +MixedModels.opt_params(::PRIMABackend) = (:rhobeg, :rhoend, :maxfeval) end # module From 095ef336816c1c83279d6295f7e5ef5eb5a07f49 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 16:21:55 -0600 Subject: [PATCH 09/27] :facepalm: --- ext/MixedModelsPRIMAExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 16a3510f9..12648cfd7 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -71,7 +71,7 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo push!(fitlog, (copy(optsum.initial), optsum.finitial)) info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; xl=optsum.lowerbd, maxfun=optsum.maxfeval, - optsum.rhoend, optsumrhobeg) + optsum.rhoend, optsum.rhobeg) ProgressMeter.finish!(prog) optsum.feval = info.nf optsum.fmin = info.fx From e250d9c8ef2555f48bbd9cb83570a60bfc9f56e0 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 16:40:41 -0600 Subject: [PATCH 10/27] more cleanup --- ext/MixedModelsPRIMAExt.jl | 19 ++----------------- src/linearmixedmodel.jl | 17 +++++++++++++++++ src/nlopt.jl | 15 --------------- 3 files changed, 19 insertions(+), 32 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 12648cfd7..2ccd85f3d 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -30,7 +30,6 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo fitlog = optsum.fitlog function obj(x) - isempty(g) || throw(ArgumentError("g should be empty for this objective")) iter += 1 val = if isone(iter) && x == optsum.initial optsum.finitial @@ -52,25 +51,11 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo return val end - try - # use explicit evaluation w/o calling opt to avoid confusing iteration count - optsum.finitial = objective!(m, optsum.initial) - catch ex - ex isa PosDefException || rethrow() - # give it one more try with a massive change in scaling - @info "Initial objective evaluation failed, rescaling initial guess and trying again." - @warn """Failure of the initial evaluation is often indicative of a model specification - that is not well supported by the data and/or a poorly scaled model. - """ - optsum.initial ./= - (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * - maximum(response(m)) - optsum.finitial = objective!(m, optsum.initial) - end empty!(fitlog) push!(fitlog, (copy(optsum.initial), optsum.finitial)) + maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; - xl=optsum.lowerbd, maxfun=optsum.maxfeval, + xl=optsum.lowerbd, maxfun, optsum.rhoend, optsum.rhobeg) ProgressMeter.finish!(prog) optsum.feval = info.nf diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index d205aa753..cdd804296 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -460,6 +460,23 @@ function StatsAPI.fit!( end optsum.REML = REML optsum.sigma = σ + + try + # use explicit evaluation w/o calling opt to avoid confusing iteration count + optsum.finitial = objective!(m, optsum.initial) + catch ex + ex isa PosDefException || rethrow() + # give it one more try with a massive change in scaling + @info "Initial objective evaluation failed, rescaling initial guess and trying again." + @warn """Failure of the initial evaluation is often indicative of a model specification + that is not well supported by the data and/or a poorly scaled model. + """ + optsum.initial ./= + (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * + maximum(response(m)) + optsum.finitial = objective!(m, optsum.initial) + end + xmin, fmin = optimize!(m; progress, thin) fitlog = optsum.fitlog ## check if small non-negative parameter values can be set to zero diff --git a/src/nlopt.jl b/src/nlopt.jl index eaa90b134..8486aafec 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -34,21 +34,6 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thi iter = 0 fitlog = optsum.fitlog - try - # use explicit evaluation w/o calling opt to avoid confusing iteration count - optsum.finitial = objective!(m, optsum.initial) - catch ex - ex isa PosDefException || rethrow() - # give it one more try with a massive change in scaling - @info "Initial objective evaluation failed, rescaling initial guess and trying again." - @warn """Failure of the initial evaluation is often indicative of a model specification - that is not well supported by the data and/or a poorly scaled model. - """ - optsum.initial ./= - (isempty(m.sqrtwts) ? 1.0 : maximum(m.sqrtwts)^2) * - maximum(response(m)) - optsum.finitial = objective!(m, optsum.initial) - end empty!(fitlog) push!(fitlog, (copy(optsum.initial), optsum.finitial)) fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial)) From 3e3cd0949f52cbd19becb0ff4215ec40a34b4c0e Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sat, 25 Jan 2025 16:59:43 -0600 Subject: [PATCH 11/27] test other PRIMA optimizers --- ext/MixedModelsPRIMAExt.jl | 5 +++-- test/prima.jl | 11 +++++++++++ 2 files changed, 14 insertions(+), 2 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 2ccd85f3d..014f34678 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -15,8 +15,8 @@ function MixedModels.prfit!(m::LinearMixedModel; end prima_optimizer!(::Val{:bobyqa}, args...; kwargs...) = PRIMA.bobyqa!(args...; kwargs...) -prima_optimizer!(::Val{:cobyla}, args...; kwargs...) = PRIMA.cobyla(args...; kwargs...) -prima_optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa(args...; kwargs...) +prima_optimizer!(::Val{:cobyla}, args...; kwargs...) = PRIMA.cobyla!(args...; kwargs...) +prima_optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa!(args...; kwargs...) push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) @@ -61,6 +61,7 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo optsum.feval = info.nf optsum.fmin = info.fx optsum.returnvalue = Symbol(info.status) + _check_prima_return(info) return optsum.final, optsum.fmin end diff --git a/test/prima.jl b/test/prima.jl index c1c9b1718..18e622c00 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -13,3 +13,14 @@ include("modelcache.jl") @test prmodel.optsum.optimizer == :bobyqa @test prmodel.optsum.backend == :prima end + +model = first(models(:sleepstudy)) +prmodel = LinearMixedModel(formula(model), dataset(:sleepstudy)) +prmodel.optsum.backend = :prima + +@testset "$optimizer" for optimizer in (:cobyla, :lincoa) + MixedModels.unfit!(prmodel) + prmodel.optsum.optimizer = optimizer + fit!(prmodel) + @test isapprox(loglikelihood(model), loglikelihood(prmodel)) +end From 9052d4758b173edbe2bea4cc2ecf8ed7b97023a0 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 00:14:58 +0000 Subject: [PATCH 12/27] nlopt for GLMM --- src/generalizedlinearmixedmodel.jl | 56 ++++++++++++------------------ src/nlopt.jl | 47 ++++++++++++++++++++++++- 2 files changed, 69 insertions(+), 34 deletions(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index d6b998b36..f1b9c2190 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -277,35 +277,10 @@ function StatsAPI.fit!( optsum.initial = vcat(β, lm.optsum.final) optsum.final = copy(optsum.initial) end - setpar! = fast ? setθ! : setβθ! - prog = ProgressUnknown(; desc="Minimizing", showspeed=true) - # start from zero for the initial call to obj before optimization - iter = 0 + + xmin, fmin = optimize!(m; progress, thin, fast, verbose) fitlog = optsum.fitlog - function obj(x, g) - isempty(g) || throw(ArgumentError("g should be empty for this objective")) - val = try - deviance(pirls!(setpar!(m, x), fast, verbose), nAGQ) - catch ex - # this allows us to recover from models where e.g. the link isn't - # as constraining as it should be - ex isa Union{PosDefException,DomainError} || rethrow() - iter == 1 && rethrow() - m.optsum.finitial - end - iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) - verbose && println(round(val; digits=5), " ", x) - progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - iter += 1 - return val - end - opt = Opt(optsum) - NLopt.min_objective!(opt, obj) - optsum.finitial = obj(optsum.initial, T[]) - empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) - fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial)) - ProgressMeter.finish!(prog) + ## check if very small parameter values bounded below by zero can be set to zero xmin_ = copy(xmin) for i in eachindex(xmin_) @@ -324,13 +299,10 @@ function StatsAPI.fit!( end end ## ensure that the parameter values saved in m are xmin - pirls!(setpar!(m, xmin), fast, verbose) - optsum.nAGQ = nAGQ - optsum.feval = opt.numevals + objective!(m, xmin; fast, verbose, nAGQ) optsum.final = xmin optsum.fmin = fmin - optsum.returnvalue = ret - _check_nlopt_return(ret) + optsum.nAGQ = nAGQ return m end @@ -540,6 +512,24 @@ function StatsAPI.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T} return accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2 end +# Base.Fix1 doesn't forward kwargs +objective!(m::GeneralizedLinearMixedModel; fast=false, kwargs...) = x -> _objective!(m, x, Val(fast); kwargs...) + +function objective!(m::GeneralizedLinearMixedModel{T}, x; fast=false, kwargs...) where {T} + return _objective!(m, x, Val(fast); kwargs...) +end + +# normally, it doesn't make sense to move a simple branch to dispatch +# HOWEVER, this wins up getting called in optimization a lot and +# moving this to a realization here allows us to avoid dynamic dispatch on setθ! / setθβ! +function _objective!(m::GeneralizedLinearMixedModel{T}, x, ::Val{true}; nAGQ=1, verbose=false) where {T} + deviance(pirls!(setθ!(m, x), true, verbose), nAGQ) +end + +function _objective!(m::GeneralizedLinearMixedModel{T}, x, ::Val{false}; nAGQ=1, verbose=false) where {T} + deviance(pirls!(setβθ!(m, x), false, verbose), nAGQ) +end + function Base.propertynames(m::GeneralizedLinearMixedModel, private::Bool=false) return ( :A, diff --git a/src/nlopt.jl b/src/nlopt.jl index 8486aafec..754fa69f0 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -2,7 +2,9 @@ push!(OPTIMIZATION_BACKENDS, :nlopt) const NLoptBackend = Val{:nlopt} -function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thin::Int=tyepmax(Int)) +function optimize!(m::LinearMixedModel, ::NLoptBackend; + progress::Bool=true, thin::Int=tyepmax(Int), + kwargs...) optsum = m.optsum opt = Opt(optsum) @@ -44,6 +46,49 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; progress::Bool=true, thi return xmin, fmin end +function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; + progress::Bool=true, thin::Int=tyepmax(Int), + fast::Bool=false, verbose::Bool=false, nAGQ=1, + kwargs...) + optsum = m.optsum + fitlog = optsum.fitlog + prog = ProgressUnknown(; desc="Minimizing", showspeed=true) + # start from zero for the initial call to obj before optimization + iter = 0 + fitlog = optsum.fitlog + function obj(x, g) + isempty(g) || throw(ArgumentError("g should be empty for this objective")) + val = try + _objective!(m, x, Val(fast); verbose, nAGQ) + catch ex + # this allows us to recover from models where e.g. the link isn't + # as constraining as it should be + ex isa Union{PosDefException,DomainError} || rethrow() + iter == 1 && rethrow() + m.optsum.finitial + end + iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + verbose && println(round(val; digits=5), " ", x) + progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) + iter += 1 + return val + end + opt = Opt(optsum) + NLopt.min_objective!(opt, obj) + optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) + empty!(fitlog) + push!(fitlog, (copy(optsum.initial), optsum.finitial)) + fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial)) + ProgressMeter.finish!(prog) + + optsum.feval = opt.numevals + optsum.returnvalue = ret + _check_nlopt_return(ret) + + return xmin, fmin + +end + function NLopt.Opt(optsum::OptSummary) lb = optsum.lowerbd From 56af378c3cf79b08818c2e07f902b55cccb1ad0f Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 00:36:39 +0000 Subject: [PATCH 13/27] prima for GLMM --- ext/MixedModelsPRIMAExt.jl | 46 +++++++++++++++++++++++++++++- src/generalizedlinearmixedmodel.jl | 7 ++++- src/linearmixedmodel.jl | 29 +++++++++---------- src/nlopt.jl | 1 - 4 files changed, 65 insertions(+), 18 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 014f34678..6fa5f830e 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -22,7 +22,8 @@ push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) const PRIMABackend = Val{:prima} -function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bool=true, thin::Int=tyepmax(Int)) +function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; + progress::Bool=true, thin::Int=tyepmax(Int)) optsum = m.optsum prog = ProgressUnknown(; desc="Minimizing", showspeed=true) # start from zero for the initial call to obj before optimization @@ -66,6 +67,49 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bo return optsum.final, optsum.fmin end +function optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; + progress::Bool=true, thin::Int=tyepmax(Int), + fast::Bool=false, verbose::Bool=false, nAGQ=1, + kwargs...) + optsum = m.optsum + fitlog = optsum.fitlog + prog = ProgressUnknown(; desc="Minimizing", showspeed=true) + # start from zero for the initial call to obj before optimization + iter = 0 + fitlog = optsum.fitlog + function obj(x) + val = try + _objective!(m, x, Val(fast); verbose, nAGQ) + catch ex + # this allows us to recover from models where e.g. the link isn't + # as constraining as it should be + ex isa Union{PosDefException,DomainError} || rethrow() + iter == 1 && rethrow() + m.optsum.finitial + end + iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + verbose && println(round(val; digits=5), " ", x) + progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) + iter += 1 + return val + end + + optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) + empty!(fitlog) + push!(fitlog, (copy(optsum.initial), optsum.finitial)) + maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) + info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; + xl=optsum.lowerbd, maxfun, + optsum.rhoend, optsum.rhobeg) + ProgressMeter.finish!(prog) + + optsum.feval = info.nf + optsum.fmin = info.fx + optsum.returnvalue = Symbol(info.status) + _check_prima_return(info) + + return xmin, fmin +end function _check_prima_return(info::PRIMA.Info) if !PRIMA.issuccess(info) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index f1b9c2190..3958e430f 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -248,6 +248,8 @@ function StatsAPI.fit!( progress::Bool=true, thin::Int=typemax(Int), init_from_lmm=Set(), + backend::Symbol=m.optsum.backend, + optimizer::Symbol=m.optsum.optimizer, ) where {T} β = copy(m.β) θ = copy(m.θ) @@ -278,7 +280,10 @@ function StatsAPI.fit!( optsum.final = copy(optsum.initial) end - xmin, fmin = optimize!(m; progress, thin, fast, verbose) + optsum.backend = backend + optsum.optimizer = optimizer + + xmin, fmin = optimize!(m; progress, thin, fast, verbose, nAGQ) fitlog = optsum.fitlog ## check if very small parameter values bounded below by zero can be set to zero diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index cdd804296..2c37fb2aa 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -205,21 +205,16 @@ function StatsAPI.fit( ) end -function StatsAPI.fit( - ::Type{LinearMixedModel}, - f::FormulaTerm, - tbl::Tables.ColumnTable; - wts=[], - contrasts=Dict{Symbol,Any}(), - progress=true, - REML=false, - σ=nothing, - thin=typemax(Int), - amalgamate=true, -) - return fit!( - LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate); progress, REML, thin - ) +function StatsAPI.fit(::Type{LinearMixedModel}, + f::FormulaTerm, + tbl::Tables.ColumnTable; + wts=[], + contrasts=Dict{Symbol,Any}(), + σ=nothing, + amalgamate=true, + kwargs...) + lmm = LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate) + return fit!(lmm; kwargs...) end function _offseterr() @@ -447,6 +442,8 @@ function StatsAPI.fit!( REML::Bool=m.optsum.REML, σ::Union{Real,Nothing}=m.optsum.sigma, thin::Int=typemax(Int), + backend::Symbol=m.optsum.backend, + optimizer::Symbol=m.optsum.optimizer, ) where {T} optsum = m.optsum # this doesn't matter for LMM, but it does for GLMM, so let's be consistent @@ -460,6 +457,8 @@ function StatsAPI.fit!( end optsum.REML = REML optsum.sigma = σ + optsum.backend = backend + optsum.optimizer = optimizer try # use explicit evaluation w/o calling opt to avoid confusing iteration count diff --git a/src/nlopt.jl b/src/nlopt.jl index 754fa69f0..801fcb8a4 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -86,7 +86,6 @@ function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; _check_nlopt_return(ret) return xmin, fmin - end function NLopt.Opt(optsum::OptSummary) From dc5ae344580319e9c63825aa397d20ed0e68ff6d Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 01:11:02 +0000 Subject: [PATCH 14/27] use scaling --- ext/MixedModelsPRIMAExt.jl | 32 ++++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 6fa5f830e..37081e20c 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -1,7 +1,9 @@ module MixedModelsPRIMAExt using MixedModels -using MixedModels: ProgressMeter, ProgressUnknown, objective! +using MixedModels: Statistics +using MixedModels: ProgressMeter, ProgressUnknown, objective!, _objective! +using LinearAlgebra: PosDefException using PRIMA: PRIMA function MixedModels.prfit!(m::LinearMixedModel; @@ -63,14 +65,13 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; optsum.fmin = info.fx optsum.returnvalue = Symbol(info.status) _check_prima_return(info) - return optsum.final, optsum.fmin end -function optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; - progress::Bool=true, thin::Int=tyepmax(Int), - fast::Bool=false, verbose::Bool=false, nAGQ=1, - kwargs...) +function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; + progress::Bool=true, thin::Int=tyepmax(Int), + fast::Bool=false, verbose::Bool=false, nAGQ=1, + kwargs...) optsum = m.optsum fitlog = optsum.fitlog prog = ProgressUnknown(; desc="Minimizing", showspeed=true) @@ -98,9 +99,24 @@ function optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; empty!(fitlog) push!(fitlog, (copy(optsum.initial), optsum.finitial)) maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) + scale = if fast + nothing + else + # scale by the standard deviation of the columns of the fixef model matrix + # when including the fixef in the nonlinear opt + sc = [map(std, eachcol(modelmatrix(m))); fill(1, length(m.θ))] + for (i, x) in enumerate(sc) + # for nearly constant things, e.g. intercept, we don't want to scale to zero... + # also, since we're scaling the _parameters_ and not the data, + # we need to invert the scale + sc[i] = ifelse(iszero(x), one(x), inv(x)) + end + sc + end info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; xl=optsum.lowerbd, maxfun, - optsum.rhoend, optsum.rhobeg) + optsum.rhoend, optsum.rhobeg, + scale) ProgressMeter.finish!(prog) optsum.feval = info.nf @@ -108,7 +124,7 @@ function optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; optsum.returnvalue = Symbol(info.status) _check_prima_return(info) - return xmin, fmin + return optsum.final, optsum.fmin end function _check_prima_return(info::PRIMA.Info) From 8d41563490033a764e6d66c91454b6f483c629a2 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 01:44:02 +0000 Subject: [PATCH 15/27] updated show methods --- src/mimeshow.jl | 20 +++++++++-------- test/mime.jl | 35 ++++++++++++++++++++--------- test/prima.jl | 60 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 96 insertions(+), 19 deletions(-) diff --git a/src/mimeshow.jl b/src/mimeshow.jl index cd6a78162..92bd3b578 100644 --- a/src/mimeshow.jl +++ b/src/mimeshow.jl @@ -179,21 +179,23 @@ function _markdown(m::MixedModel) end function _markdown(s::OptSummary) + optimizer_settings = [["Optimizer", "`$(s.optimizer)`"], + ["Backend", "`$(s.backend)`"], + ["Lower bounds", string(s.lowerbd)]] + + for param in opt_params(Val(s.backend)) + push!(optimizer_settings, [string(param), string(getfield(s, param))]) + end + rows = [ ["", ""], ["**Initialization**", ""], ["Initial parameter vector", string(s.initial)], ["Initial objective value", string(s.finitial)], ["**Optimizer settings** ", ""], - ["Optimizer (from NLopt)", "`$(s.optimizer)`"], - ["Lower bounds", string(s.lowerbd)], - ["`ftol_rel`", string(s.ftol_rel)], - ["`ftol_abs`", string(s.ftol_abs)], - ["`xtol_rel`", string(s.xtol_rel)], - ["`xtol_abs`", string(s.xtol_abs)], - ["`initial_step`", string(s.initial_step)], - ["`maxfeval`", string(s.maxfeval)], - ["`maxtime`", string(s.maxtime)], + optimizer_settings..., + ["xtol_zero_abs", string(s.xtol_zero_abs)], + ["ftol_zero_abs", string(s.ftol_zero_abs)], ["**Result**", ""], ["Function evaluations", string(s.feval)], ["Final parameter vector", "$(round.(s.final; digits=4))"], diff --git a/test/mime.jl b/test/mime.jl index 21c6d188b..fd8cd2e09 100644 --- a/test/mime.jl +++ b/test/mime.jl @@ -120,19 +120,34 @@ lrt = likelihoodratiotest(fm0, fm1) fm1.optsum.finitial = 1784.642296192471 fm1.optsum.final = [0.9292, 0.0182, 0.2226] fm1.optsum.fmin =1751.9393444647023 - out = sprint(show, mime, fm1.optsum) + out = sprint(show, mime, fm1.optsum) @test startswith(out,""" -| | | -|:------------------------ |:--------------------------- | -| **Initialization** | | -| Initial parameter vector | [1.0, 0.0, 1.0] | -| Initial objective value | 1784.642296192471 | -| **Optimizer settings** | | -| Optimizer (from NLopt) | `LN_BOBYQA` | -| Lower bounds | [0.0, -Inf, 0.0] |""") + | | | + |:------------------------ |:--------------------------- | + | **Initialization** | | + | Initial parameter vector | [1.0, 0.0, 1.0] | + | Initial objective value | 1784.642296192471 | + | **Optimizer settings** | | + | Optimizer | `LN_BOBYQA` | + | Backend | `nlopt` | + | Lower bounds | [0.0, -Inf, 0.0] | + | ftol_rel | 1.0e-12 | + | ftol_abs | 1.0e-8 | + | xtol_rel | 0.0 | + | xtol_abs | [1.0e-10, 1.0e-10, 1.0e-10] | + | initial_step | [0.75, 1.0, 0.75] | + | maxfeval | -1 | + | maxtime | -1.0 | + | xtol_zero_abs | 0.001 | + | ftol_zero_abs | 1.0e-5 | + | **Result** | | + | Function evaluations | 1 | + | Final parameter vector | [0.9292, 0.0182, 0.2226] | + | Final objective value | 1751.9393 | + | Return code | `FTOL_REACHED` | + """) end - @testset "varcorr" begin @test sprint(show, mime, VarCorr(fm1)) == """ diff --git a/test/prima.jl b/test/prima.jl index 18e622c00..38b8556d8 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -24,3 +24,63 @@ prmodel.optsum.backend = :prima fit!(prmodel) @test isapprox(loglikelihood(model), loglikelihood(prmodel)) end + +@testset "optsum show" begin + model = first(models(:contra)) + prmodel = unfit!(deepcopy(model)) + fit!(prmodel; optimizer=:bobyqa, backend=:prima) + @test isapprox(loglikelihood(model), loglikelihood(prmodel)) + + optsum = deepcopy(prmodel.optsum) + optsum.final = [0.2612] + optsum.finitial = 2595.85 + optsum.fmin = 2486.42 + optsum.feval = 17 + + out = sprint(show, MIME("text/plain"), optsum) + expected = """ + Initial parameter vector: [1.0] + Initial objective value: 2595.85 + + Backend: prima + Optimizer: bobyqa + Lower bounds: [0.0] + rhobeg: 1.0 + rhoend: 1.0e-6 + maxfeval: -1 + + Function evaluations: 17 + xtol_zero_abs: 0.001 + ftol_zero_abs: 1.0e-5 + Final parameter vector: [0.2612] + Final objective value: 2486.42 + Return code: SMALL_TR_RADIUS + """ + + @test startswith(out, expected) + + out = sprint(show, MIME("text/markdown"), optsum) + expected = """ + | | | + |:------------------------ |:----------------- | + | **Initialization** | | + | Initial parameter vector | [1.0] | + | Initial objective value | 2595.85 | + | **Optimizer settings** | | + | Optimizer | `bobyqa` | + | Backend | `prima` | + | Lower bounds | [0.0] | + | rhobeg | 1.0 | + | rhoend | 1.0e-6 | + | maxfeval | -1 | + | xtol_zero_abs | 0.001 | + | ftol_zero_abs | 1.0e-5 | + | **Result** | | + | Function evaluations | 17 | + | Final parameter vector | [0.2612] | + | Final objective value | 2486.42 | + | Return code | `SMALL_TR_RADIUS` | + """ + + @test startswith(out, expected) +end From b136447cd568158a3a2d966a4c8b684a2f396ee0 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 02:32:38 +0000 Subject: [PATCH 16/27] cleanup --- ext/MixedModelsPRIMAExt.jl | 1 - src/generalizedlinearmixedmodel.jl | 2 +- src/nlopt.jl | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 37081e20c..49b6c8d7b 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -97,7 +97,6 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) scale = if fast nothing diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 3958e430f..5b07b9005 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -295,7 +295,7 @@ function StatsAPI.fit!( end loglength = length(fitlog) if xmin ≠ xmin_ - if (zeroobj = obj(xmin_, T[])) ≤ (fmin + optsum.ftol_zero_abs) + if (zeroobj = objective!(m, xmin_; nAGQ, fast, verbose)) ≤ (fmin + optsum.ftol_zero_abs) fmin = zeroobj copyto!(xmin, xmin_) elseif length(fitlog) > loglength diff --git a/src/nlopt.jl b/src/nlopt.jl index 801fcb8a4..95cf96c40 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -77,7 +77,6 @@ function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; NLopt.min_objective!(opt, obj) optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial)) ProgressMeter.finish!(prog) From 4613eb2bc37fd1711c3be462fe89a2973b4f5fb2 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 02:39:41 +0000 Subject: [PATCH 17/27] test fix --- test/prima.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/test/prima.jl b/test/prima.jl index 38b8556d8..373560779 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -1,11 +1,8 @@ using PRIMA -using MixedModels: prfit! -using MixedModels: dataset +using MixedModels: prfit!, unfit!, dataset include("modelcache.jl") -# model = first(models(:sleepstudy)) - @testset "formula($model)" for model in models(:sleepstudy) prmodel = prfit!(LinearMixedModel(formula(model), dataset(:sleepstudy)); progress=false) @@ -25,10 +22,14 @@ prmodel.optsum.backend = :prima @test isapprox(loglikelihood(model), loglikelihood(prmodel)) end -@testset "optsum show" begin - model = first(models(:contra)) +@testset "GLMM + optsum show" begin + model = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), + dataset(:contra), Binomial(); progress=false) prmodel = unfit!(deepcopy(model)) - fit!(prmodel; optimizer=:bobyqa, backend=:prima) + fit!(prmodel; optimizer=:bobyqa, backend=:prima, progress=false) + @test isapprox(loglikelihood(model), loglikelihood(prmodel)) + refit!(prmodel; fast=true, progress=false) + refit!(model; fast=true, progress=false) @test isapprox(loglikelihood(model), loglikelihood(prmodel)) optsum = deepcopy(prmodel.optsum) From 1c947f3afbab4cc6e7c51bf524b3c58b6bca83e2 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 03:03:39 +0000 Subject: [PATCH 18/27] fix?? --- test/pirls.jl | 2 +- test/pls.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/test/pirls.jl b/test/pirls.jl index 1e90d9a56..391e61f0e 100644 --- a/test/pirls.jl +++ b/test/pirls.jl @@ -37,7 +37,7 @@ end thin = 5 gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin) fitlog = gm0.optsum.fitlog - @test length(fitlog) == (div(gm0.optsum.feval, thin) + 1) # for the initial value + @test length(fitlog) == length(0:thin:gm0.optsum.feval) @test first(fitlog) == (gm0.optsum.initial, gm0.optsum.finitial) @test gm0.lowerbd == zeros(1) @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001) diff --git a/test/pls.jl b/test/pls.jl index 95b1e877b..3b1ccfe6e 100644 --- a/test/pls.jl +++ b/test/pls.jl @@ -160,7 +160,7 @@ end @test d == length(first(fitlogtbl.θ)) thin = 2 refit!(fm1; REML=false, progress=false, thin) - @test length(fitlog) == (div(fm1.optsum.feval, thin) + 1) # for the initial value + @test length(fitlog) == length(0:thin:fm1.optsum.feval, thin) @test first(fitlog) == (fm1.optsum.initial, fm1.optsum.finitial) end @testset "profile" begin From af5f896fa5b34df06944fd537e8421b243e052b8 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 03:08:46 +0000 Subject: [PATCH 19/27] :facepalm: --- test/pls.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pls.jl b/test/pls.jl index 3b1ccfe6e..e630ffd8e 100644 --- a/test/pls.jl +++ b/test/pls.jl @@ -160,7 +160,7 @@ end @test d == length(first(fitlogtbl.θ)) thin = 2 refit!(fm1; REML=false, progress=false, thin) - @test length(fitlog) == length(0:thin:fm1.optsum.feval, thin) + @test length(fitlog) == length(0:thin:fm1.optsum.feval) @test first(fitlog) == (fm1.optsum.initial, fm1.optsum.finitial) end @testset "profile" begin From 705542c3188517720413f47e8f42f05da3a7d69c Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 05:07:11 +0000 Subject: [PATCH 20/27] make thinning a noop --- NEWS.md | 4 +++- ext/MixedModelsPRIMAExt.jl | 28 ++++++++-------------- src/bootstrap.jl | 4 ++-- src/generalizedlinearmixedmodel.jl | 13 +++++------ src/linearmixedmodel.jl | 23 +++++++++++-------- src/nlopt.jl | 37 ++++++++++++------------------ src/optsummary.jl | 11 ++++----- src/profile/fixefpr.jl | 2 +- src/profile/profile.jl | 2 +- src/profile/sigmapr.jl | 4 ++-- test/pirls.jl | 12 ++++++---- test/pls.jl | 21 ++++++++--------- 12 files changed, 75 insertions(+), 86 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5d8c8d614..4969b8b37 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,9 @@ MixedModels v4.30.0 Release Notes ============================== - Refactor calls to backend optimizer to make it easier to add and use different optimization backends. - The structure of `OptSummary` has been expanded and `prfit!` has been updated to use this new structure. [#802] + The structure of `OptSummary` has been accordlingly expanded and `prfit!` has been updated to use this new structure. [#802] +- Make the `thin` argument to `fit!` a no-op. It complicated several bits of logic without having any real performance benefit in the majority of cases. This argument has been replaced with a `fitlog::Bool=false` that determines whether a log is kept.[#802] + MixedModels v4.29.1 Release Notes ============================== diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 49b6c8d7b..dfb30b422 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -25,16 +25,14 @@ push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) const PRIMABackend = Val{:prima} function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; - progress::Bool=true, thin::Int=tyepmax(Int)) + progress::Bool=true, fitlog::Bool=false, kwargs...) optsum = m.optsum prog = ProgressUnknown(; desc="Minimizing", showspeed=true) - # start from zero for the initial call to obj before optimization - iter = 0 - fitlog = optsum.fitlog + fitlog && empty!(optsum.fitlog) function obj(x) - iter += 1 - val = if isone(iter) && x == optsum.initial + val = if x == optsum.initial + # fast path since we've already evaluated the initial value optsum.finitial else try @@ -50,12 +48,10 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; end end progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - !isone(iter) && iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + fitlog && push!(optsum.fitlog, (copy(x), val)) return val end - empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) info = prima_optimizer!(Val(optsum.optimizer), obj, optsum.final; xl=optsum.lowerbd, maxfun, @@ -69,15 +65,13 @@ function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; end function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; - progress::Bool=true, thin::Int=tyepmax(Int), + progress::Bool=true, fitlog::Bool=false, fast::Bool=false, verbose::Bool=false, nAGQ=1, kwargs...) optsum = m.optsum - fitlog = optsum.fitlog prog = ProgressUnknown(; desc="Minimizing", showspeed=true) - # start from zero for the initial call to obj before optimization - iter = 0 - fitlog = optsum.fitlog + fitlog && empty!(opstum.fitlog) + function obj(x) val = try _objective!(m, x, Val(fast); verbose, nAGQ) @@ -85,18 +79,16 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; # this allows us to recover from models where e.g. the link isn't # as constraining as it should be ex isa Union{PosDefException,DomainError} || rethrow() - iter == 1 && rethrow() + x == optsum.initial && rethrow() m.optsum.finitial end - iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + fitlog && push!(optsum.fitlog, (copy(x), val)) verbose && println(round(val; digits=5), " ", x) progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - iter += 1 return val end optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) - empty!(fitlog) maxfun = optsum.maxfeval > 0 ? optsum.maxfeval : 500 * length(optsum.initial) scale = if fast nothing diff --git a/src/bootstrap.jl b/src/bootstrap.jl index e59809bab..0bd711387 100644 --- a/src/bootstrap.jl +++ b/src/bootstrap.jl @@ -234,7 +234,7 @@ function parametricbootstrap( end β = convert(Vector{T}, β) θ = convert(Vector{T}, θ) - # scratch -- note that this is the length of the unpivoted coef vector + # scratch -- note that this is the length of the unpivoted coef vector βsc = coef(morig) θsc = zeros(ftype, length(θ)) p = length(βsc) @@ -254,7 +254,7 @@ function parametricbootstrap( ) samp = replicate(n; progress) do simulate!(rng, m; β, σ, θ) - refit!(m; progress=false) + refit!(m; progress=false, fitlog=false) ( objective=ftype.(m.objective), σ=ismissing(m.σ) ? missing : ftype(m.σ), diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 5b07b9005..b9fd880fd 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -219,7 +219,8 @@ and it may not be shown at all for models that are optimized quickly. If `verbose` is `true`, then both the intermediate results of both the nonlinear optimization and PIRLS are also displayed on standard output. -At every `thin`th iteration is recorded in `fitlog`, optimization progress is saved in `m.optsum.fitlog`. +The `thin` argument is ignored: it had no impact on the final model fit and the logic around +thinning the `fitlog` was needlessly complicated for a trivial performance gain. By default, the starting values for model fitting are taken from a (non mixed, i.e. marginal ) GLM fit. Experience with larger datasets (many thousands of @@ -247,6 +248,7 @@ function StatsAPI.fit!( nAGQ::Integer=1, progress::Bool=true, thin::Int=typemax(Int), + fitlog::Bool=false, init_from_lmm=Set(), backend::Symbol=m.optsum.backend, optimizer::Symbol=m.optsum.optimizer, @@ -283,8 +285,7 @@ function StatsAPI.fit!( optsum.backend = backend optsum.optimizer = optimizer - xmin, fmin = optimize!(m; progress, thin, fast, verbose, nAGQ) - fitlog = optsum.fitlog + xmin, fmin = optimize!(m; progress, fitlog, fast, verbose, nAGQ) ## check if very small parameter values bounded below by zero can be set to zero xmin_ = copy(xmin) @@ -293,16 +294,14 @@ function StatsAPI.fit!( xmin_[i] = zero(T) end end - loglength = length(fitlog) if xmin ≠ xmin_ if (zeroobj = objective!(m, xmin_; nAGQ, fast, verbose)) ≤ (fmin + optsum.ftol_zero_abs) fmin = zeroobj copyto!(xmin, xmin_) - elseif length(fitlog) > loglength - # remove unused extra log entry - pop!(fitlog) + fitlog && push!(optsum.fitlog, (copy(xmin), fmin)) end end + ## ensure that the parameter values saved in m are xmin objective!(m, xmin; fast, verbose, nAGQ) optsum.final = xmin diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 2c37fb2aa..65b2e0d56 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -427,14 +427,15 @@ end """ fit!(m::LinearMixedModel; progress::Bool=true, REML::Bool=m.optsum.REML, σ::Union{Real, Nothing}=m.optsum.sigma, - thin::Int=typemax(Int)) + thin::Int=typemax(Int), + fitlog::Bool=true) Optimize the objective of a `LinearMixedModel`. When `progress` is `true` a `ProgressMeter.ProgressUnknown` display is shown during the optimization of the objective, if the optimization takes more than one second or so. -At every `thin`th iteration is recorded in `fitlog`, optimization progress is -saved in `m.optsum.fitlog`. +The `thin` argument is ignored: it had no impact on the final model fit and the logic around +thinning the `fitlog` was needlessly complicated for a trivial performance gain. """ function StatsAPI.fit!( m::LinearMixedModel{T}; @@ -442,6 +443,7 @@ function StatsAPI.fit!( REML::Bool=m.optsum.REML, σ::Union{Real,Nothing}=m.optsum.sigma, thin::Int=typemax(Int), + fitlog::Bool=false, backend::Symbol=m.optsum.backend, optimizer::Symbol=m.optsum.optimizer, ) where {T} @@ -476,8 +478,8 @@ function StatsAPI.fit!( optsum.finitial = objective!(m, optsum.initial) end - xmin, fmin = optimize!(m; progress, thin) - fitlog = optsum.fitlog + xmin, fmin = optimize!(m; progress, fitlog) + ## check if small non-negative parameter values can be set to zero xmin_ = copy(xmin) lb = optsum.lowerbd @@ -486,18 +488,19 @@ function StatsAPI.fit!( xmin_[i] = zero(T) end end - loglength = length(fitlog) if xmin ≠ xmin_ if (zeroobj = objective!(m, xmin_)) ≤ (fmin + optsum.ftol_zero_abs) fmin = zeroobj copyto!(xmin, xmin_) - elseif length(fitlog) > loglength - # remove unused extra log entry - pop!(fitlog) + fitlog && push!(optsum.fitlog, (copy(xmin), fmin)) end end + + # unlike GLMM we don't need to populate optsum.finitial here + # because we do that during the initial guess and rescale check + ## ensure that the parameter values saved in m are xmin - updateL!(setθ!(m, xmin)) + objective!(m)(xmin) optsum.final = xmin optsum.fmin = fmin diff --git a/src/nlopt.jl b/src/nlopt.jl index 95cf96c40..8df9078ed 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -3,16 +3,16 @@ push!(OPTIMIZATION_BACKENDS, :nlopt) const NLoptBackend = Val{:nlopt} function optimize!(m::LinearMixedModel, ::NLoptBackend; - progress::Bool=true, thin::Int=tyepmax(Int), + progress::Bool=true, fitlog::Bool=false, kwargs...) optsum = m.optsum - opt = Opt(optsum) - prog = ProgressUnknown(; desc="Minimizing", showspeed=true) + fitlog && empty!(optsum.fitlog) + function obj(x, g) isempty(g) || throw(ArgumentError("g should be empty for this objective")) - iter += 1 - val = if isone(iter) && x == optsum.initial + val = if x == optsum.initial + # fast path since we've already evaluated the initial value optsum.finitial else try @@ -28,16 +28,12 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; end end progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - !isone(iter) && iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + fitlog && push!(optsum.fitlog, (copy(x), val)) return val end - NLopt.min_objective!(opt, obj) - # start from zero for the initial call to obj before optimization - iter = 0 - fitlog = optsum.fitlog - empty!(fitlog) - push!(fitlog, (copy(optsum.initial), optsum.finitial)) + opt = Opt(optsum) + NLopt.min_objective!(opt, obj) fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial)) ProgressMeter.finish!(prog) optsum.feval = opt.numevals @@ -47,15 +43,13 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; end function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; - progress::Bool=true, thin::Int=tyepmax(Int), + progress::Bool=true, fitlog::Bool=false, fast::Bool=false, verbose::Bool=false, nAGQ=1, kwargs...) optsum = m.optsum - fitlog = optsum.fitlog prog = ProgressUnknown(; desc="Minimizing", showspeed=true) - # start from zero for the initial call to obj before optimization - iter = 0 - fitlog = optsum.fitlog + fitlog && empty!(optsum.fitlog) + function obj(x, g) isempty(g) || throw(ArgumentError("g should be empty for this objective")) val = try @@ -64,19 +58,18 @@ function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; # this allows us to recover from models where e.g. the link isn't # as constraining as it should be ex isa Union{PosDefException,DomainError} || rethrow() - iter == 1 && rethrow() - m.optsum.finitial + x == optsum.initial && rethrow() + optsum.finitial end - iszero(rem(iter, thin)) && push!(fitlog, (copy(x), val)) + fitlog && push!(optsum.fitlog, (copy(x), val)) verbose && println(round(val; digits=5), " ", x) progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) - iter += 1 return val end + opt = Opt(optsum) NLopt.min_objective!(opt, obj) optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) - empty!(fitlog) fmin, xmin, ret = NLopt.optimize(opt, copyto!(optsum.final, optsum.initial)) ProgressMeter.finish!(prog) diff --git a/src/optsummary.jl b/src/optsummary.jl index f5c1dd930..c9e79f978 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -76,7 +76,7 @@ Base.@kwdef mutable struct OptSummary{T<:AbstractFloat} rhoend::T = rhobeg / 1_000_000 # not SVector because we would need to parameterize on size (which breaks GLMM) - fitlog::Vector{Tuple{Vector{T},T}} = [(initial, fmin)] + fitlog::Vector{Tuple{Vector{T},T}} = Vector{Tuple{Vector{T},T}}() nAGQ::Int = 1 REML::Bool = false sigma::Union{T,Nothing} = nothing @@ -97,13 +97,10 @@ end Return `s.fitlog` as a `Tables.columntable`. When `stack` is false (the default), there will be 3 columns in the result: -- `iter`: the sample number +- `iter`: the iteration number - `objective`: the value of the objective at that sample - `θ`: the parameter vector at that sample -(The term `sample` here refers to the fact that when the `thin` argument to the `fit` or -`refit!` call is greater than 1 only a subset of the iterations have results recorded.) - When `stack` is true, there will be 4 columns: `iter`, `objective`, `par`, and `value` where `value` is the stacked contents of the `θ` vectors (the equivalent of `vcat(θ...)`) and `par` is a vector of parameter numbers. @@ -174,8 +171,8 @@ with the second argument of type `::Val{:backend_name}` and adding `:backend_nam be defined, which returns the optimization parameters (e.g. `xtol_abs` or `rho_end`) used by the backend. -Common keyword arguments are `progress` to show a progress meter and `thin` to control thinning -of the fitlog. +Common keyword arguments are `progress` to show a progress meter as well as +`nAQG` and `fast` for GLMMs. """ function optimize! end diff --git a/src/profile/fixefpr.jl b/src/profile/fixefpr.jl index 98ab98166..ead371321 100644 --- a/src/profile/fixefpr.jl +++ b/src/profile/fixefpr.jl @@ -54,7 +54,7 @@ function betaprofile!( pr::FeProfile{T}, tc::TableColumns{T}, βⱼ::T, j::Integer, obj::T, neg::Bool ) where {T} prm = pr.m - refit!(prm, mul!(copyto!(prm.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false) + refit!(prm, mul!(copyto!(prm.y, pr.y₀), pr.xⱼ, βⱼ, -1, 1); progress=false, fitlog=false) (; positions, v) = tc v[1] = (-1)^neg * sqrt(prm.objective - obj) getθ!(view(v, positions[:θ]), prm) diff --git a/src/profile/profile.jl b/src/profile/profile.jl index 82d71944c..1364a2c4c 100644 --- a/src/profile/profile.jl +++ b/src/profile/profile.jl @@ -33,7 +33,7 @@ Profiling starts at the parameter estimate and continues until reaching a parame value of ζ exceeds `threshold`. """ function profile(m::LinearMixedModel; threshold=4) - isfitted(m) || refit!(m) + isfitted(m) || refit!(m; progress=false, fitlog=false) final = copy(m.optsum.final) profile = try tc = TableColumns(m) diff --git a/src/profile/sigmapr.jl b/src/profile/sigmapr.jl index 2fbf9cf22..ef8c05f7d 100644 --- a/src/profile/sigmapr.jl +++ b/src/profile/sigmapr.jl @@ -13,7 +13,7 @@ function refitσ!( m::LinearMixedModel{T}, σ, tc::TableColumns{T}, obj::T, neg::Bool ) where {T} m.optsum.sigma = σ - refit!(m; progress=false) + refit!(m; progress=false, fitlog=false) return mkrow!(tc, m, (neg ? -one(T) : one(T)) * sqrt(m.objective - obj)) end @@ -26,7 +26,7 @@ function _facsz(m::LinearMixedModel{T}, σ::T, obj::T) where {T} i64 = T(inv(64)) expi64 = exp(i64) # help the compiler infer it is a constant m.optsum.sigma = σ * expi64 - return exp(i64 / (2 * sqrt(refit!(m; progress=false).objective - obj))) + return exp(i64 / (2 * sqrt(refit!(m; progress=false, fitlog=false).objective - obj))) end """ diff --git a/test/pirls.jl b/test/pirls.jl index 391e61f0e..a00e4b21c 100644 --- a/test/pirls.jl +++ b/test/pirls.jl @@ -34,11 +34,15 @@ end @testset "contra" begin contra = dataset(:contra) - thin = 5 - gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, thin) + gm0 = fit(MixedModel, first(gfms[:contra]), contra, Bernoulli(); fast=true, progress=false, fitlog=true) fitlog = gm0.optsum.fitlog - @test length(fitlog) == length(0:thin:gm0.optsum.feval) - @test first(fitlog) == (gm0.optsum.initial, gm0.optsum.finitial) + @test length(fitlog) == gm0.optsum.feval + @test first(fitlog)[1] == gm0.optsum.initial + @test first(fitlog)[2] ≈ gm0.optsum.finitial + # XXX this should be exact equality and it is indeed when stepping through manually + # but not when run via Pkg.test(). I have no idea why. + @test last(fitlog)[1] ≈ gm0.optsum.final + @test last(fitlog)[2] ≈ gm0.optsum.fmin @test gm0.lowerbd == zeros(1) @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001) @test !issingular(gm0) diff --git a/test/pls.jl b/test/pls.jl index e630ffd8e..ae6409547 100644 --- a/test/pls.jl +++ b/test/pls.jl @@ -143,8 +143,7 @@ end # since we're caching the fits, we should get it back to being correctly fitted # we also take this opportunity to test fitlog @testset "fitlog" begin - thin = 1 - refit!(fm1; REML=false, progress=false, thin) + refit!(fm1; REML=false, progress=false, fitlog=true) fitlog = fm1.optsum.fitlog fitlogtbl = columntable(fm1.optsum) @test length(fitlogtbl) == 3 @@ -158,10 +157,6 @@ end d, r = divrem(length(first(fitlogstackedtbl)), length(first(fitlogtbl))) @test iszero(r) @test d == length(first(fitlogtbl.θ)) - thin = 2 - refit!(fm1; REML=false, progress=false, thin) - @test length(fitlog) == length(0:thin:fm1.optsum.feval) - @test first(fitlog) == (fm1.optsum.initial, fm1.optsum.finitial) end @testset "profile" begin dspr01 = profile(only(models(:dyestuff))) @@ -518,14 +513,18 @@ end @test "BlkDiag" in Set(split(String(take!(io)), r"\s+")) @testset "optsumJSON" begin - fm = last(models(:sleepstudy)) + fm = refit!(last(models(:sleepstudy)), progress=false, fitlog=true) # using a IOBuffer for saving JSON saveoptsum(seekstart(io), fm) m = LinearMixedModel(fm.formula, MixedModels.dataset(:sleepstudy)) restoreoptsum!(m, seekstart(io)) - @test loglikelihood(fm) ≈ loglikelihood(m) - @test bic(fm) ≈ bic(m) - @test coef(fm) ≈ coef(m) + @test length(fm.optsum.fitlog) == length(m.optsum.fitlog) + + # try it out with an empty fitlog + empty!(fm.optsum.fitlog) + saveoptsum(seekstart(io), fm) + restoreoptsum!(m, seekstart(io)) + @test length(fm.optsum.fitlog) == length(m.optsum.fitlog) == 0 fm_mod = deepcopy(fm) fm_mod.optsum.fmin += 1 @@ -702,7 +701,7 @@ end contrasts = Dict(:item => Grouping(), :subj => Grouping(), :prec => EffectsCoding(; base="maintain"), :spkr => EffectsCoding(), :load => EffectsCoding()) kbf03 = @formula rt_trunc ~ 1+prec+spkr+load+(1+prec|item)+(1|subj) - kbpr03 = profile(fit(MixedModel, kbf03, MixedModels.dataset(:kb07); contrasts, thin=1, progress=false)) + kbpr03 = profile(fit(MixedModel, kbf03, MixedModels.dataset(:kb07); contrasts, progress=false)) @test length(Tables.columnnames(kbpr03.tbl)) == 15 @test length(Tables.rows(kbpr03.tbl)) > 200 end From 92e2f1757d328ec24fc80b4cec7a28a33c48f723 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 05:07:40 +0000 Subject: [PATCH 21/27] Update NEWS.md Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 4969b8b37..fbace7377 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,7 @@ MixedModels v4.30.0 Release Notes ============================== - Refactor calls to backend optimizer to make it easier to add and use different optimization backends. - The structure of `OptSummary` has been accordlingly expanded and `prfit!` has been updated to use this new structure. [#802] + The structure of `OptSummary` has been accordingly expanded and `prfit!` has been updated to use this new structure. [#802] - Make the `thin` argument to `fit!` a no-op. It complicated several bits of logic without having any real performance benefit in the majority of cases. This argument has been replaced with a `fitlog::Bool=false` that determines whether a log is kept.[#802] From 49d547899412bffb7ec6b8d8410d79543010ccfc Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 05:15:29 +0000 Subject: [PATCH 22/27] BlueStyle --- src/generalizedlinearmixedmodel.jl | 19 +++++++++++++------ src/linearmixedmodel.jl | 14 +++++++------- src/mimeshow.jl | 4 ++-- src/nlopt.jl | 15 ++++++++------- src/optsummary.jl | 3 +-- src/serialization.jl | 1 - 6 files changed, 31 insertions(+), 25 deletions(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index b9fd880fd..3679dfc32 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -295,7 +295,8 @@ function StatsAPI.fit!( end end if xmin ≠ xmin_ - if (zeroobj = objective!(m, xmin_; nAGQ, fast, verbose)) ≤ (fmin + optsum.ftol_zero_abs) + if (zeroobj = objective!(m, xmin_; nAGQ, fast, verbose)) ≤ + (fmin + optsum.ftol_zero_abs) fmin = zeroobj copyto!(xmin, xmin_) fitlog && push!(optsum.fitlog, (copy(xmin), fmin)) @@ -517,7 +518,9 @@ function StatsAPI.loglikelihood(m::GeneralizedLinearMixedModel{T}) where {T} end # Base.Fix1 doesn't forward kwargs -objective!(m::GeneralizedLinearMixedModel; fast=false, kwargs...) = x -> _objective!(m, x, Val(fast); kwargs...) +function objective!(m::GeneralizedLinearMixedModel; fast=false, kwargs...) + return x -> _objective!(m, x, Val(fast); kwargs...) +end function objective!(m::GeneralizedLinearMixedModel{T}, x; fast=false, kwargs...) where {T} return _objective!(m, x, Val(fast); kwargs...) @@ -526,12 +529,16 @@ end # normally, it doesn't make sense to move a simple branch to dispatch # HOWEVER, this wins up getting called in optimization a lot and # moving this to a realization here allows us to avoid dynamic dispatch on setθ! / setθβ! -function _objective!(m::GeneralizedLinearMixedModel{T}, x, ::Val{true}; nAGQ=1, verbose=false) where {T} - deviance(pirls!(setθ!(m, x), true, verbose), nAGQ) +function _objective!( + m::GeneralizedLinearMixedModel{T}, x, ::Val{true}; nAGQ=1, verbose=false +) where {T} + return deviance(pirls!(setθ!(m, x), true, verbose), nAGQ) end -function _objective!(m::GeneralizedLinearMixedModel{T}, x, ::Val{false}; nAGQ=1, verbose=false) where {T} - deviance(pirls!(setβθ!(m, x), false, verbose), nAGQ) +function _objective!( + m::GeneralizedLinearMixedModel{T}, x, ::Val{false}; nAGQ=1, verbose=false +) where {T} + return deviance(pirls!(setβθ!(m, x), false, verbose), nAGQ) end function Base.propertynames(m::GeneralizedLinearMixedModel, private::Bool=false) diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 65b2e0d56..6d571a6b5 100644 --- a/src/linearmixedmodel.jl +++ b/src/linearmixedmodel.jl @@ -206,13 +206,13 @@ function StatsAPI.fit( end function StatsAPI.fit(::Type{LinearMixedModel}, - f::FormulaTerm, - tbl::Tables.ColumnTable; - wts=[], - contrasts=Dict{Symbol,Any}(), - σ=nothing, - amalgamate=true, - kwargs...) + f::FormulaTerm, + tbl::Tables.ColumnTable; + wts=[], + contrasts=Dict{Symbol,Any}(), + σ=nothing, + amalgamate=true, + kwargs...) lmm = LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate) return fit!(lmm; kwargs...) end diff --git a/src/mimeshow.jl b/src/mimeshow.jl index 92bd3b578..31db00694 100644 --- a/src/mimeshow.jl +++ b/src/mimeshow.jl @@ -180,8 +180,8 @@ end function _markdown(s::OptSummary) optimizer_settings = [["Optimizer", "`$(s.optimizer)`"], - ["Backend", "`$(s.backend)`"], - ["Lower bounds", string(s.lowerbd)]] + ["Backend", "`$(s.backend)`"], + ["Lower bounds", string(s.lowerbd)]] for param in opt_params(Val(s.backend)) push!(optimizer_settings, [string(param), string(getfield(s, param))]) diff --git a/src/nlopt.jl b/src/nlopt.jl index 8df9078ed..4f41c3951 100644 --- a/src/nlopt.jl +++ b/src/nlopt.jl @@ -3,8 +3,8 @@ push!(OPTIMIZATION_BACKENDS, :nlopt) const NLoptBackend = Val{:nlopt} function optimize!(m::LinearMixedModel, ::NLoptBackend; - progress::Bool=true, fitlog::Bool=false, - kwargs...) + progress::Bool=true, fitlog::Bool=false, + kwargs...) optsum = m.optsum prog = ProgressUnknown(; desc="Minimizing", showspeed=true) fitlog && empty!(optsum.fitlog) @@ -43,9 +43,9 @@ function optimize!(m::LinearMixedModel, ::NLoptBackend; end function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; - progress::Bool=true, fitlog::Bool=false, - fast::Bool=false, verbose::Bool=false, nAGQ=1, - kwargs...) + progress::Bool=true, fitlog::Bool=false, + fast::Bool=false, verbose::Bool=false, nAGQ=1, + kwargs...) optsum = m.optsum prog = ProgressUnknown(; desc="Minimizing", showspeed=true) fitlog && empty!(optsum.fitlog) @@ -101,7 +101,6 @@ function NLopt.Opt(optsum::OptSummary) return opt end - const _NLOPT_FAILURE_MODES = [ :FAILURE, :INVALID_ARGS, @@ -118,4 +117,6 @@ function _check_nlopt_return(ret, failure_modes=_NLOPT_FAILURE_MODES) end end -opt_params(::NLoptBackend) = (:ftol_rel, :ftol_abs, :xtol_rel, :xtol_abs, :initial_step, :maxfeval, :maxtime) +function opt_params(::NLoptBackend) + return (:ftol_rel, :ftol_abs, :xtol_rel, :xtol_abs, :initial_step, :maxfeval, :maxtime) +end diff --git a/src/optsummary.jl b/src/optsummary.jl index c9e79f978..e9c3ef11d 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -129,7 +129,7 @@ function Base.show(io::IO, ::MIME"text/plain", s::OptSummary) for param in opt_params(Val(s.backend)) println(io, rpad(string(param, ":"), length("Initial parameter vector: ")), - getfield(s, param)) + getfield(s, param)) end println(io) println(io, "Function evaluations: ", s.feval) @@ -176,7 +176,6 @@ Common keyword arguments are `progress` to show a progress meter as well as """ function optimize! end - """ opt_params(::Val{backend}) diff --git a/src/serialization.jl b/src/serialization.jl index e0e44d9f0..ba470c07f 100644 --- a/src/serialization.jl +++ b/src/serialization.jl @@ -123,7 +123,6 @@ function restoreoptsum!(ops::OptSummary{T}, dict::AbstractDict) where {T} return ops end - StructTypes.StructType(::Type{<:OptSummary}) = StructTypes.Mutable() StructTypes.excludes(::Type{<:OptSummary}) = (:lowerbd,) From d105dea3d7e14e1091a64918140b468e7b9f3dad Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 05:38:59 +0000 Subject: [PATCH 23/27] docs update --- docs/src/optimization.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/optimization.md b/docs/src/optimization.md index 100b18914..f73a8e887 100644 --- a/docs/src/optimization.md +++ b/docs/src/optimization.md @@ -190,11 +190,11 @@ DisplayAs.Text(ans) # hide ``` More detailed information about the intermediate steps of the nonlinear optimizer can be obtained the `fitlog` field. -By default, `fitlog` contains entries for only the initial and final steps, but additional information about every nth step can be obtained with the `thin` keyword-argument to `fit`, `fit!` and `refit!`: +By default, `fitlog` is not populated, but passing the keyword argument `fitlog=true` to `fit!` and `refit!` will result in it being populated with the values obtained at each step of optimization: ```@example Main -refit!(fm2; thin=1) -fm2.optsum.fitlog[1:10] +refit!(fm2; fitlog=true) +first(fm2.optsum.fitlog, 5) DisplayAs.Text(ans) # hide ``` From fe30063bbb8f6bdb58567989c49d697936393040 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 05:47:09 +0000 Subject: [PATCH 24/27] more coverage --- ext/MixedModelsPRIMAExt.jl | 2 +- test/prima.jl | 10 +++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index dfb30b422..8366db2bb 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -120,7 +120,7 @@ end function _check_prima_return(info::PRIMA.Info) if !PRIMA.issuccess(info) - @warn "PRIMA optimization failure: $(ret)\n$(PRIMA.reason(info))" + @warn "PRIMA optimization failure: $(info.status)\n$(PRIMA.reason(info))" end return nothing diff --git a/test/prima.jl b/test/prima.jl index 373560779..bdc37d257 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -18,10 +18,18 @@ prmodel.optsum.backend = :prima @testset "$optimizer" for optimizer in (:cobyla, :lincoa) MixedModels.unfit!(prmodel) prmodel.optsum.optimizer = optimizer - fit!(prmodel) + fit!(prmodel; progress=false, fitlog=false) @test isapprox(loglikelihood(model), loglikelihood(prmodel)) end +@testset "failure" begin + MixedModels.unfit!(prmodel) + prmodel.optsum.optimizer = :bobyqa + prmodel.optsum.maxfeval = 5 + @test_logs((:warn, r"PRIMA optimization failure"), + fit!(prmodel; progress=false, fitlog=false)) +end + @testset "GLMM + optsum show" begin model = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), dataset(:contra), Binomial(); progress=false) From 8549f8a16aa0740ff54425f8e0e2be1f128eb460 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 18:18:24 +0000 Subject: [PATCH 25/27] Update src/generalizedlinearmixedmodel.jl --- src/generalizedlinearmixedmodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generalizedlinearmixedmodel.jl b/src/generalizedlinearmixedmodel.jl index 3679dfc32..7b9d9bdcb 100644 --- a/src/generalizedlinearmixedmodel.jl +++ b/src/generalizedlinearmixedmodel.jl @@ -527,7 +527,7 @@ function objective!(m::GeneralizedLinearMixedModel{T}, x; fast=false, kwargs...) end # normally, it doesn't make sense to move a simple branch to dispatch -# HOWEVER, this wins up getting called in optimization a lot and +# HOWEVER, this winds up getting called in optimization a lot and # moving this to a realization here allows us to avoid dynamic dispatch on setθ! / setθβ! function _objective!( m::GeneralizedLinearMixedModel{T}, x, ::Val{true}; nAGQ=1, verbose=false From e437477522b1139351c09080e5670498958ba39a Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 12:20:48 -0600 Subject: [PATCH 26/27] move push! to init --- ext/MixedModelsPRIMAExt.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 8366db2bb..1fb55f11d 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -6,6 +6,13 @@ using MixedModels: ProgressMeter, ProgressUnknown, objective!, _objective! using LinearAlgebra: PosDefException using PRIMA: PRIMA +function __init__() + push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) + return nothing +end + +const PRIMABackend = Val{:prima} + function MixedModels.prfit!(m::LinearMixedModel; kwargs...) @@ -20,10 +27,6 @@ prima_optimizer!(::Val{:bobyqa}, args...; kwargs...) = PRIMA.bobyqa!(args...; kw prima_optimizer!(::Val{:cobyla}, args...; kwargs...) = PRIMA.cobyla!(args...; kwargs...) prima_optimizer!(::Val{:lincoa}, args...; kwargs...) = PRIMA.lincoa!(args...; kwargs...) -push!(MixedModels.OPTIMIZATION_BACKENDS, :prima) - -const PRIMABackend = Val{:prima} - function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; progress::Bool=true, fitlog::Bool=false, kwargs...) optsum = m.optsum From 4bbd1f51f7d82a4ddbe964b2040adad4aa003574 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Sun, 26 Jan 2025 12:37:58 -0600 Subject: [PATCH 27/27] rm redundant newline --- NEWS.md | 1 - 1 file changed, 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index fbace7377..bcd740e34 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,7 +4,6 @@ MixedModels v4.30.0 Release Notes The structure of `OptSummary` has been accordingly expanded and `prfit!` has been updated to use this new structure. [#802] - Make the `thin` argument to `fit!` a no-op. It complicated several bits of logic without having any real performance benefit in the majority of cases. This argument has been replaced with a `fitlog::Bool=false` that determines whether a log is kept.[#802] - MixedModels v4.29.1 Release Notes ============================== - Populate `optsum` in `prfit!` call. [#801]