diff --git a/NEWS.md b/NEWS.md index c34ee21d1..bcd740e34 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +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 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] @@ -595,3 +601,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" 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 ``` diff --git a/ext/MixedModelsPRIMAExt.jl b/ext/MixedModelsPRIMAExt.jl index 6e80d8995..1fb55f11d 100644 --- a/ext/MixedModelsPRIMAExt.jl +++ b/ext/MixedModelsPRIMAExt.jl @@ -1,27 +1,41 @@ module MixedModelsPRIMAExt -using MixedModels: MixedModels, LinearMixedModel, objective! -using MixedModels: ProgressMeter, ProgressUnknown +using MixedModels +using MixedModels: Statistics +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; - 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...) + + MixedModels.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...) + +function MixedModels.optimize!(m::LinearMixedModel, ::PRIMABackend; + 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 = empty!(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 @@ -37,19 +51,84 @@ 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)) + fitlog && push!(optsum.fitlog, (copy(x), val)) + return val + end + + 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 optsum.final, optsum.fmin +end + +function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend; + 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!(opstum.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() + x == optsum.initial && rethrow() + m.optsum.finitial end + fitlog && push!(optsum.fitlog, (copy(x), val)) + verbose && println(round(val; digits=5), " ", x) + progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) return val end + optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) + 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, + scale) 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 + _check_prima_return(info) + + return optsum.final, optsum.fmin +end + +function _check_prima_return(info::PRIMA.Info) + if !PRIMA.issuccess(info) + @warn "PRIMA optimization failure: $(info.status)\n$(PRIMA.reason(info))" + end + + return nothing end +MixedModels.opt_params(::PRIMABackend) = (:rhobeg, :rhoend, :maxfeval) + end # module 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/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 d6b998b36..7b9d9bdcb 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,7 +248,10 @@ 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, ) where {T} β = copy(m.β) θ = copy(m.θ) @@ -277,35 +281,12 @@ 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 - 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) + + optsum.backend = backend + optsum.optimizer = optimizer + + 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) for i in eachindex(xmin_) @@ -313,24 +294,20 @@ function StatsAPI.fit!( xmin_[i] = zero(T) end 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 - # 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 - 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 +517,30 @@ 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 +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...) +end + +# normally, it doesn't make sense to move a simple branch to dispatch +# 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 +) 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} + return deviance(pirls!(setβθ!(m, x), false, verbose), nAGQ) +end + function Base.propertynames(m::GeneralizedLinearMixedModel, private::Bool=false) return ( :A, diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl index 14f297835..6d571a6b5 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, @@ -206,21 +205,16 @@ function StatsAPI.fit( ) end -function StatsAPI.fit( - ::Type{LinearMixedModel}, +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 - ) + kwargs...) + lmm = LinearMixedModel(f, tbl; contrasts, wts, σ, amalgamate) + return fit!(lmm; kwargs...) end function _offseterr() @@ -433,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}; @@ -448,6 +443,9 @@ 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} optsum = m.optsum # this doesn't matter for LMM, but it does for GLMM, so let's be consistent @@ -459,55 +457,29 @@ 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 - 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) + optsum.backend = backend + optsum.optimizer = optimizer + try # use explicit evaluation w/o calling opt to avoid confusing iteration count - optsum.finitial = objective(updateL!(setθ!(m, optsum.initial))) + 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. - """ + 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))) + 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) + + xmin, fmin = optimize!(m; progress, fitlog) + ## check if small non-negative parameter values can be set to zero xmin_ = copy(xmin) lb = optsum.lowerbd @@ -516,24 +488,22 @@ function StatsAPI.fit!( xmin_[i] = zero(T) end 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 - # 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.feval = opt.numevals optsum.final = xmin optsum.fmin = fmin - optsum.returnvalue = ret - _check_nlopt_return(ret) return m end @@ -828,25 +798,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/mimeshow.jl b/src/mimeshow.jl index cd6a78162..31db00694 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/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..4f41c3951 --- /dev/null +++ b/src/nlopt.jl @@ -0,0 +1,122 @@ +push!(OPTIMIZATION_BACKENDS, :nlopt) + +const NLoptBackend = Val{:nlopt} + +function optimize!(m::LinearMixedModel, ::NLoptBackend; + progress::Bool=true, fitlog::Bool=false, + kwargs...) + optsum = m.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")) + val = if x == optsum.initial + # fast path since we've already evaluated the initial value + 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)]) + fitlog && push!(optsum.fitlog, (copy(x), val)) + return val + end + + 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 + optsum.returnvalue = ret + _check_nlopt_return(ret) + return xmin, fmin +end + +function optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend; + 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) + + 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() + x == optsum.initial && rethrow() + optsum.finitial + end + fitlog && push!(optsum.fitlog, (copy(x), val)) + verbose && println(round(val; digits=5), " ", x) + progress && ProgressMeter.next!(prog; showvalues=[(:objective, val)]) + return val + end + + opt = Opt(optsum) + NLopt.min_objective!(opt, obj) + optsum.finitial = _objective!(m, optsum.initial, Val(fast); verbose, nAGQ) + 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 + +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 864f09cd6..e9c3ef11d 100644 --- a/src/optsummary.jl +++ b/src/optsummary.jl @@ -1,61 +1,82 @@ """ 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`. +* `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-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 * `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 + fitlog::Vector{Tuple{Vector{T},T}} = Vector{Tuple{Vector{T},T}}() nAGQ::Int = 1 REML::Bool = false sigma::Union{T,Nothing} = nothing @@ -76,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. @@ -105,66 +123,65 @@ 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) -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 + println(io, "Return code: ", s.returnvalue) + return nothing 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...) + +""" + 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 parameters (e.g. `xtol_abs` or `rho_end`) used +by the backend. + +Common keyword arguments are `progress` to show a progress meter as well as +`nAQG` and `fast` for GLMMs. +""" +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 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/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/src/serialization.jl b/src/serialization.jl index 708504c18..ba470c07f 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,9 @@ 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) 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/pirls.jl b/test/pirls.jl index 1e90d9a56..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) == (div(gm0.optsum.feval, thin) + 1) # for the initial value - @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 95b1e877b..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) == (div(fm1.optsum.feval, thin) + 1) # for the initial value - @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 diff --git a/test/prima.jl b/test/prima.jl index ff32cf03d..bdc37d257 100644 --- a/test/prima.jl +++ b/test/prima.jl @@ -1,13 +1,95 @@ 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) @test isapprox(loglikelihood(model), loglikelihood(prmodel)) + @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; 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) + prmodel = unfit!(deepcopy(model)) + 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) + 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