From ecd0b3f07a3b962e7b36db0c4dd7279362dcc167 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Tue, 27 Jul 2021 19:49:58 -0500 Subject: [PATCH 1/3] remove dead path for pre julia 1.6 support --- src/linalg.jl | 138 ++++++++++++++++---------------------------------- 1 file changed, 43 insertions(+), 95 deletions(-) diff --git a/src/linalg.jl b/src/linalg.jl index db4e9146e..8e90849d9 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -45,104 +45,52 @@ 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 - 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.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::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::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::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 +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 From 3e2827bc8877d16f85d34aca6688aaf4f2ebc502 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Tue, 27 Jul 2021 20:08:52 -0500 Subject: [PATCH 2/3] pivoted qr on 1.7 without deprecation warnings --- src/likelihoodratiotest.jl | 2 +- src/linalg.jl | 6 ++++++ src/linalg/pivot.jl | 6 +++--- 3 files changed, 10 insertions(+), 4 deletions(-) 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 diff --git a/src/linalg.jl b/src/linalg.jl index 8e90849d9..5cf62bce9 100644 --- a/src/linalg.jl +++ b/src/linalg.jl @@ -94,3 +94,9 @@ function LinearAlgebra.rdiv!( 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 From 0afce6ccf7a8d6aa8e581d1aa949e02820c05aef Mon Sep 17 00:00:00 2001 From: Douglas Bates Date: Wed, 28 Jul 2021 10:19:49 -0500 Subject: [PATCH 3/3] Remove vestigial _zerocorr while code cleaning --- src/linearmixedmodel.jl | 42 ----------------------------------------- 1 file changed, 42 deletions(-) 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))