experimental harness for PRIMA (#799)

* experimental harness for PRIMA

---------

Co-authored-by: Douglas Bates <[email protected]>

Author: palday

NEWS.md

@@ -1,3 +1,7 @@
+MixedModels v4.29.0 Release Notes
+==============================
+- Testbed for experimental support for using PRIMA as an optimization backend introduced via the experimental `prfit!` function. [#799]
+
 MixedModels v4.28.0 Release Notes
 ==============================
 - `GeneralizedLinearMixedModel` now attempts to fall back to very constrained variance values when the default initial values result in a non positive semidefinite covariance matrix. [#792]
@@ -585,3 +589,4 @@ Package dependencies
 [#791]: https://github.com/JuliaStats/MixedModels.jl/issues/791
 [#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

Project.toml

@@ -30,6 +30,12 @@ StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
+[weakdeps]
+PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
+
+[extensions]
+MixedModelsPRIMAExt = ["PRIMA"]
+
 [compat]
 Aqua = "0.8"
 Arrow = "1, 2"
@@ -46,6 +52,7 @@ LinearAlgebra = "1"
 Markdown = "1"
 MixedModelsDatasets = "0.1"
 NLopt = "0.6, 1"
+PRIMA = "0.2"
 PooledArrays = "0.5, 1"
 PrecompileTools = "1"
 ProgressMeter = "1.7"
@@ -63,16 +70,17 @@ Suppressor = "0.2"
 Tables = "1"
 Test = "1"
 TypedTables = "1"
-julia = "1.8"
+julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test"]
+test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "PRIMA", "StableRNGs", "Suppressor", "Test"]

ext/MixedModelsPRIMAExt.jl

@@ -0,0 +1,11 @@
+module MixedModelsPRIMAExt
+
+using MixedModels: MixedModels, LinearMixedModel, objective!
+using PRIMA: PRIMA
+
+function MixedModels.prfit!(m::LinearMixedModel)
+    PRIMA.bobyqa!(objective!(m), copy(m.optsum.initial); xl=m.optsum.lowerbd)
+    return m
+end
+
+end # module

src/MixedModels.jl

@@ -205,6 +205,7 @@ include("grouping.jl")
 include("mimeshow.jl")
 include("serialization.jl")
 include("profile/profile.jl")
+include("prima.jl")
 
 # COV_EXCL_START
 @setup_workload begin

src/prima.jl

@@ -0,0 +1,21 @@
+
+"""
+    prfit!(m::LinearMixedModel)
+
+Fit a mixed model using the [PRIMA](https://github.com/libprima/PRIMA.jl) implementation
+of the BOBYQA optimizer.
+
+!!! warning "Experimental feature"
+    This function is an experimental feature that will go away in the future.
+    Do **not** rely on it, unless you are willing to pin the precise MixedModels.jl
+    version. The purpose of the function is to provide the MixedModels developers
+    a chance to explore the performance of the PRIMA implementation without the large
+    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 "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
+    by the user.
+"""
+function prfit! end

test/prima.jl

@@ -0,0 +1,13 @@
+using PRIMA
+using MixedModels: prfit!
+using MixedModels: dataset
+
+include("modelcache.jl")
+
+# model = first(models(:sleepstudy))
+
+@testset "formula($model)" for model in models(:sleepstudy)
+    prmodel = prfit!(LinearMixedModel(formula(model), dataset(:sleepstudy)))
+
+    @test isapprox(loglikelihood(model), loglikelihood(prmodel))
+end

test/runtests.jl

@@ -7,8 +7,6 @@ using Test
 import InteractiveUtils: versioninfo
 import LinearAlgebra: BLAS
 
-using Base: splat # necessary for Julia 1.8 compat
-
 # there seem to be processor-specific issues and knowing this is helpful
 @info sprint(versioninfo)
 @info BLAS.get_config()
@@ -45,3 +43,5 @@ include("mime.jl")
 include("optsummary.jl")
 include("predict.jl")
 include("sigma.jl")
+
+@testset "PRIMA" include("prima.jl")