Skip to content

Commit

Permalink
clean code
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Dec 28, 2022
1 parent 6096e01 commit a0cc527
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

0 comments on commit a0cc527

Please sign in to comment.