Skip to content

Commit

Permalink
Merge pull request #28 from PharmCat/dev14
Browse files Browse the repository at this point in the history 
clean code
  • Loading branch information
@PharmCat
PharmCat authored Dec 28, 2022
2 parents b3bd87a + a0cc527 commit 959dc81
Show file tree
Hide file tree
Showing 21 changed files with 741 additions and 935 deletions.
485 changes: 485 additions & 0 deletions src/deprecated.jl
View file

Large diffs are not rendered by default.

12 changes: 0 additions & 12 deletions src/dof_satter.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@ function gradc(lmm::LMM{T}, theta) where T
if !isnothing(lmm.result.grc) return lmm.result.grc end
vloptf(x) = sweep_β_cov(lmm, lmm.dv, x, lmm.result.beta)
chunk = ForwardDiff.Chunk{min(10, length(theta))}()
#chunk = ForwardDiff.Chunk{1}()
jcfg = ForwardDiff.JacobianConfig(vloptf, theta, chunk)
jic = ForwardDiff.jacobian(vloptf, theta, jcfg)
grad = Vector{Matrix{T}}(undef, thetalength(lmm))
Expand Down Expand Up @@ -74,15 +73,11 @@ function dof_satter(lmm::LMM{T}, l::AbstractVector) where T
grad = gradc(lmm, theta)
g = Vector{T}(undef, length(grad))
for i = 1:length(grad)
#g[i] = (l' * grad[i] * l)[1]
#g[i] = mulαtβα(l, grad[i])
g[i] = dot(l, grad[i], l)
end
#d = g' * A * g
d = dot(g, A, g)
#d = mulαtβα(g, A)
df = 2*(dot(l, lmm.result.c, l))^2 / d
#df = 2*(mulαtβα(l, lmm.result.c))^2 / d
if df < 1.0 return 1.0 elseif df > dof_residual(lmm) return dof_residual(lmm) else return df end
end
"""
Expand Down Expand Up @@ -112,14 +107,10 @@ function dof_satter(lmm::LMM{T}) where T
l[gi] = 1
g = Vector{T}(undef, length(grad))
for i = 1:length(grad)
#g[i] = (l' * grad[i] * l)[1]
#g[i] = mulαtβα(l, grad[i])
g[i] = dot(l, grad[i], l)
end
#d = g' * A * g
#d = mulαtβα(g, A)
d = dot(g, A, g)
#df = 2*(mulαtβα(l, lmm.result.c))^2 / d
df = 2*(dot(l, lmm.result.c, l))^2 / d
if df < 1.0 dof[gi] = 1.0 elseif df > dof_residual(lmm) dof[gi] = dof_residual(lmm) else dof[gi] = df end
end
Expand Down Expand Up @@ -160,12 +151,9 @@ function dof_satter(lmm::LMM{T}, l::AbstractMatrix) where T
for i = 1:lclr
plm = pl[i,:]
for i2 = 1:length(grad)
#g[i2] = (plm' * grad[i2] * plm)[1]
#g[i2] = mulαtβα(plm, grad[i2])
g[i2] = dot(plm, grad[i2], plm)
end
#d = g' * A * g
#d = mulαtβα(g, A)
d = dot(g, A, g)
vm[i] = 2*lcle.values[i]^2 / d
if vm[i] > 2.0 em += vm[i] / (vm[i] - 2.0) end
Expand Down
2 changes: 0 additions & 2 deletions src/estimate.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ Estimate table for l vector. Satter DF used.
"""
function estimate(lmm, l::AbstractVector; level = 0.95, name = "Estimate")
est = coef(lmm)'*l
#se = sqrt(mulαtβα(l, vcov(lmm)))
se = sqrt(dot(l, vcov(lmm), l))
df = dof_satter(lmm, l)
t = abs(est/se)
Expand Down Expand Up @@ -47,7 +46,6 @@ function estimate(lmm; level = 0.95)
fill!(l, 0)
l[i] = 1
vest[i] = coe'*l
#vse[i] = sqrt(mulαtβα(l, vcov(lmm)))
vse[i] = sqrt(dot(l, vcov(lmm), l))
vdf[i] = dof_satter(lmm, l)
vt[i] = abs(vest[i]/vse[i])
Expand Down
4 changes: 0 additions & 4 deletions src/fit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,6 @@ function fit!(lmm::LMM{T}; kwargs...) where T
end
else
initθ = sqrt(initvar(lmm.data.yv, lmm.mm.m)[1])/(length(lmm.covstr.random)+1)
#θ .= initθ
for i = 1:length(θ)
if lmm.covstr.ct[i] == :var
θ[i] = initθ
Expand All @@ -147,7 +146,6 @@ function fit!(lmm::LMM{T}; kwargs...) where T
end
# Initial step with modified Newton method
chunk = ForwardDiff.Chunk{min(8, length(θ))}()
#chunk = ForwardDiff.Chunk{1}()
if isa(aifirst, Bool)
if aifirst aifirst == :ai else aifirst == :default end
end
Expand Down Expand Up @@ -203,7 +201,6 @@ function fit!(lmm::LMM{T}; kwargs...) where T
# SE
if !any(x -> x < 0.0, diag(lmm.result.c))
diag!(sqrt, lmm.result.se, lmm.result.c)
#lmm.result.se = sqrt.(diag(lmm.result.c))
if any(x-> x < singtol, lmm.result.se) && minimum(lmm.result.se)/maximum(lmm.result.se) < singtol lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Some of the SE parameters is suspicious.")) end
lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Model fitted."))
lmm.result.fit = true
Expand Down Expand Up @@ -252,7 +249,6 @@ function optstep!(lmm, data, θ; method::Symbol = :ai, maxopt::Int=10)
reml, beta, θs₂, θ₃, rt = reml_sweep_β(lmm, data, θ)
rt || error("Wrong initial conditions.")
chunk = ForwardDiff.Chunk{min(length(θ), 8)}()
#chunk = ForwardDiff.Chunk{1}()
aif(x) = ai_func(lmm, data, x, beta)
grf(x) = reml_sweep_β(lmm, data, x, beta)[1]
aihcfg = ForwardDiff.HessianConfig(aif, θ, chunk)
Expand Down
132 changes: 0 additions & 132 deletions src/gmat.jl
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
###############################################################################
# G MATRIX FUNCTIONS
################################################################################

################################################################################
@noinline function gmatvec(θ::Vector{T}, covstr) where T
gt = [Symmetric(zeros(T, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn]
Expand All @@ -26,29 +25,7 @@ end
end
mx
end

################################################################################
#=
@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!(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)
end
end
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 Expand Up @@ -248,112 +225,3 @@ function tpnum(m, n, s)
end
b -= s - n
end
#=
function tpnum(m, n, s)
div(m*(2s - 1 - m), 2) - s + n
end
=#

################################################################################
# Grads

function gmat_g!(mx, θ, g::Int, ct::AbstractCovarianceType)
T = ForwardDiff.Dual{Nothing, Float64, 1}
gθ = Vector{T}(undef, length(θ))
for i = 1:length(θ)
if i == g gθ[i] = ForwardDiff.Dual(θ[i], 1) else gθ[i] = ForwardDiff.Dual(θ[i], 0) end
end
gmx = zeros(T, size(mx))
Metida.gmat!(gmx, gθ, ct)
for n = 1:size(gmx, 2)
for m = 1:n
mx[m, n] = ForwardDiff.partials(gmx[m, n])[1]
end
end
mx
end
#ZERO
function gmat_g!(mx, ::Any, ::Int, ::ZERO)
mx
end
#SI
function gmat_g!(mx, θ, ::Int, ::SI_)
val = 2θ[1]
@inbounds @simd for i = 1:size(mx, 1)
mx[i, i] = val
end
mx
end
#DIAG
function gmat_g!(mx, θ, g::Int, ::DIAG_)
for i = 1:size(mx, 1)
if i == g mx[i, i] = 2θ[i] end
end
mx
end

"""
Number of parameters correspondig to random effects (G) in theta vector.
"""
function thetagnum(covstr)
covstr.tl - covstr.t[end]
end
function grad_gmatvec(θ::Vector{T}, covstr) where T
gi = thetagnum(covstr) # total number of grads for G side
gt = Vector{Symmetric{T, Matrix{T}}}(undef, gi)
gn = 1
for r = 1:covstr.rn
if covstr.random[r].covtype.z
for g = 1:covstr.t[r]
gt[gn] = Symmetric(zeros(T, covstr.q[r], covstr.q[r]))
gmat_g!(gt[gn].data, view(θ, covstr.tr[r]), g, covstr.random[r].covtype.s)
gn += 1
end
end
end
gt
end

function grad_zgz_base_inc!(mx::AbstractArray, G, gn, covstr, bi)
if gn > thetagnum(covstr) return mx end
block = covstr.vcovblock[bi]
if covstr.random[1].covtype.z
r = covstr.emap[gn]
zblock = view(covstr.z, block, covstr.zrndur[r])
for i = 1:subjn(covstr, r, bi)
sb = getsubj(covstr, r, bi, i)
mulαβαtinc!(view(mx, sb, sb), view(zblock, sb, :), G[gn])
end
end
mx
end

@noinline function grad_rmat_base_inc!(mx, θ, gn, covstr, bi)
tgn = thetagnum(covstr)
if gn <= tgn return mx end
rn = gn - tgn
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_g!(view(mx, sb, sb), θ, view(zblock, sb, :), rn, covstr.repeated.covtype.s)
end
mx
end

function grad_vmatrix(G, rθ::AbstractVector{T}, lmm, i::Int) where T
gv = Vector{Symmetric{T, Matrix{T}}}(undef, lmm.covstr.tl)
for gi = 1:lmm.covstr.tl

q = length(lmm.covstr.vcovblock[i])

gv[gi] = Symmetric(zeros(T, q, q))

grad_zgz_base_inc!(gv[gi].data, G, gi, lmm.covstr, i)

grad_rmat_base_inc!(gv[gi].data, rθ, gi, lmm.covstr, i)
end
gv
end

Loading

2 comments on commit 959dc81

@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/74733

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.3 -m "<description of version>" 959dc811e52db3f1c0300a1fe2526cf3e486d350
git push origin v0.14.3

Please sign in to comment.