Merge branch 'main' of github.com:JuliaStats/MixedModels.jl into pa/b…

…enchmark_ci2

.github/workflows/CompatHelper.yml

@@ -9,8 +9,8 @@ jobs:
     steps:
       - name: Pkg.add("CompatHelper")
         run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
-      - name: CompatHelper.main()
+      - name: CompatHelper.main(; bump_version=true)
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
-        run: julia -e 'using CompatHelper; CompatHelper.main()'
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

.github/workflows/spellcheck.yml

@@ -0,0 +1,21 @@
+# adapted from https://github.com/JuliaDocs/Documenter.jl/blob/master/.github/workflows/SpellCheck.yml
+# see docs at https://github.com/crate-ci/typos
+name: Spell Check
+on: [pull_request]
+
+jobs:
+  typos-check:
+    name: Spell Check with Typos
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout Actions Repository
+        uses: actions/checkout@v4
+      - name: Check spelling
+        uses: crate-ci/typos@master
+        with:
+            config: _typos.toml
+            write_changes: true
+      - uses: reviewdog/action-suggester@v1
+        with:
+          tool_name: Typos
+          fail_on_error: true

NEWS.md

@@ -1,8 +1,23 @@
+MixedModels v4.22.5 Release Notes
+==============================
+* Use `muladd` where possible to enable fused multiply-add (FMA) on architectures with hardware support. FMA will generally improve computational speed and gives more accurate rounding. [#740]
+* Replace broadcasted lambda with explicit loop and use `one`. This may result in a small performance improvement. [#738]
+
+MixedModels v4.22.4 Release Notes
+==============================
+* Switch to explicit imports from all included packages (i.e. replace `using Foo` by `using Foo: Foo, bar, baz`) [#748]
+* Reset parameter values before a `deepcopy` in a test (doesn't change test result) [#744]
+
+MixedModels v4.22.3 Release Notes
+==============================
+* Comment out calls to `@debug` [#733]
+* Update package versions in compat and change `Aqua.test_all` argument name [#733]
+
 MixedModels v4.22.0 Release Notes
 ==============================
 * Support for equal-tail confidence intervals for `MixedModelBootstrap`. [#715]
 * Basic `show` methods for `MixedModelBootstrap` and `MixedModelProfile`. [#715]
-* The `hide_progress` keyword argument to `parametricbootstrap` is now deprecated. Users should isntead use `progress` (which is consistent with e.g. `fit`). [#717]
+* The `hide_progress` keyword argument to `parametricbootstrap` is now deprecated. Users should instead use `progress` (which is consistent with e.g. `fit`). [#717]
 
 MixedModels v4.21.0 Release Notes
 ==============================
@@ -484,3 +499,8 @@ Package dependencies
 [#709]: https://github.com/JuliaStats/MixedModels.jl/issues/709
 [#715]: https://github.com/JuliaStats/MixedModels.jl/issues/715
 [#717]: https://github.com/JuliaStats/MixedModels.jl/issues/717
+[#733]: https://github.com/JuliaStats/MixedModels.jl/issues/733
+[#738]: https://github.com/JuliaStats/MixedModels.jl/issues/738
+[#740]: https://github.com/JuliaStats/MixedModels.jl/issues/740
+[#744]: https://github.com/JuliaStats/MixedModels.jl/issues/744
+[#748]: https://github.com/JuliaStats/MixedModels.jl/issues/748

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>", "Jose Bayoan Santiago Calderon <[email protected]>"]
-version = "4.22.1"
+version = "4.22.5"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
@@ -30,31 +30,47 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
 
 [compat]
+Aqua = "0.8"
 Arrow = "1, 2"
-BSplineKit = "0.14, 0.15, 0.16"
+BSplineKit = "0.14, 0.15, 0.16, 0.17"
 DataAPI = "1"
+DataFrames = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
+ExplicitImports = "1.3"
 GLM = "1.8.2"
+InteractiveUtils = "1"
 JSON3 = "1"
+LinearAlgebra = "1"
+Markdown = "1"
 MixedModelsDatasets = "0.1"
 NLopt = "0.6, 1"
 PooledArrays = "0.5, 1"
 PrecompileTools = "1"
 ProgressMeter = "1.7"
+Random = "1"
+SparseArrays = "1"
+StableRNGs = "0.1, 1"
 StaticArrays = "0.11, 0.12, 1"
 Statistics = "1"
 StatsAPI = "1.5"
 StatsBase = "0.31, 0.32, 0.33, 0.34"
 StatsFuns = "0.8, 0.9, 1"
 StatsModels = "0.7"
 StructTypes = "1"
+Suppressor = "0.2"
 Tables = "1"
+Test = "1"
 TypedTables = "1"
 julia = "1.8"
 
 [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"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "DataFrames"]
+test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "StableRNGs", "Suppressor", "Test"]

_typos.toml

@@ -0,0 +1,20 @@
+# https://github.com/crate-ci/typos#false-positives
+[default]
+
+[default.extend-identifiers]
+Lik = "Lik"
+missings = "missings"
+
+[default.extend-words]
+Lik = "Lik"
+missings = "missings"
+
+[type.package_toml]
+# Don't check spellings in these files
+extend-glob = ["Manifest.toml", "Project.toml"]
+check-file = false
+
+[type.bib]
+# contain lots of names, which are a great spot for false positives
+extend-glob = ["*.bib"]
+check-file = false

src/MixedModels.jl

@@ -1,35 +1,47 @@
 module MixedModels
 
-using Arrow
-using BSplineKit
-using DataAPI
-using Distributions
-using GLM
-using JSON3
-using LinearAlgebra
-using Markdown
-using NLopt
-using Random
-using PooledArrays
-using ProgressMeter
-using SparseArrays
-using StaticArrays
-using Statistics
-using StatsAPI
-using StatsBase
-using StatsModels
-using StructTypes
-using Tables
-using TypedTables
-
-using LinearAlgebra: BlasFloat, BlasReal, HermOrSym, PosDefException, copytri!
+using Arrow: Arrow
 using Base: Ryu, require_one_based_indexing
-using DataAPI: levels, refpool, refarray, refvalue
-using GLM: Link, canonicallink, linkfun, linkinv, dispersion, dispersion_parameter
+using BSplineKit: BSplineKit, BSplineOrder, Natural, Derivative, SplineInterpolation
+using BSplineKit: interpolate
+using DataAPI: DataAPI, levels, refpool, refarray, refvalue
+using Distributions: Distributions, Bernoulli, Binomial, Chisq, Distribution, Gamma
+using Distributions: InverseGaussian, Normal, Poisson, ccdf, estimate
+using GLM: GLM, GeneralizedLinearModel, IdentityLink, InverseLink, LinearModel
+using GLM: Link, LogLink, LogitLink, ProbitLink, SqrtLink
+using GLM: canonicallink, glm, linkinv, dispersion, dispersion_parameter
+using JSON3: JSON3
+using LinearAlgebra: LinearAlgebra, Adjoint, BLAS, BlasFloat, ColumnNorm
+using LinearAlgebra: Diagonal, Hermitian, HermOrSym, I, LAPACK, LowerTriangular
+using LinearAlgebra: PosDefException, SVD, SymTridiagonal, Symmetric
+using LinearAlgebra: UpperTriangular, cond, diag, diagind, dot, eigen, isdiag
+using LinearAlgebra: ldiv!, lmul!, logdet, mul!, norm, normalize, normalize!, qr
+using LinearAlgebra: rank, rdiv!, rmul!, svd, tr, tril!
+using Markdown: Markdown
 using MixedModelsDatasets: dataset, datasets
-using NLopt: Opt
-using StatsModels: TableRegressionModel
+using NLopt: NLopt, Opt, ftol_abs, ftol_rel, initial_step, maxtime, xtol_abs, xtol_rel
+using PooledArrays: PooledArrays, PooledArray
+using PrecompileTools: PrecompileTools, @setup_workload, @compile_workload
+using ProgressMeter: ProgressMeter, Progress, ProgressUnknown, finish!, next!
+using Random: Random, AbstractRNG, randn!
+using SparseArrays: SparseArrays, SparseMatrixCSC, SparseVector, dropzeros!, nnz
+using SparseArrays: nonzeros, nzrange, rowvals, sparse
+using StaticArrays: StaticArrays, SVector
+using Statistics: Statistics, mean, quantile, std
+using StatsAPI: StatsAPI, aic, aicc, bic, coef, coefnames, coeftable, confint, deviance
+using StatsAPI: dof, dof_residual, fit, fit!, fitted, isfitted, islinear, leverage
+using StatsAPI: loglikelihood, meanresponse, modelmatrix, nobs, predict, r2, residuals
+using StatsAPI: response, responsename, stderror, vcov, weights
+using StatsBase: StatsBase, CoefTable, model_response, summarystats
 using StatsFuns: log2π, normccdf
+using StatsModels: StatsModels, AbstractContrasts, AbstractTerm, CategoricalTerm
+using StatsModels: ConstantTerm, DummyCoding, EffectsCoding, FormulaTerm, FunctionTerm
+using StatsModels: HelmertCoding, HypothesisCoding, InteractionTerm, InterceptTerm
+using StatsModels: MatrixTerm, SeqDiffCoding, TableRegressionModel, Term
+using StatsModels: apply_schema, drop_term, formula, modelcols, term, @formula
+using StructTypes: StructTypes
+using Tables: Tables, columntable, rows
+using TypedTables: TypedTables, DictTable, FlexTable, Table
 
 export @formula,
     AbstractReMat,
@@ -77,6 +89,8 @@ export @formula,
     condVar,
     condVartables,
     confint,
+    dataset,
+    datasets,
     deviance,
     dispersion,
     dispersion_parameter,
@@ -190,8 +204,6 @@ include("mimeshow.jl")
 include("serialization.jl")
 include("profile/profile.jl")
 
-using PrecompileTools
-
 @setup_workload begin
     # Putting some things in `setup` can reduce the size of the
     # precompile file and potentially make loading faster.

src/Xymat.jl

@@ -48,11 +48,11 @@ end
 
 Convenience constructor for a sparse [`FeTerm`](@ref) assuming full rank, identity pivot and unit weights.
 
-Note: automatic rank deficiency handling may be added to this method in the future, as discused in
+Note: automatic rank deficiency handling may be added to this method in the future, as discussed in
 the vignette "[Rank deficiency in mixed-effects models](@ref)" for general `FeTerm`.
 """
 function FeTerm(X::SparseMatrixCSC, cnms::AbstractVector{String})
-    @debug "Full rank is assumed for sparse fixed-effect matrices."
+    #@debug "Full rank is assumed for sparse fixed-effect matrices."
     rank = size(X, 2)
     return FeTerm{eltype(X),typeof(X)}(X, collect(1:rank), rank, collect(cnms))
 end

src/bootstrap.jl

@@ -75,7 +75,7 @@ end
 
 Restore replicates from `f`, using `m` to create the desired subtype of [`MixedModelFitCollection`](@ref).
 
-`f` can be any entity suppored by `Arrow.Table`. `m` does not have to be fitted, but it must have
+`f` can be any entity supported by `Arrow.Table`. `m` does not have to be fitted, but it must have
 been constructed with the same structure as the source of the saved replicates.
 
 The two-argument method constructs a [`MixedModelBootstrap`](@ref) with the same eltype as `m`.
@@ -92,7 +92,7 @@ end
 # why this weird second method? it allows us to define custom types and write methods
 # to load into those types directly. For example, we could define a `PowerAnalysis <: MixedModelFitCollection`
 # in MixedModelsSim and then overload this method to get a convenient object. 
-# Also, this allows us to write `restorereplicateS(f, m, ::Type{<:MixedModelNonparametricBoostrap})` for
+# Also, this allows us to write `restorereplicateS(f, m, ::Type{<:MixedModelNonparametricBootstrap})` for
 # entities in MixedModels bootstrap
 function restorereplicates(
     f, m::MixedModel, ctype::Type{<:MixedModelFitCollection{T}}
@@ -193,7 +193,7 @@ performance benefits.
 - `β`, `σ`, and `θ` are the values of `m`'s parameters for simulating the responses.
 - `σ` is only valid for `LinearMixedModel` and `GeneralizedLinearMixedModel` for
 families with a dispersion parameter.
-- `progress` controls whehter the progress bar is shown. Note that the progress
+- `progress` controls whether the progress bar is shown. Note that the progress
 bar is automatically disabled for non-interactive (i.e. logging) contexts.
 - `optsum_overrides` is used to override values of [OptSummary](@ref) in the models
 fit during the bootstrapping process. For example, `optsum_overrides=(;ftol_rel=1e-08)`

src/grouping.jl

@@ -44,10 +44,10 @@ function StatsModels.apply_schema(
     context::AbstractTerm,
 )
     aliased = drop_term(context, t)
-    @debug "$t in context of $context: aliases $aliased\n  seen already: $(schema.already)"
+    #@debug "$t in context of $context: aliases $aliased\n  seen already: $(schema.already)"
     for seen in schema.already
         if StatsModels.symequal(aliased, seen)
-            @debug "  aliased term already present: $seen"
+            #@debug "  aliased term already present: $seen"
             return t
         end
     end
@@ -57,6 +57,6 @@ function StatsModels.apply_schema(
     # repair:
     new_contrasts = StatsModels.ContrastsMatrix(Grouping(), t.contrasts.levels)
     t = CategoricalTerm(t.sym, new_contrasts)
-    @debug "  aliased term absent, repairing: $t"
+    #@debug "  aliased term absent, repairing: $t"
     return t
 end

src/linalg.jl

@@ -18,7 +18,7 @@ function LinearAlgebra.mul!(
             αbnz = α * bnz[ib]
             jj = brv[ib]
             for ia in nzrange(A, j)
-                C[arv[ia], jj] += anz[ia] * αbnz
+                C[arv[ia], jj] = muladd(anz[ia], αbnz, C[arv[ia], jj])
             end
         end
     end

src/linalg/rankUpdate.jl

@@ -22,7 +22,7 @@ function MixedModels.rankUpdate!(
     Cdiag = C.data.diag
     Adiag = A.diag
     @inbounds for idx in eachindex(Cdiag, Adiag)
-        Cdiag[idx] = β * Cdiag[idx] + α * abs2(Adiag[idx])
+        Cdiag[idx] = muladd(β, Cdiag[idx], α * abs2(Adiag[idx]))
     end
     return C
 end
@@ -52,7 +52,7 @@ function _columndot(rv, nz, rngi, rngj)
     while i ≤ ni && j ≤ nj
         @inbounds ri, rj = rv[rngi[i]], rv[rngj[j]]
         if ri == rj
-            @inbounds accum += nz[rngi[i]] * nz[rngj[j]]
+            @inbounds accum = muladd(nz[rngi[i]], nz[rngj[j]], accum)
             i += 1
             j += 1
         elseif ri < rj
@@ -80,17 +80,17 @@ function rankUpdate!(C::HermOrSym{T,S}, A::SparseMatrixCSC{T}, α, β) where {T,
                 rvj = rv[j]
                 for i in k:lenrngjj
                     kk = rangejj[i]
-                    Cd[rv[kk], rvj] += nz[kk] * anzj
+                    Cd[rv[kk], rvj] = muladd(nz[kk], anzj, Cd[rv[kk], rvj])
                 end
             end
         end
     else
         @inbounds for j in axes(C, 2)
             rngj = nzrange(A, j)
             for i in 1:(j - 1)
-                Cd[i, j] += α * _columndot(rv, nz, nzrange(A, i), rngj)
+                Cd[i, j] = muladd(α, _columndot(rv, nz, nzrange(A, i), rngj), Cd[i, j])
             end
-            Cd[j, j] += α * sum(i -> abs2(nz[i]), rngj)
+            Cd[j, j] = muladd(α, sum(i -> abs2(nz[i]), rngj), Cd[j, j])
         end
     end
     return C
@@ -109,7 +109,7 @@ function rankUpdate!(
     isone(β) || rmul!(Cdiag, β)
 
     @inbounds for i in eachindex(Cdiag)
-        Cdiag[i] += α * sum(abs2, view(A, i, :))
+        Cdiag[i] = muladd(α, sum(abs2, view(A, i, :)), Cdiag[i])
     end
 
     return C
@@ -132,9 +132,9 @@ function rankUpdate!(
             AtAij = 0
             for idx in axes(A, 2)
                 # because the second multiplicant is from A', swap index order
-                AtAij += A[iind, idx] * A[jind, idx]
+                AtAij = muladd(A[iind, idx], A[jind, idx], AtAij)
             end
-            Cdat[i, j, k] += α * AtAij
+            Cdat[i, j, k] = muladd(α, AtAij, Cdat[i, j, k])
         end
     end
 
@@ -152,7 +152,7 @@ function rankUpdate!(
         throw(ArgumentError("Columns of A must have exactly 1 nonzero"))
 
     for (r, nz) in zip(rowvals(A), nonzeros(A))
-        dd[r] += α * abs2(nz)
+        dd[r] = muladd(α, abs2(nz), dd[r])
     end
 
     return C