Skip to content

Commit

Permalink
Merge pull request #24 from PharmCat/dev14
Browse files Browse the repository at this point in the history 
internal structure changes
  • Loading branch information
@PharmCat
PharmCat authored Dec 13, 2022
2 parents f7954ec + 059a464 commit 34c5c98
Show file tree
Hide file tree
Showing 9 changed files with 99 additions and 46 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
DiffResults = "1"
ForwardDiff = "0.10"
LineSearches = "7"
MetidaBase = "0.10.1"
MetidaBase = "0.10.1, 0.10.2"
Optim = "1"
ProgressMeter = "1"
StatsBase = "0.29, 0.30, 0.31, 0.32, 0.33"
Expand Down
9 changes: 5 additions & 4 deletions change.log
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,11 @@
V0.14.2
* responce name

v0.14.1
* docs
* test
* minor changes
* minor optimization
* internal structure changes
* responce name
* raneff initial realization
* bootstrap changes (output, two-stage botstrap, api)

v0.14
* lmmformula
Expand Down
7 changes: 4 additions & 3 deletions src/dof_contain.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,10 +5,11 @@ function zmatrix(lmm::LMM, i)
sn = 1
for b = 1:length(lmm.covstr.vcovblock)
zblock = view(lmm.covstr.z, lmm.covstr.vcovblock[b], lmm.covstr.zrndur[i])
for s = 1:length(lmm.covstr.sblock[b][i])
zi = view(zblock, lmm.covstr.sblock[b][i][s], :)
for s = 1:subjn(lmm.covstr, i, b)
suji = getsubj(lmm.covstr, i, b, s)
zi = view(zblock, suji, :)
l[sn] = 1
copyto!(view(rzm, lmm.covstr.vcovblock[b][lmm.covstr.sblock[b][i][s]], :), kron(l, zi))
copyto!(view(rzm, lmm.covstr.vcovblock[b][suji], :), kron(l, zi))
l[sn] = 0
sn += 1
end
Expand Down
23 changes: 13 additions & 10 deletions src/gmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,37 +3,37 @@
################################################################################

################################################################################
function gmatvec(θ::Vector{T}, covstr) where T
gt = [Symmetric(zeros(T, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn] #<: fix to vec
@noinline function gmatvec(θ::Vector{T}, covstr) where T
gt = [Symmetric(zeros(T, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn]
for r = 1:covstr.rn
#fill!(gt[r].data, zero(T))
if covstr.random[r].covtype.z
gmat!(gt[r].data, view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
end
end
gt
end
# Main
function zgz_base_inc!(mx::AbstractArray, G, θ, covstr, block, sblock)
@noinline function zgz_base_inc!(mx::AbstractArray, G, covstr, bi)
block = covstr.vcovblock[bi]
if covstr.random[1].covtype.z
for r = 1:covstr.rn
zblock = view(covstr.z, block, covstr.zrndur[r])
@inbounds for i = 1:length(sblock[r])
mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i]), view(zblock, sblock[r][i], :), G[r])
@inbounds for i = 1:subjn(covstr, r, bi)
sb = getsubj(covstr, r, bi, i)
mulαβαtinc!(view(mx, sb, sb), view(zblock, sb, :), G[r])
end
end
end
mx
end


################################################################################

function zgz_base_inc!(mx::AbstractArray, θ::AbstractArray{T}, covstr, block, sblock) where T
#=
@noinline function zgz_base_inc!(mx::AbstractArray, θ::AbstractArray{T}, covstr, block, sblock) where T
if covstr.random[1].covtype.z
for r = 1:covstr.rn
G = fill!(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])), zero(T))
gmat_switch!(G.data, θ, covstr, r)
gmat!(G.data, view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
zblock = view(covstr.z, block, covstr.zrndur[r])
@inbounds for i = 1:length(sblock[r])
mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i]), view(zblock, sblock[r][i], :), G)
Expand All @@ -42,10 +42,13 @@ function zgz_base_inc!(mx::AbstractArray, θ::AbstractArray{T}, covstr, block, s
end
mx
end
=#
#=
function gmat_switch!(G, θ, covstr, r)
gmat!(G, view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
G
end
=#
################################################################################
#SI
function gmat!(::Any, ::Any, ::AbstractCovarianceType)
Expand Down
6 changes: 4 additions & 2 deletions src/random.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,10 +27,11 @@ function rand!(rng::AbstractRNG, v::AbstractVector, lmm::LMM{T}, theta::Abstract
v = Vector{T}(undef, nobs(lmm))
tv = Vector{T}(undef, lmm.maxvcbl)
gvec = gmatvec(theta, lmm.covstr)
rtheta = theta[lmm.covstr.tr[end]]
for i = 1:n
q = length(lmm.covstr.vcovblock[i])
V = zeros(q, q)
vmatrix!(V, gvec, theta, lmm, i)
vmatrix!(V, gvec, rtheta, lmm, i)
if length(tv) != q resize!(tv, q) end
Distributions.rand!(rng, MvNormal(Symmetric(V)), tv)
v[lmm.covstr.vcovblock[i]] .= tv
Expand All @@ -54,6 +55,7 @@ function rand!(rng::AbstractRNG, v::AbstractVector, lmm::LMM{T}, theta::Abstract
tv = Vector{T}(undef, lmm.maxvcbl)
m = Vector{T}(undef, lmm.maxvcbl)
gvec = gmatvec(theta, lmm.covstr)
rtheta = theta[lmm.covstr.tr[end]]
for i = 1:n
q = length(lmm.covstr.vcovblock[i])
if length(tv) != q
Expand All @@ -62,7 +64,7 @@ function rand!(rng::AbstractRNG, v::AbstractVector, lmm::LMM{T}, theta::Abstract
end
mul!(m, lmm.dv.xv[i], beta)
V = zeros(q, q)
vmatrix!(V, gvec, theta, lmm, i)
vmatrix!(V, gvec, rtheta, lmm, i)
Distributions.rand!(rng, MvNormal(m, Symmetric(V)), tv)
v[lmm.covstr.vcovblock[i]] .= tv
end
Expand Down
8 changes: 4 additions & 4 deletions src/reml.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T
#---------------------------------------------------------------------------
#logdetθ₂ = zero(T)
gvec = gmatvec(θ, lmm.covstr)
rθ = θ[lmm.covstr.tr[end]] # R part of θ
noerror = true
ncore = min(num_cores(), n, 16)
accθ₁ = zeros(T, ncore)
Expand Down Expand Up @@ -80,7 +81,7 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T
fillzeroutri!(Vc)
#-------------------------------------------------------------------
# Make V matrix
vmatrix!(V, gvec, θ, lmm, i)
vmatrix!(V, gvec, rθ, lmm, i)
#-----------------------------------------------------------------------
if length(swtw[t]) != qswm resize!(swtw[t], qswm) end
swm, swr, ne = sweepb!(swtw[t], Vp, 1:q; logdet = true, syrkblas = syrkblas)
Expand Down Expand Up @@ -132,10 +133,9 @@ function core_sweep_β(lmm, data, θ::Vector{T}, β, n) where T
accθ₃ = zeros(T, ncore)
erroracc = trues(ncore)
gvec = gmatvec(θ, lmm.covstr)
rθ = θ[lmm.covstr.tr[end]] # R part of θ
d, r = divrem(n, ncore)
Base.Threads.@threads for t = 1:ncore
#@batch for t = 1:ncore
#for t = 1:ncore
offset = min(t-1, r) + (t-1)*d
accθ₂[t] = zeros(T, lmm.rankx, lmm.rankx)
@inbounds for j ∈ 1:d+(t ≤ r)
Expand All @@ -149,7 +149,7 @@ function core_sweep_β(lmm, data, θ::Vector{T}, β, n) where T
fillzeroutri!(V)
copyto!(Vx, data.xv[i])
fillzeroutri!(Vc)
vmatrix!(V, gvec, θ, lmm, i)
vmatrix!(V, gvec, rθ, lmm, i)
#-----------------------------------------------------------------------
swm, swr, ne = sweepb!(Vector{T}(undef, qswm), Vp, 1:q; logdet = true)
#-----------------------------------------------------------------------
Expand Down
10 changes: 10 additions & 0 deletions src/rmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,16 @@ function rmat_base_inc!(mx, θ, covstr, block, sblock)
end
mx
end
@noinline function rmat_base_inc!(mx, θ, covstr, bi)
en = covstr.rn + 1
block = covstr.vcovblock[bi]
zblock = view(covstr.rz, block, :)
@simd for i = 1:subjn(covstr, en, bi)
sb = getsubj(covstr, en, bi, i)
rmat!(view(mx, sb, sb), θ, view(zblock, sb, :), covstr.repeated.covtype.s)
end
mx
end
################################################################################
function rmat!(::Any, ::Any, ::Any, ::AbstractCovarianceType)
error("No rmat! method defined for thit structure!")
Expand Down
29 changes: 18 additions & 11 deletions src/utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,7 +243,7 @@ function rmatrix(lmm::LMM{T}, i::Int) where T
if i > length(lmm.covstr.vcovblock) error("Invalid block number: $(i)!") end
q = length(lmm.covstr.vcovblock[i])
R = zeros(T, q, q)
rmat_base_inc!(R, lmm.result.theta[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
rmat_base_inc!(R, lmm.result.theta[lmm.covstr.tr[end]], lmm.covstr, i)
Symmetric(R)
end
"""
Expand All @@ -253,16 +253,21 @@ Update variance-covariance matrix V for i bolock. Upper triangular updated.
"""
function vmatrix!(V, θ, lmm::LMM, i::Int) # pub API
gvec = gmatvec(θ, lmm.covstr)
zgz_base_inc!(V, gvec, θ, lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
zgz_base_inc!(V, gvec, lmm.covstr, i)
rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, i)
end

# !!! Main function REML used
function vmatrix!(V, G, rθ, lmm::LMM, i::Int)
zgz_base_inc!(V, G, lmm.covstr, i)
rmat_base_inc!(V, rθ, lmm.covstr, i)
end
#=
function vmatrix!(V, G, θ, lmm::LMM, i::Int)
zgz_base_inc!(V, G, θ, lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
end

=#
"""
vmatrix(lmm::LMM, i::Int)
Expand All @@ -275,18 +280,18 @@ end
function vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int) where T
V = zeros(T, length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i]))
gvec = gmatvec(θ, lmm.covstr)
vmatrix!(V, gvec, θ, lmm, i)
vmatrix!(V, gvec, θ[lmm.covstr.tr[end]], lmm, i)
Symmetric(V)
end

#deprecated
# For Multiple Imputation
function vmatrix(θ::Vector, covstr::CovStructure, i::Int)
V = zeros(length(covstr.vcovblock[i]), length(covstr.vcovblock[i]))
gvec = gmatvec(θ, covstr)
zgz_base_inc!(V, gvec, θ, covstr, covstr.vcovblock[i], covstr.sblock[i])
rmat_base_inc!(V, θ[covstr.tr[end]], covstr, covstr.vcovblock[i], covstr.sblock[i])
zgz_base_inc!(V, gvec, covstr, i)
rmat_base_inc!(V, θ[covstr.tr[end]], covstr, i)
Symmetric(V)
end

#=
function grad_vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int)
V = zeros(T, length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i]))
Expand Down Expand Up @@ -319,7 +324,8 @@ function blockgmatrix(lmm::LMM{T}, v) where T
end

function blockzmatrix(lmm::LMM{T}, i) where T
sn = length.(lmm.covstr.sblock[i])[1:end-1]
#sn = length.(lmm.covstr.sblock[i])[1:end-1]
sn = raneflenv(lmm.covstr, i)
mx = Vector{Matrix{T}}(undef, raneffn(lmm))
s = 1
for j = 1:raneffn(lmm)
Expand All @@ -337,7 +343,8 @@ Vector of random effect coefficients for block `i`.
"""
function raneff(lmm::LMM{T}, i) where T
if raneffn(lmm) == 0 return nothing end
sn = length.(lmm.covstr.sblock[i])[1:end-1]
sn = raneflenv(lmm.covstr, i)
#sn = length.(lmm.covstr.esb.sblock[i, :])[1:end-1]
G = blockgmatrix(lmm, sn)
Z = blockzmatrix(lmm, i)
G * Z' * inv(vmatrix(lmm, i)) * (lmm.dv.yv[i] - lmm.dv.xv[i]*lmm.result.beta)
Expand Down
51 changes: 40 additions & 11 deletions src/varstruct.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,27 @@ end

tabcols(data, symbs) = Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs)


struct EffectSubjectBlock
sblock::Matrix{Vector{Tuple{Vector{Int}, Int}}}
snames::Vector
end
function getsubj(covstr, effn, block, sbjn)
covstr.esb.sblock[block, effn][sbjn][1]
end
function subjn(covstr, effn, block)
length(covstr.esb.sblock[block, effn])
end
function raneflenv(covstr, block)
l = size(covstr.esb.sblock, 2) - 1
v = Vector{Int}(undef, l)
for i = 1:l
v[i] = length(covstr.esb.sblock[block, i])
end
v
end
"""
Covarince structure.
"""
struct CovStructure{T} <: AbstractCovarianceStructure
# Random effects
random::Vector{VarEffect}
Expand All @@ -159,7 +179,9 @@ struct CovStructure{T} <: AbstractCovarianceStructure
z::Matrix{T}
#subjz::Vector{BitArray{2}}
# Blocks for each blocking subject, each effect, each effect subject sblock[block][rand eff][subj]
sblock::Vector{Vector{Vector{Vector{Int}}}}
#sblock::Vector{Vector{Vector{Vector{Int}}}}
#
esb::EffectSubjectBlock
#unit range z column range for each random effect
zrndur::Vector{UnitRange{Int}}
# repeated effect parametrization matrix
Expand Down Expand Up @@ -276,18 +298,23 @@ struct CovStructure{T} <: AbstractCovarianceStructure
end
blocks = collect(values(subjblockdict))
#end

sblock = Vector{Vector{Vector{Vector{Int}}}}(undef, length(blocks))
########################################################################
sblock = Matrix{Vector{Tuple{Vector{Int}, Int}}}(undef, length(blocks), alleffl)
nblock = []
#######################################################################
#######################################################################
nli = 1
@inbounds for i = 1:length(blocks) # i - block number
sblock[i] = Vector{Vector{Vector{Int}}}(undef, alleffl)
@inbounds for s = 1:alleffl # s - effect number
sblock[i][s] = Vector{Vector{Int}}(undef, 0)

tempv = Vector{Tuple{Vector{Int}, Int}}(undef, 0)
for (k, v) in dicts[s]
fa = findall(x-> x in v, blocks[i]) # Try to optimize it
if length(fa) > 0 push!(sblock[i][s], fa) end
if length(fa) > 0
push!(tempv, (fa, nli))
push!(nblock, k)
nli += 1
end
end
sblock[i, s] = tempv
end
end
#
Expand All @@ -296,10 +323,12 @@ struct CovStructure{T} <: AbstractCovarianceStructure
lvcb = length(i)
if lvcb > maxn maxn = lvcb end
end
new{eltype(z)}(random, repeated, schema, rcnames, blocks, rn, z, sblock, zrndur, rz, q, t, tr, tl, ct, sn, maxn)
esb = EffectSubjectBlock(sblock, nblock)
#######################################################################
new{eltype(z)}(random, repeated, schema, rcnames, blocks, rn, z, esb, zrndur, rz, q, t, tr, tl, ct, sn, maxn)
end
end
################################################################################
###############################################################################
function fillur!(ur, i, v)
if i > 1
ur[i] = UnitRange(sum(v[1:i-1]) + 1, sum(v[1:i-1]) + v[i])
Expand Down

2 comments on commit 34c5c98

@PharmCat
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/74055

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.14.1 -m "<description of version>" 34c5c9883fa3060401a70b9ecf589a8a0f1630eb
git push origin v0.14.1

Please sign in to comment.