Skip to content

Commit

Permalink
v0.2.2
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Feb 2, 2021
1 parent 8d9a24f commit f624fff
Show file tree
Hide file tree
Showing 12 changed files with 235 additions and 86 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Metida"
uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
authors = ["Vladimir Arnautov <[email protected]>"]
version = "0.2.1"
version = "0.2.2"

[deps]
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
Expand Down
8 changes: 8 additions & 0 deletions change.log
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,11 @@
- 0.2.2
* random effect zero fix
* ai-like algo fix
* dof-satter
* ranx
* rholinkf
* docs

- 0.2.1
* StatsBase methods
* ARMA
Expand Down
24 changes: 24 additions & 0 deletions docs/src/api.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,25 @@ Metida.CompoundSymmetry
```@docs
Metida.HeterogeneousCompoundSymmetry
```
### StatsBase.dof
```@docs
StatsBase.dof
```

### StatsBase.dof_residual
```@docs
StatsBase.dof_residual
```

### Metida.dof_satter
```@docs
Metida.dof_satter
```

### Metida.caic
```@docs
Metida.caic
```

### Metida.gmatrix
```@docs
Expand All @@ -59,3 +78,8 @@ Metida.gmatrix
```@docs
Metida.rmatrix
```

### Metida.vmatrix!
```@docs
Metida.vmatrix!
```
3 changes: 2 additions & 1 deletion src/Metida.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ AR, Autoregressive,
ARH, HeterogeneousAutoregressive,
CS, CompoundSymmetry,
CSH, HeterogeneousCompoundSymmetry,
fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength, dof_satter
fit!, LMM, VarEffect, theta, logreml, m2logreml, thetalength, dof_satter, rankx, caic,
gmatrix, rmatrix, vmatrix!

include("abstracttype.jl")
include("sweep.jl")
Expand Down
72 changes: 35 additions & 37 deletions src/fit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,12 +96,17 @@ function fit!(lmm::LMM{T};
end
#dfc = TwiceDifferentiableConstraints(lx, ux)
#Initial step with modified Newton method
chunk = ForwardDiff.Chunk{1}()
############################################################################
if aifirst
aif(x) = optfunc(lmm, x)[4]

beta = reml_sweep_β_c(lmm, θ)
aif(x) = reml_sweep_ai(lmm, x, beta)
grf(x) = optfunc(lmm, x)[1]
ai = ForwardDiff.hessian(aif, θ)
gr = ForwardDiff.gradient(grf, θ)
aihcfg = ForwardDiff.HessianConfig(aif, θ, chunk)
aigcfg = ForwardDiff.GradientConfig(grf, θ, chunk)
ai = ForwardDiff.hessian(aif, θ, aihcfg, Val{false}())
gr = ForwardDiff.gradient(grf, θ, aigcfg, Val{false}())
initθ = copy(θ)
try
θ .-= (inv(ai) ./4 )*gr
Expand All @@ -127,11 +132,12 @@ function fit!(lmm::LMM{T};
#Twice differentiable object

vloptf(x) = optfunc(lmm, varlinkvecapply(x, lmm.covstr.ct; rholinkf = rholinkf))[1]
chunk = ForwardDiff.Chunk{1}()
gcfg = ForwardDiff.GradientConfig(vloptf, θ, chunk)
hcfg = ForwardDiff.HessianConfig(vloptf, θ, chunk)
gfunc!(g, x) = ForwardDiff.gradient!(g, vloptf, x, gcfg, Val{false}())
hfunc!(h, x) = ForwardDiff.hessian!(h, vloptf, x, hcfg, Val{false}())
hfunc!(h, x) = begin
ForwardDiff.hessian!(h, vloptf, x, hcfg, Val{false}())
end
td = TwiceDifferentiable(vloptf, gfunc!, hfunc!, θ)
#td = TwiceDifferentiable(x ->optfunc(lmm, varlinkvecapply2(x, lmm.covstr.ct))[1], θ; autodiff = :forward)

Expand All @@ -147,43 +153,35 @@ function fit!(lmm::LMM{T};
#Theta (θ) vector
lmm.result.theta = varlinkvecapply!(deepcopy(Optim.minimizer(lmm.result.optim)), lmm.covstr.ct; rholinkf = rholinkf)
lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Resulting θ: "*string(lmm.result.theta)))
#try
if hcalck
#Hessian
lmm.result.h = hessian(lmm, lmm.result.theta)
#H positive definite check
if !isposdef(lmm.result.h)
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Hessian is not positive definite."))
end
qrd = qr(lmm.result.h, Val(true))
for i = 1:length(lmm.result.theta)
if abs(qrd.R[i,i]) < 1E-10
if lmm.covstr.ct[qrd.jpvt[i]] == :var
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Variation QR.R diagonal value ($(qrd.jpvt[i])) is less than 1e-10."))
elseif lmm.covstr.ct[qrd.jpvt[i]] == :rho
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Rho QR.R diagonal value ($(qrd.jpvt[i])) is less than 1e-10."))
end
end
end
end
try

#-2 LogREML, β, iC
lmm.result.reml, lmm.result.beta, iC, θ₃ = optfunc(lmm, lmm.result.theta)
lmm.result.reml, lmm.result.beta, iC, θ₃ = optfunc(lmm, lmm.result.theta)
#Variance-vovariance matrix of β
lmm.result.c = pinv(iC)
lmm.result.c = pinv(iC)
#SE
lmm.result.se = sqrt.(diag(lmm.result.c))
lmm.result.se = sqrt.(diag(lmm.result.c))
#Fit true
lmm.result.fit = true
lmm.result.fit = true

lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Model fitted."))
catch
#-2 LogREML, β, iC
lmm.result.reml, lmm.result.beta, iC, θ₃ = optfunc(lmm, lmm.result.theta)
#Fit false
lmm.result.fit = false

lmmlog!(io, lmm, verbose, LMMLogMsg(:ERROR, "Model not fitted!"))
if hcalck
#Hessian
lmm.result.h = hessian(lmm, lmm.result.theta)
#H positive definite check
if !isposdef(lmm.result.h)
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Hessian is not positive definite."))
end
qrd = qr(lmm.result.h, Val(true))
for i = 1:length(lmm.result.theta)
if abs(qrd.R[i,i]) < 1E-8
if lmm.covstr.ct[qrd.jpvt[i]] == :var
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Variation QR.R diagonal value ($(qrd.jpvt[i])) is less than 1e-10."))
elseif lmm.covstr.ct[qrd.jpvt[i]] == :rho
lmmlog!(io, lmm, verbose, LMMLogMsg(:WARN, "Rho QR.R diagonal value ($(qrd.jpvt[i])) is less than 1e-10."))
end
end
end
end

lmmlog!(io, lmm, verbose, LMMLogMsg(:INFO, "Model fitted."))
lmm
end
44 changes: 16 additions & 28 deletions src/gmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,12 @@
function gmat_base(θ::Vector{T}, covstr) where T
q = size(covstr.z, 2)
mx = zeros(T, q, q)
for i = 1:length(covstr.random)
vmxr = (1 + sum(covstr.q[1:i]) - covstr.q[i]):sum(covstr.q[1:i])
vmx = view(mx, vmxr, vmxr)
gmat_switch!(vmx, θ, covstr, i)
if covstr.random[1].covtype.s != :ZERO
for i = 1:length(covstr.random)
vmxr = (1 + sum(covstr.q[1:i]) - covstr.q[i]):sum(covstr.q[1:i])
vmx = view(mx, vmxr, vmxr)
gmat_switch!(vmx, θ, covstr, i)
end
end
mx
end
Expand All @@ -27,47 +29,33 @@ function gmat_switch!(G, θ, covstr, r)
gmat_cs!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype)
elseif covstr.random[r].covtype.s == :ARMA
gmat_arma!(G, θ[covstr.tr[r]], covstr.q[r], covstr.random[r].covtype)
elseif covstr.random[r].covtype.s == :ZERO
gmat_zero!(G, similar(θ, 0), covstr.q[r], covstr.random[r].covtype)
else
error("Unknown covariance structure!")
end
G
end
################################################################################
#=
function gmat_base_z!(mx, θ::Vector{T}, covstr) where T
q = sum(length.(covstr.block[1]))
for r = 1:length(covstr.random)
G = zeros(T, covstr.q[r], covstr.q[r])
gmat_switch!(G, θ, covstr, r)
for i = 1:length(covstr.block[r])
mulαβαtinc!(view(mx, covstr.block[r][i], covstr.block[r][i]), view(covstr.z, covstr.block[r][i], covstr.zrndur[r]), G)
end
end
mx
end
=#
################################################################################
function zgz_base_inc!(mx, θ::Vector{T}, covstr, block, sblock) where T
q = sum(length.(covstr.block[1]))
for r = 1:length(covstr.random)
G = zeros(T, covstr.q[r], covstr.q[r])
gmat_switch!(G, θ, covstr, r)
zblock = view(covstr.z, block, covstr.zrndur[r])
for i = 1:length(sblock[r])
mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i]), view(zblock, sblock[r][i], :), G)
if covstr.random[1].covtype.s != :ZERO
for r = 1:length(covstr.random)
G = zeros(T, covstr.q[r], covstr.q[r])
gmat_switch!(G, θ, covstr, r)
zblock = view(covstr.z, block, covstr.zrndur[r])
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_zero!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
mx .= zero(T)
nothing
end
=#
function gmat_si!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
val = θ[1] ^ 2
for i = 1:size(mx, 1)
Expand Down
2 changes: 1 addition & 1 deletion src/linearalgebra.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ end
"""
(y - X * β)' * V * (y - X * β)
"""
function mulθ₃(y::AbstractVector, X::AbstractMatrix, β::AbstractVector, V::AbstractMatrix{T})::T where T
function mulθ₃(y, X, β, V::AbstractMatrix{T})::T where T
q = size(V, 1)
p = size(X, 2)
θ = zero(T)
Expand Down
105 changes: 100 additions & 5 deletions src/reml.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ function reml_sweep_β_b(lmm, θ::Vector{T}) where T <: Number
try
θ₁ += logdet(cholesky(V))
catch
θ₁ = Inf
return (Inf, nothing, nothing, Inf)
end
sweep!(Vp, 1:q)
V⁻¹[i] = Symmetric(utriaply!(x -> -x, V))
Expand All @@ -46,11 +46,10 @@ function reml_sweep_β_b(lmm, θ::Vector{T}) where T <: Number
try
logdetθ₂ = logdet(θ₂)
catch
logdetθ₂ = Inf
return (Inf, nothing, nothing, Inf)
end
return θ₁ + logdetθ₂ + θ₃ + c, β, θ₂, θ₃
end

function reml_sweep_β(lmm, θ::Vector{T}) where T <: Number
n = length(lmm.data.block)
N = length(lmm.data.yv)
Expand All @@ -76,8 +75,9 @@ function reml_sweep_β(lmm, θ::Vector{T}) where T <: Number
V = view(Vp, 1:q, 1:q)
Vx = view(Vp, 1:q, q+1:q+lmm.rankx)
Vx .= view(lmm.data.xv, lmm.data.block[i],:)
zgz_base_inc!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])
rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])
vmatrix!(V, θ, lmm, i)
#zgz_base_inc!(V, θ, lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])
#rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.data.block[i], lmm.covstr.sblock[i])
#-----------------------------------------------------------------------
if isposdef(V)
θ₁ += logdet(cholesky(V))
Expand All @@ -101,3 +101,98 @@ function reml_sweep_β(lmm, θ::Vector{T}) where T <: Number
end
################################################################################
################################################################################
function reml_sweep_β(lmm, θ::Vector{T}, β) where T <: Number
n = length(lmm.data.block)
N = length(lmm.data.yv)
c = (N - lmm.rankx)*log(2π)
#---------------------------------------------------------------------------
# Vector log determinant of V matrix
θ₁ = zero(T)
θ₂ = zeros(T, lmm.rankx, lmm.rankx)
θ₃ = zero(T)
#---------------------------------------------------------------------------
q = zero(Int)
qswm = zero(Int)
Vp = Matrix{T}(undef, 0, 0)
logdetθ₂ = zero(T)
@inbounds for i = 1:n
q = length(lmm.data.block[i])
qswm = q + lmm.rankx
Vp = zeros(T, q + lmm.rankx, q + lmm.rankx)
V = view(Vp, 1:q, 1:q)
Vx = view(Vp, 1:q, q+1:q+lmm.rankx)
Vx .= view(lmm.data.xv, lmm.data.block[i],:)
vmatrix!(V, θ, lmm, i)
#-----------------------------------------------------------------------
if isposdef(V)
θ₁ += logdet(cholesky(V))
else
return (Inf, nothing, nothing, Inf)
end
sweep!(Vp, 1:q)
V⁻¹ = Symmetric(utriaply!(x -> -x, V))
#-----------------------------------------------------------------------
qswm = size(Vp, 1)
θ₂ .-= Symmetric(view(Vp, q + 1:qswm, q + 1:qswm))
#-----------------------------------------------------------------------
@inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹)
end
logdetθ₂ = logdet(θ₂)
return θ₁ + logdetθ₂ + θ₃ + c, θ₂, θ₃ #REML, iC, θ₃
end
function reml_sweep_ai(lmm, θ::Vector{T}, β) where T <: Number
n = length(lmm.data.block)
N = length(lmm.data.yv)
c = (N - lmm.rankx)*log(2π)
θ₃ = zero(T)
q = zero(Int)
qswm = zero(Int)
Vp = Matrix{T}(undef, 0, 0)
@inbounds for i = 1:n
q = length(lmm.data.block[i])
qswm = q + lmm.rankx
Vp = zeros(T, q + lmm.rankx, q + lmm.rankx)
V = view(Vp, 1:q, 1:q)
Vx = view(Vp, 1:q, q+1:q+lmm.rankx)
Vx .= view(lmm.data.xv, lmm.data.block[i],:)
vmatrix!(V, θ, lmm, i)
sweep!(Vp, 1:q)
V⁻¹ = Symmetric(utriaply!(x -> -x, V))
@inbounds θ₃ += mulθ₃(view(lmm.data.yv, lmm.data.block[i]), view(lmm.data.xv, lmm.data.block[i],:), β, V⁻¹)
end
return θ₃
end
function reml_sweep_β_c(lmm, θ::Vector{T}) where T <: Number
n = length(lmm.data.block)
N = length(lmm.data.yv)
c = (N - lmm.rankx)*log(2π)
#---------------------------------------------------------------------------
V⁻¹ = Vector{Matrix{T}}(undef, n)
# Vector log determinant of V matrix
θ₂ = zeros(T, lmm.rankx, lmm.rankx)
βm = zeros(T, lmm.rankx)
β = Vector{T}(undef, lmm.rankx)
#---------------------------------------------------------------------------
q = zero(Int)
qswm = zero(Int)
Vp = Matrix{T}(undef, 0, 0)
@inbounds for i = 1:n
q = length(lmm.data.block[i])
qswm = q + lmm.rankx
Vp = zeros(T, q + lmm.rankx, q + lmm.rankx)
V = view(Vp, 1:q, 1:q)
Vx = view(Vp, 1:q, q+1:q+lmm.rankx)
Vx .= view(lmm.data.xv, lmm.data.block[i],:)
vmatrix!(V, θ, lmm, i)
#-----------------------------------------------------------------------
sweep!(Vp, 1:q)
V⁻¹[i] = Symmetric(utriaply!(x -> -x, V))
#-----------------------------------------------------------------------
qswm = size(Vp, 1)
θ₂ .-= Symmetric(view(Vp, q + 1:qswm, q + 1:qswm))
mulαtβinc!(βm, view(Vp, 1:q, q + 1:qswm), view(lmm.data.yv, lmm.data.block[i]))
#-----------------------------------------------------------------------
end
mul!(β, inv(θ₂), βm)
return β
end
Loading

2 comments on commit f624fff

@PharmCat
Copy link
Owner

Choose a reason for hiding this comment

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

@JuliaRegistrator register

Release notes:

  • 0.2.2
    • random effect zero fix
    • ai-like algo fix
    • dof-satter
    • ranx
    • rholinkf
    • docs

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

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.2.2 -m "<description of version>" f624fff05fa7431d60f5ff6179425d022b733cfe
git push origin v0.2.2

Please sign in to comment.