remove lowerbd field from OptSummary

.github/workflows/current.yml

@@ -19,12 +19,13 @@ jobs:
   ci:
     runs-on: ${{ matrix.os }}
     strategy:
+      fail-fast: false
       matrix:
         julia-version: [1]
         # julia-arch: [x64]
         os:
-         - ubuntu-22.04
-         - macOS-14 # apple silicon!
+         - ubuntu-24.04
+         - macOS-15 # apple silicon!
     steps:
       - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v2

NEWS.md

@@ -2,7 +2,7 @@ MixedModels v5.0.0 Release Notes
 ==============================
 - Optimization is now performed _without constraints_. In a post-fitting step, the Cholesky factor is canonicalized to have non-negative diagonal elements. [#840]
 - The default optimizer has changed to NLopt's implementation of NEWUOA where possible. NLopt's implementation fails on 1-dimensional problems, so in the case of a single, scalar random effect, BOBYQA is used instead. In the future, the default optimizer backend will likely change to PRIMA and NLopt support will be moved to an extension. Blocking this change in backend is an issue with PRIMA.jl when running in VSCode's built-in REPL on Linux. [#840]
-- [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`.
+- [BREAKING] Support for constrained optimization has been completely removed, i.e. the field `lowerbd` has been removed from `OptSummary`. [#849]
 - [BREAKING] A fitlog is always kept -- the deprecated keyword argument `thin` has been removed as has the `fitlog` keyword argument. [#850]
 - The fitlog is now stored as Tables.jl-compatible column table. [#850]
 - Internal code around the default optimizer has been restructured. In particular, the NLopt backend has been moved to a submodule, which will make it easier to move it to an extension if we promote another backend to the default. [#853]

ext/MixedModelsPRIMAExt.jl

@@ -105,7 +105,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::PRIMABackend;
         sc
     end
     info = _optimizer!(Val(optsum.optimizer), obj, optsum.final;
-        xl=optsum.lowerbd, maxfun,
+        maxfun,
         optsum.rhoend, optsum.rhobeg,
         scale)
     ProgressMeter.finish!(prog)

src/MixedModelsNLoptExt.jl

@@ -97,8 +97,7 @@ function MixedModels.optimize!(m::GeneralizedLinearMixedModel, ::NLoptBackend;
 end
 
 function NLopt.Opt(optsum::OptSummary)
-    lb = optsum.lowerbd
-    n = length(lb)
+    n = length(optsum.initial)
 
     if optsum.optimizer == :LN_NEWUOA && isone(n) # :LN_NEWUOA doesn't allow n == 1
         optsum.optimizer = :LN_BOBYQA
@@ -110,11 +109,10 @@ function NLopt.Opt(optsum::OptSummary)
     if length(optsum.xtol_abs) == n  # 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)   # use unconstrained optimization even for :LN_BOBYQA
     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))
+        optsum.initial_step = NLopt.initial_step(opt, optsum.initial)
     else
         NLopt.initial_step!(opt, optsum.initial_step)
     end

src/bootstrap.jl

@@ -12,7 +12,7 @@ Object returned by `parametericbootstrap` with fields
 - `fits`: the parameter estimates from the bootstrap replicates as a vector of named tuples.
 - `λ`: `Vector{LowerTriangular{T,Matrix{T}}}` containing copies of the λ field from `ReMat` model terms
 - `inds`: `Vector{Vector{Int}}` containing copies of the `inds` field from `ReMat` model terms
-- `lowerbd`: `Vector{T}` containing the vector of lower bounds (corresponds to the identically named field of [`OptSummary`](@ref))
+- `lowerbd`: `Vector{T}` containing the vector of lower bounds (corresponds to the `lowerbd(model)` of the original model)
 - `fcnames`: NamedTuple whose keys are the grouping factor names and whose values are the column names
 
 The schema of `fits` is, by default,
@@ -34,7 +34,7 @@ struct MixedModelBootstrap{T<:AbstractFloat} <: MixedModelFitCollection{T}
     fits::Vector
     λ::Vector{Union{LowerTriangular{T},Diagonal{T}}}
     inds::Vector{Vector{Int}}
-    lowerbd::Vector{T}
+    lowerbd::Vector{T} # we need to store this explicitly because we no longer have access to the ReMats
     fcnames::NamedTuple
 end
 
@@ -112,7 +112,7 @@ function restorereplicates(
         samp,
         map(vv -> T.(vv), m.λ), # also does a deepcopy if no type conversion is necessary
         getfield.(m.reterms, :inds),
-        T.(m.optsum.lowerbd[1:length(first(samp).θ)]),
+        T.(lowerbd(m)),
         NamedTuple{Symbol.(fnames(m))}(map(t -> Tuple(t.cnames), m.reterms)),
     )
 end
@@ -175,6 +175,8 @@ function Base.show(io::IO, mime::MIME"text/plain", x::MixedModelBootstrap)
     return nothing
 end
 
+lowerbd(x::MixedModelFitCollection) = x.lowerbd
+
 """
     parametricbootstrap([rng::AbstractRNG], nsamp::Integer, m::MixedModel{T}, ftype=T;
         β = fixef(m), σ = m.σ, θ = m.θ, progress=true, optsum_overrides=(;))
@@ -243,8 +245,6 @@ function parametricbootstrap(
     for (key, val) in pairs(optsum_overrides)
         setfield!(m.optsum, key, val)
     end
-    # this seemed to slow things down?!
-    # _copy_away_from_lowerbd!(m.optsum.initial, morig.optsum.final, m.lowerbd; incr=0.05)
 
     β_names = Tuple(Symbol.(coefnames(morig)))
 
@@ -267,7 +267,7 @@ function parametricbootstrap(
         samp,
         map(vv -> ftype.(vv), morig.λ), # also does a deepcopy if no type conversion is necessary
         getfield.(morig.reterms, :inds),
-        ftype.(morig.optsum.lowerbd[1:length(first(samp).θ)]),
+        ftype.(lowerbd(morig)),
         NamedTuple{Symbol.(fnames(morig))}(map(t -> Tuple(t.cnames), morig.reterms)),
     )
 end
@@ -421,7 +421,7 @@ function issingular(
     bsamp::MixedModelFitCollection; atol::Real=0, rtol::Real=atol > 0 ? 0 : √eps()
 )
     return map(bsamp.θ) do θ
-        return _issingular(bsamp.lowerbd, θ; atol, rtol)
+        return _issingular(lowerbd(bsamp), θ; atol, rtol)
     end
 end
 

src/generalizedlinearmixedmodel.jl

@@ -30,7 +30,6 @@
 - `theta`: synonym for `θ`
 - `beta`: synonym for `β`
 - `σ` or `sigma`: common scale parameter (value is `NaN` for distributions without a scale parameter)
-- `lowerbd`: vector of lower bounds on the combined elements of `β` and `θ`
 - `formula`, `trms`, `A`, `L`, and `optsum`: fields of the `LMM` field
 - `X`: fixed-effects model matrix
 - `y`: response vector
@@ -280,7 +279,6 @@
     end
 
     if !fast
-        optsum.lowerbd = vcat(fill!(similar(β), T(-Inf)), optsum.lowerbd)
         optsum.initial = vcat(β, lm.optsum.final)
         optsum.final = copy(optsum.initial)
     end
@@ -296,8 +294,9 @@
 
     ## check if very small parameter values bounded below by zero can be set to zero
     xmin_ = copy(xmin)
+    lb = fast ? lowerbd(m) : vcat(zero(β), lowerbd(m))
     for i in eachindex(xmin_)
-        if iszero(optsum.lowerbd[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
+        if iszero(lb[i]) && zero(T) < xmin_[i] < optsum.xtol_zero_abs
             xmin_[i] = zero(T)
         end
     end
@@ -524,6 +523,8 @@
     return accum - (mapreduce(u -> sum(abs2, u), +, m.u) + logdet(m)) / 2
 end
 
+lowerbd(m::GeneralizedLinearMixedModel) = lowerbd(m.LMM)
+
 # Base.Fix1 doesn't forward kwargs
 function objective!(m::GeneralizedLinearMixedModel; fast=false, kwargs...)
     return x -> _objective!(m, x, Val(fast); kwargs...)
@@ -560,7 +561,6 @@
         :sigma,
         :X,
         :y,
-        :lowerbd,
         :objective,
         :σρs,
         :σs,
@@ -733,7 +733,7 @@
     io::IO, ::MIME"text/plain", m::GeneralizedLinearMixedModel{T,D}
 ) where {T,D}
     if m.optsum.feval < 0
         @warn("Model has not been fit")
         return nothing
     end
     nAGQ = m.LMM.optsum.nAGQ
@@ -785,11 +785,6 @@
 
     reterms = model.LMM.reterms
     optsum = model.LMM.optsum
-    # we need to reset optsum so that it
-    # plays nice with the modifications fit!() does
-    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)      # probably don't need this anymore - now trivial with all elements = -Inf
-    # for variances (bounded at zero), we have ones, while
-    # for everything else (bounded at -Inf), we have zeros
     optsum.initial = map(x -> T(x[2] == x[3]), model.LMM.parmap)
     optsum.final = copy(optsum.initial)
     optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
@@ -839,7 +834,7 @@
 end
 
 # delegate GLMM method to LMM field
-for f in (:feL, :fetrm, :fixefnames, :(LinearAlgebra.logdet), :lowerbd, :PCA, :rePCA)
+for f in (:feL, :fetrm, :fixefnames, :(LinearAlgebra.logdet), :PCA, :rePCA)
     @eval begin
         $f(m::GeneralizedLinearMixedModel) = $f(m.LMM)
     end

src/linearmixedmodel.jl

@@ -24,7 +24,6 @@ Linear mixed-effects model representation
 * `σ` or `sigma`: current value of the standard deviation of the per-observation noise
 * `b`: random effects on the original scale, as a vector of matrices
 * `u`: random effects on the orthogonal scale, as a vector of matrices
-* `lowerbd`: lower bounds on the elements of θ
 * `X`: the fixed-effects model matrix
 * `y`: the response vector
 """
@@ -173,9 +172,8 @@ function LinearMixedModel(
     reweight!.(reterms, Ref(sqrtwts))
     reweight!(Xy, sqrtwts)
     A, L = createAL(reterms, Xy)
-    lbd = foldl(vcat, lowerbd(c) for c in reterms)
     θ = foldl(vcat, getθ(c) for c in reterms)
-    optsum = OptSummary(θ, lbd)
+    optsum = OptSummary(θ)
     optsum.sigma = isnothing(σ) ? nothing : T(σ)
     return LinearMixedModel(
         form,
@@ -662,8 +660,6 @@ function Base.getproperty(m::LinearMixedModel{T}, s::Symbol) where {T}
         stderror(m)
     elseif s == :u
         ranef(m; uscale=true)
-    elseif s == :lowerbd
-        m.optsum.lowerbd
     elseif s == :X
         modelmatrix(m)
     elseif s == :y
@@ -791,7 +787,7 @@ function StatsAPI.loglikelihood(m::LinearMixedModel)
     return -objective(m) / 2
 end
 
-lowerbd(m::LinearMixedModel) = m.optsum.lowerbd
+lowerbd(m::LinearMixedModel) = foldl(vcat, lowerbd(c) for c in m.reterms)
 
 function mkparmap(reterms::Vector{<:AbstractReMat{T}}) where {T}
     parmap = NTuple{3,Int}[]
@@ -856,7 +852,6 @@ function Base.propertynames(m::LinearMixedModel, private::Bool=false)
         :sigmarhos,
         :b,
         :u,
-        :lowerbd,
         :X,
         :y,
         :corr,

src/mimeshow.jl

@@ -184,7 +184,7 @@ 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))])

src/optsummary.jl

@@ -8,7 +8,6 @@ Summary of an optimization
 ## 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
 * `final`: a copy of the final parameter values from the optimization
 * `fmin`: the final value of the objective
 * `feval`: the number of function evaluations
@@ -65,7 +64,6 @@ Similarly, the list of applicable optimization parameters can be inspected with
 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
@@ -100,11 +98,9 @@ end
 
 function OptSummary(
     initial::Vector{T},
-    lowerbd::Vector{S},
     optimizer::Symbol=:LN_NEWUOA; kwargs...,
-) where {T<:AbstractFloat,S<:AbstractFloat}
-    TS = promote_type(T, S)
-    return OptSummary{TS}(; initial, lowerbd, optimizer, kwargs...)
+) where {T<:AbstractFloat}
+    return OptSummary{T}(; initial, optimizer, kwargs...)
 end
 
 """
@@ -141,7 +137,6 @@ function Base.show(io::IO, ::MIME"text/plain", s::OptSummary)
     println(io)
     println(io, "Backend:                  ", s.backend)
     println(io, "Optimizer:                ", s.optimizer)
-    println(io, "Lower bounds:             ", s.lowerbd)
 
     for param in opt_params(Val(s.backend))
         println(io, rpad(string(param, ":"), length("Initial parameter vector: ")),

src/profile/fixefpr.jl

@@ -44,9 +44,7 @@ function FeProfile(m::LinearMixedModel, tc::TableColumns, j::Integer)
         LinearMixedModel(y₀ - xⱼ * m.β[j], feterm, reterms, m.formula); progress=false
     )
     # not sure this next call makes sense - should the second argument be m.optsum.final?
-    _copy_away_from_lowerbd!(
-        mnew.optsum.initial, mnew.optsum.final, mnew.lowerbd; incr=0.05
-    )
+    copyto!(mnew.optsum.initial, mnew.optsum.final)
     return FeProfile(mnew, tc, y₀, xⱼ, j)
 end
 

src/profile/sigmapr.jl

@@ -43,7 +43,7 @@ function profileσ(m::LinearMixedModel{T}, tc::TableColumns{T}; threshold=4) whe
         throw(ArgumentError("Can't profile σ, which is fixed at $(optsum.sigma)"))
     θ = copy(optsum.final)
     θinitial = copy(optsum.initial)
-    _copy_away_from_lowerbd!(optsum.initial, optsum.final, optsum.lowerbd)
+    copyto!(optsum.initial, optsum.final)
     obj = optsum.fmin
     σ = m.σ
     pnm = (p=:σ,)

src/profile/thetapr.jl

@@ -9,7 +9,6 @@ Return an `OptSummary` with the `j`'th component of the parameter omitted.
 function optsumj(os::OptSummary, j::Integer)
     return OptSummary(
         deleteat!(copy(os.final), j),
-        deleteat!(copy(os.lowerbd), j),
         os.optimizer;
         os.backend,
     )

src/profile/utilities.jl

@@ -145,35 +145,6 @@ function ρvals!(
     return v
 end
 
-"""
-    _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
-
-Replace `sink[i]` by `max(source[i], bd[i] + incr)`.  When `bd[i] == -Inf` this simply copies `source[i]`.
-"""
-function _copy_away_from_lowerbd!(sink, source, bd; incr=0.01)
-    for i in eachindex(sink, source, bd)
-        @inbounds sink[i] = max(source[i], bd[i] + incr)
-    end
-    return sink
-end
-
-#=  # It appears that this method is not used
-"""
-    stepsize(tbl::Vector{NamedTuple}, resp::Symbol, pred::Symbol; rev::Bool=false)
-
-Return the stepsize from the last value of `tbl.pred` to increase `resp` by approximately 0.5
-"""
-function stepsize(tbl::Vector{<:NamedTuple}, resp::Symbol, pred::Symbol)
-    ntbl = length(tbl)
-    lm1tbl = tbl[ntbl - 1]
-    x1 = getproperty(lm1tbl, pred)
-    y1 = getproperty(lm1tbl, resp)
-    x2 = getproperty(last(tbl), pred)
-    y2 = getproperty(last(tbl), resp)
-    return (x2 - x1) / (2 * (y2 - y1))
-end
-=#
-
 function isnondecreasing(spl::SplineInterpolation)
     return all(≥(0), (Derivative(1) * spl).(spl.x))
 end

src/profile/vcpr.jl

@@ -36,10 +36,11 @@ function profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T
     m = val.m
     (; λ, σ, β, optsum, parmap, reterms) = m
     isnothing(optsum.sigma) || throw(ArgumentError("Can't profile vc's when σ is fixed"))
-    (; initial, final, fmin, lowerbd) = optsum
-    lowerbd .+= T(nzlb)                       # lower bounds must be > 0 b/c θ's occur in denominators
+    (; initial, final, fmin) = optsum
+    # lowerbd .+= T(nzlb)                       # lower bounds must be > 0 b/c θ's occur in denominators
     saveinitial = copy(initial)
-    copyto!(initial, max.(final, lowerbd))
+    # copyto!(initial, max.(final, lowerbd))
+    copyto!(initial, final)
     zetazero = mkrow!(tc, m, zero(T))         # parameter estimates
     vcnms = filter(keys(first(val.tbl))) do sym
         str = string(sym)
@@ -84,7 +85,6 @@ function profileσs!(val::NamedTuple, tc::TableColumns{T}; nzlb=1.0e-8) where {T
     end
     copyto!(final, initial)
     copyto!(initial, saveinitial)
-    lowerbd .-= T(nzlb)
     optsum.sigma = nothing
     updateL!(setθ!(m, final))
     return val