Repository: JuliaStats/MixedModels.jl (julia)

Issue:
Adjustments for julia-v1.7.0
When the workaround for https://github.com/JuliaLang/julia/issues/40048 is merged we can start testing on julia-v1.7.0-DEV

In `test/runtests.jl` we should replace the call to determine the BLAS in use with `BLAS.get_config()` for v1.7.0-DEV

A couple of instances of `qr(x, Val(true))` need to be replaced by `qr(x, ColumnNorm())` for recent versions of v1.7.0-DEV.

Pull request:
Closes #532 

- deleted a dead path now that we've dropped pre 1.6 support
- added an internal wrapper method to dodge some deprecation warnings in 1.7+

Patch:
diff --git a/src/linalg.jl b/src/linalg.jl
index db4e9146e..5cf62bce9 100644
--- a/src/linalg.jl
+++ b/src/linalg.jl
@@ -45,104 +45,58 @@ function LinearAlgebra.mul!(
     return mul!(C, adjoint(adjA.parent.cscmat), B, α, β)
 end
 
-@static if VERSION < v"1.6.0-DEV.1468"
-    function LinearAlgebra.ldiv!(
-        adjA::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}}, B::StridedVector{T}
-    ) where {T}
-        A = adjA.parent
-        length(B) == size(A, 2) || throw(DimensionMismatch(""))
-        Adat = A.data.data
-        m, n, k = size(Adat)
-        bb = reshape(B, (n, k))
-        for j in axes(Adat, 3)
-            ldiv!(adjoint(LowerTriangular(view(Adat, :, :, j))), view(bb, :, j))
-        end
-        return B
-    end
-
-    function LinearAlgebra.rdiv!(
-        A::Matrix{T}, adjB::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}}
-    ) where {T}
-        m, n = size(A)
-        Bd = adjB.parent.data
-        Bdd = Bd.data
-        r, s, blk = size(Bdd)
-        n == size(Bd, 1) && r == s || throw(DimensionMismatch())
-        for b in axes(Bd.data, 3)
-            coloffset = (b - 1) * s
-            rdiv!(
-                view(A, :, (coloffset + 1):(coloffset + s)),
-                adjoint(LowerTriangular(view(Bdd, :, :, b))),
-            )
-        end
-        return A
+function LinearAlgebra.ldiv!(
+    A::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}, B::StridedVector{T}
+) where {T}
+    adjA = A.data
+    length(B) == size(A, 2) || throw(DimensionMismatch(""))
+    Adat = adjA.parent.data
+    m, n, k = size(Adat)
+    bb = reshape(B, (n, k))
+    for j in axes(Adat, 3)
+        ldiv!(UpperTriangular(adjoint(view(Adat, :, :, j))), view(bb, :, j))
     end
+    return B
+end
 
-    function LinearAlgebra.rdiv!(
-        A::BlockedSparse{T,S,P}, B::Adjoint{T,<:LowerTriangular{T,UniformBlockDiagonal{T}}}
-    ) where {T,S,P}
-        Bpd = B.parent.data
-        Bdat = Bpd.data
-        j, k, l = size(Bdat)
-        cbpt = A.colblkptr
-        nzv = A.cscmat.nzval
-        P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
-        for j in axes(Bdat, 3)
-            rdiv!(
-                reshape(view(nzv, cbpt[j]:(cbpt[j + 1] - 1)), :, P),
-                adjoint(LowerTriangular(view(Bdat, :, :, j))),
-            )
-        end
-        return A
-    end
-else
-    function LinearAlgebra.ldiv!(
-        A::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}, B::StridedVector{T}
-    ) where {T}
-        adjA = A.data
-        length(B) == size(A, 2) || throw(DimensionMismatch(""))
-        Adat = adjA.parent.data
-        m, n, k = size(Adat)
-        bb = reshape(B, (n, k))
-        for j in axes(Adat, 3)
-            ldiv!(UpperTriangular(adjoint(view(Adat, :, :, j))), view(bb, :, j))
-        end
-        return B
+function LinearAlgebra.rdiv!(
+    A::Matrix{T}, B::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}
+) where {T}
+    m, n = size(A)
+    Bd = B.data.parent
+    Bdd = Bd.data
+    r, s, blk = size(Bdd)
+    n == size(Bd, 1) && r == s || throw(DimensionMismatch())
+    for b in axes(Bd.data, 3)
+        coloffset = (b - 1) * s
+        rdiv!(
+            view(A, :, (coloffset + 1):(coloffset + s)),
+            UpperTriangular(adjoint(view(Bdd, :, :, b))),
+        )
     end
+    return A
+end
 
-    function LinearAlgebra.rdiv!(
-        A::Matrix{T}, B::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}
-    ) where {T}
-        m, n = size(A)
-        Bd = B.data.parent
-        Bdd = Bd.data
-        r, s, blk = size(Bdd)
-        n == size(Bd, 1) && r == s || throw(DimensionMismatch())
-        for b in axes(Bd.data, 3)
-            coloffset = (b - 1) * s
-            rdiv!(
-                view(A, :, (coloffset + 1):(coloffset + s)),
-                UpperTriangular(adjoint(view(Bdd, :, :, b))),
-            )
-        end
-        return A
+function LinearAlgebra.rdiv!(
+    A::BlockedSparse{T,S,P}, B::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}
+) where {T,S,P}
+    Bpd = B.data.parent
+    Bdat = Bpd.data
+    j, k, l = size(Bdat)
+    cbpt = A.colblkptr
+    nzv = A.cscmat.nzval
+    P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
+    for j in axes(Bdat, 3)
+        rdiv!(
+            reshape(view(nzv, cbpt[j]:(cbpt[j + 1] - 1)), :, P),
+            UpperTriangular(adjoint(view(Bdat, :, :, j))),
+        )
     end
+    return A
+end
 
-    function LinearAlgebra.rdiv!(
-        A::BlockedSparse{T,S,P}, B::UpperTriangular{T,<:Adjoint{T,UniformBlockDiagonal{T}}}
-    ) where {T,S,P}
-        Bpd = B.data.parent
-        Bdat = Bpd.data
-        j, k, l = size(Bdat)
-        cbpt = A.colblkptr
-        nzv = A.cscmat.nzval
-        P == j == k && length(cbpt) == (l + 1) || throw(DimensionMismatch(""))
-        for j in axes(Bdat, 3)
-            rdiv!(
-                reshape(view(nzv, cbpt[j]:(cbpt[j + 1] - 1)), :, P),
-                UpperTriangular(adjoint(view(Bdat, :, :, j))),
-            )
-        end
-        return A
-    end
+@static if VERSION < v"1.7.0-DEV.1188" # julialang sha e0ecc557a24eb3338b8dc672d02c98e8b31111fa
+    pivoted_qr(A; kwargs...) = qr(A, Val(true); kwargs...)
+else
+    pivoted_qr(A; kwargs...) = qr(A, ColumnNorm(); kwargs...)
 end
diff --git a/src/linalg/pivot.jl b/src/linalg/pivot.jl
index d0195a088..bf1f70e16 100644
--- a/src/linalg/pivot.jl
+++ b/src/linalg/pivot.jl
@@ -4,7 +4,7 @@
 Return the numerical column rank and a pivot vector.
 
 The rank is determined from the absolute values of the diagonal of R from
-a pivoted QR decomposition, relative to the first (and, hence, largest) 
+a pivoted QR decomposition, relative to the first (and, hence, largest)
 element of this vector.
 
 In the full-rank case the pivot vector is `collect(axes(x, 2))`.
@@ -15,7 +15,7 @@ function statsrank(x::AbstractMatrix{T}; ranktol=1e-8) where {T<:AbstractFloat}
 
     iszero(n) && return (rank=n, piv=piv)
 
-    qrpiv = qr(x, Val(true))
+    qrpiv = pivoted_qr(x)
     dvec = abs.(diag(qrpiv.R))
     fdv = first(dvec)
     cmp = fdv * ranktol
@@ -28,7 +28,7 @@ function statsrank(x::AbstractMatrix{T}; ranktol=1e-8) where {T<:AbstractFloat}
     if all(isone, v1) && first(piv) ≠ 1
         # make sure the first column isn't moved by inflating v1
         v1 .*= (fdv + one(fdv)) / sqrt(m)
-        qrpiv = qr(x, Val(true))
+        qrpiv = pivoted_qr(x)
         piv = qrpiv.p
         fill!(v1, one(T))    # restore the contents of the first column
     end
diff --git a/src/linearmixedmodel.jl b/src/linearmixedmodel.jl
index 097386221..abdc1e7de 100644
--- a/src/linearmixedmodel.jl
+++ b/src/linearmixedmodel.jl
@@ -1177,45 +1177,3 @@ function StatsBase.weights(m::LinearMixedModel)
     rtwts = m.sqrtwts
     return isempty(rtwts) ? ones(eltype(rtwts), nobs(m)) : abs2.(rtwts)
 end
-
-"""
-    _zerocorr!(m::LinearMixedModel[, trmnms::Vector{Symbol}])
-
-Rewrite the random effects specification for the grouping factors in `trmnms` to zero correlation parameter.
-
-The default for `trmnms` is all the names of random-effects terms.
-
-A random effects term is in the zero correlation parameter configuration when the off-diagonal elements of
-λ are all zero - hence there are no correlation parameters in that term being estimated.
-
-Note that this is numerically equivalent to specifying a formula with `zerocorr` around each random effects
-term, but the `formula`  fields in the resulting model will differ. In particular, `zerocorr!` will **not**
-change the original `formula`'s terms to be of type of `ZeroCorr` because this would involve changing
-immutable types.  This may have implications for software that manipulates the formula of a fitted model.
-
-This is an internal function and may disappear in a future release without being considered a breaking change.
-"""
-function _zerocorr!(m::LinearMixedModel{T}, trmns) where {T}
-    reterms = m.reterms
-    for trm in reterms
-        if fname(trm) in trmns
-            zerocorr!(trm)
-        end
-    end
-    newparmap = mkparmap(reterms)
-    copyto!(m.parmap, newparmap)
-    resize!(m.parmap, length(newparmap))
-    optsum = m.optsum
-    optsum.lowerbd = mapfoldl(lowerbd, vcat, reterms)
-    optsum.initial = mapfoldl(getθ, vcat, reterms)
-    optsum.final = copy(optsum.initial)
-    optsum.xtol_abs = fill!(copy(optsum.initial), 1.0e-10)
-    optsum.initial_step = T[]
-
-    # the model is no longer fitted
-    optsum.feval = -1
-
-    return m
-end
-
-_zerocorr!(m::LinearMixedModel) = _zerocorr!(m, fnames(m))


Test patch:
diff --git a/src/likelihoodratiotest.jl b/src/likelihoodratiotest.jl
index 5023d64e9..ce4c1a181 100644
--- a/src/likelihoodratiotest.jl
+++ b/src/likelihoodratiotest.jl
@@ -353,7 +353,7 @@ function _isnested(x::AbstractMatrix, y::AbstractMatrix; rtol=1e-8, ranktol=1e-8
 
     qy = qr(y).Q
 
-    qrx = qr(x, Val(true))
+    qrx = pivoted_qr(x)
     dvec = abs.(diag(qrx.R))
     fdv = first(dvec)
     cmp = fdv * ranktol