From b5a6a5501704eb9e42b007b46d5361541047eee2 Mon Sep 17 00:00:00 2001 From: willtebbutt Date: Sun, 21 Jul 2019 13:15:31 -0400 Subject: [PATCH 1/3] Refactor invlink + link --- src/Bijectors.jl | 30 ++++++++++-------------------- 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/Bijectors.jl b/src/Bijectors.jl index 6d7c6f6e..79dd0445 100644 --- a/src/Bijectors.jl +++ b/src/Bijectors.jl @@ -323,36 +323,26 @@ end const PDMatDistribution = Union{InverseWishart, Wishart} -function link(d::PDMatDistribution, X::AbstractMatrix{T}) where {T<:Real} - Y = cholesky(X).L - @inbounds @simd for m in 1:size(Y, 1) - Y[m, m] = log(Y[m, m]) - end - return Matrix(Y) +function link(d::PDMatDistribution, X::AbstractMatrix{<:Real}) + Y = Matrix(cholesky(X).L) + Y[diagind(Y)] .= log.(view(Y, diagind(Y))) + return Y end -function invlink(d::PDMatDistribution, Y::AbstractMatrix{T}) where {T<:Real} - X, dim = copy(Y), size(Y) - @inbounds @simd for m in 1:size(X, 1) - X[m, m] = exp(X[m, m]) - end - Z = similar(X) - return mul!(Z, LowerTriangular(X), LowerTriangular(X)') +function invlink(d::PDMatDistribution, Y::AbstractMatrix{<:Real}) + X = copy(Y) + X[diagind(X)] .= exp.(view(X, diagind(X))) + return LowerTriangular(X) * LowerTriangular(X)' end -function logpdf_with_trans( - d::PDMatDistribution, - X::AbstractMatrix{<:Real}, - transform::Bool -) - T = eltype(X) +function logpdf_with_trans(d::PDMatDistribution, X::AbstractMatrix{<:Real}, transform::Bool) lp = logpdf(d, X) if transform && isfinite(lp) U = cholesky(X).U @inbounds @simd for i in 1:dim(d) lp += (dim(d) - i + 2) * log(U[i, i]) end - lp += dim(d) * log(T(2)) + lp += dim(d) * log(eltype(X)(2)) end return lp end From 3224ee5a9a037079adea3cda77a4219465e8416d Mon Sep 17 00:00:00 2001 From: willtebbutt Date: Sun, 21 Jul 2019 13:39:50 -0400 Subject: [PATCH 2/3] Tidy up logpdf_with_trans --- src/Bijectors.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/Bijectors.jl b/src/Bijectors.jl index 79dd0445..7d2b1ad1 100644 --- a/src/Bijectors.jl +++ b/src/Bijectors.jl @@ -339,10 +339,8 @@ function logpdf_with_trans(d::PDMatDistribution, X::AbstractMatrix{<:Real}, tran lp = logpdf(d, X) if transform && isfinite(lp) U = cholesky(X).U - @inbounds @simd for i in 1:dim(d) - lp += (dim(d) - i + 2) * log(U[i, i]) - end - lp += dim(d) * log(eltype(X)(2)) + lp += sum(i->(dim(d) - i + 2) * log(U[i, i]), 1:dim(d)) + lp += dim(d) * log(2) end return lp end From 74ba08cf036cde5764ea63f60a28faa1c2b08832 Mon Sep 17 00:00:00 2001 From: willtebbutt Date: Sun, 21 Jul 2019 14:00:30 -0400 Subject: [PATCH 3/3] Minor tweak --- src/Bijectors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Bijectors.jl b/src/Bijectors.jl index 7d2b1ad1..d1d47e43 100644 --- a/src/Bijectors.jl +++ b/src/Bijectors.jl @@ -339,7 +339,7 @@ function logpdf_with_trans(d::PDMatDistribution, X::AbstractMatrix{<:Real}, tran lp = logpdf(d, X) if transform && isfinite(lp) U = cholesky(X).U - lp += sum(i->(dim(d) - i + 2) * log(U[i, i]), 1:dim(d)) + lp += sum((dim(d) .- (1:dim(d)) .+ 2) .* log.(view(U, diagind(U)))) lp += dim(d) * log(2) end return lp