Skip to content

Commit

Permalink
compat fixes for Julia 1.7 (#547)
Browse files Browse the repository at this point in the history
* remove dead path for pre julia 1.6 support

* pivoted qr on 1.7 without deprecation warnings

* Remove vestigial _zerocorr while code cleaning

Co-authored-by: Douglas Bates <dmbates@gmail.com>
  • Loading branch information
palday and dmbates authored Jul 28, 2021
1 parent 1edd541 commit fcb4357
Show file tree
Hide file tree
Showing 4 changed files with 52 additions and 140 deletions.
2 changes: 1 addition & 1 deletion src/likelihoodratiotest.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
142 changes: 48 additions & 94 deletions src/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
6 changes: 3 additions & 3 deletions src/linalg/pivot.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))`.
Expand All @@ -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
Expand All @@ -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
Expand Down
42 changes: 0 additions & 42 deletions src/linearmixedmodel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))

0 comments on commit fcb4357

Please sign in to comment.