populate optsum in prfit!

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v4.29.1 Release Notes
+==============================
+- Populate `optsum` in `prfit!` call. [#801]
+
 MixedModels v4.29.0 Release Notes
 ==============================
 - Testbed for experimental support for using PRIMA as an optimization backend introduced via the experimental `prfit!` function. [#799]
@@ -590,3 +594,4 @@ Package dependencies
 [#792]: https://github.com/JuliaStats/MixedModels.jl/issues/792
 [#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

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>"]
-version = "4.29.0"
+version = "4.29.1"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

ext/MixedModelsPRIMAExt.jl

@@ -1,10 +1,54 @@
 module MixedModelsPRIMAExt
 
 using MixedModels: MixedModels, LinearMixedModel, objective!
+using MixedModels: ProgressMeter, ProgressUnknown
 using PRIMA: PRIMA
 
-function MixedModels.prfit!(m::LinearMixedModel)
-    PRIMA.bobyqa!(objective!(m), copy(m.optsum.initial); xl=m.optsum.lowerbd)
+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)
+
+    prog = ProgressUnknown(; desc="Minimizing", showspeed=true)
+    # start from zero for the initial call to obj before optimization
+    iter = 0
+    fitlog = empty!(optsum.fitlog)
+    function obj(x)
+        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)])
+        if isone(iter) || iszero(rem(iter, thin))
+            push!(fitlog, (copy(x), val))
+        end
+        return val
+    end
+
+    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
 end
 

src/prima.jl

@@ -13,6 +13,25 @@ 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

test/prima.jl

@@ -7,7 +7,7 @@ include("modelcache.jl")
 # model = first(models(:sleepstudy))
 
 @testset "formula($model)" for model in models(:sleepstudy)
-    prmodel = prfit!(LinearMixedModel(formula(model), dataset(:sleepstudy)))
+    prmodel = prfit!(LinearMixedModel(formula(model), dataset(:sleepstudy)); progress=false)
 
     @test isapprox(loglikelihood(model), loglikelihood(prmodel))
 end