UN

PharmCat · Jul 29, 2022 · cf0cd1d · cf0cd1d
1 parent 0a67bc7
commit cf0cd1d
Show file tree

Hide file tree

Showing 12 changed files with 122 additions and 76 deletions.
diff --git a/Project.toml b/Project.toml
@@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
 keywords = ["lenearmodel", "mixedmodel"]
 desc = "Mixed-effects models with flexible covariance structure."
 authors = ["Vladimir Arnautov <[email protected]>"]
-version = "0.13.0"
+version = "0.14.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

diff --git a/change.log b/change.log
@@ -1,5 +1,9 @@
-v0.13
+v0.14
+  * Unstructured Covariance (UN)
+  * initial functions for analytical gradient vector
+  * docs
 
+v0.13
   * Bootstrap
   * MI
   * Changes in types

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -27,6 +27,7 @@ Implemented covariance structures:
   * Spatial Exponential (SPEXP)
   * Spatial Power (SPPOW)
   * Spatial Gaussian (SPGAU)
+  * Unstructured (UN) 
   * Custom Covariance Type
 
 ## Limitations

diff --git a/src/Metida.jl b/src/Metida.jl
@@ -16,7 +16,7 @@ import MetidaBase.PrettyTables: TextFormat, pretty_table, tf_borderless, ft_prin
 import LinearAlgebra:checksquare, BlasFloat
 import StatsModels: @formula, termvars, ModelFrame, FunctionTerm, AbstractTerm, CategoricalTerm, AbstractContrasts, ConstantTerm, InterceptTerm, Term, InteractionTerm, FormulaTerm, ModelMatrix, schema, apply_schema, MatrixTerm, modelcols
 import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof, loglikelihood, aic, bic, aicc, isfitted, vcov, mean, var, stderror, modelmatrix, response, CoefTable, coeftable
-import Base:show, rand, ht_keyindex
+import Base:show, rand, ht_keyindex, getproperty
 import Random: default_rng, AbstractRNG, rand!
 
 export @formula, @covstr,

diff --git a/src/gmat.jl b/src/gmat.jl
@@ -8,7 +8,7 @@ function gmat_switch!(G, θ, covstr, r)
     G
 end
 function gmatvec(θ::Vector{T}, covstr) where T
-    gt = Tuple(Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn)
+    gt = [Symmetric(Matrix{T}(undef, covstr.q[r], covstr.q[r])) for r in 1:covstr.rn] #<: fix to vec
     for r = 1:covstr.rn
         fill!(gt[r].data, zero(T))
         if covstr.random[r].covtype.z
@@ -18,9 +18,8 @@ function gmatvec(θ::Vector{T}, covstr) where T
     gt
 end
 ################################################################################
-function zgz_base_inc!(mx::AbstractArray{T}, θ::AbstractArray{T}, covstr, block, sblock) where T
+function zgz_base_inc!(mx::AbstractArray, θ::AbstractArray{T}, covstr, block, sblock) where T
     if covstr.random[1].covtype.z
-        #length of random to covstr
         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)
@@ -30,24 +29,17 @@ function zgz_base_inc!(mx::AbstractArray{T}, θ::AbstractArray{T}, covstr, block
             end
         end
     end
-    #    fill!(mx, zero(T))
-    #end
     mx
 end
-function zgz_base_inc!(mx::AbstractArray{T}, G, θ::AbstractArray{T}, covstr, block, sblock) where T
+function zgz_base_inc!(mx::AbstractArray, G, θ, covstr, block, sblock)
     if covstr.random[1].covtype.z
-        #length of random to covstr
         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)
             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])
             end
         end
     end
-    #    fill!(mx, zero(T))
-    #end
     mx
 end
 ################################################################################
@@ -64,18 +56,15 @@ Base.@propagate_inbounds function gmat!(mx, θ, ::SI_)
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = val
     end
-    nothing
+    mx
 end
 #DIAG
 function gmat!(mx, θ, ::DIAG_)
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = θ[i] ^ 2
     end
-    nothing
+    mx
 end
-#function gmat_vc!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
-#    nothing
-#end
 #AR
 function gmat!(mx, θ, ::AR_)
     de  = θ[1] ^ 2
@@ -90,7 +79,7 @@ function gmat!(mx, θ, ::AR_)
             end
         end
     end
-    nothing
+    mx
 end
 #ARH
 function gmat!(mx, θ, ::ARH_)
@@ -108,7 +97,7 @@ function gmat!(mx, θ, ::ARH_)
     @inbounds @simd for m = 1:s
         mx[m, m] = mx[m, m] * mx[m, m]
     end
-    nothing
+    mx
 end
 #CS
 function gmat!(mx, θ, ::CS_)
@@ -121,7 +110,7 @@ function gmat!(mx, θ, ::CS_)
             end
         end
     end
-    nothing
+    mx
 end
 #CSH
 function gmat!(mx, θ, ::CSH_)
@@ -139,7 +128,7 @@ function gmat!(mx, θ, ::CSH_)
     @inbounds @simd for m = 1:s
         mx[m, m] = mx[m, m] * mx[m, m]
     end
-    nothing
+    mx
 end
 ################################################################################
 #ARMA
@@ -156,7 +145,7 @@ function gmat!(mx, θ, ::ARMA_)
             end
         end
     end
-    nothing
+    mx
 end
 #TOEP
 function gmat!(mx, θ, ::TOEP_)
@@ -172,7 +161,7 @@ function gmat!(mx, θ, ::TOEP_)
             end
         end
     end
-    nothing
+    mx
 end
 function gmat!(mx, θ, ct::TOEPP_)
     de  = θ[1] ^ 2    #diagonal element
@@ -187,7 +176,7 @@ function gmat!(mx, θ, ct::TOEPP_)
             end
         end
     end
-    nothing
+    mx
 end
 #TOEPH
 function gmat!(mx, θ, ::TOEPH_)
@@ -205,7 +194,7 @@ function gmat!(mx, θ, ::TOEPH_)
     @inbounds @simd for m = 1:s
         mx[m, m] = mx[m, m] * mx[m, m]
     end
-    nothing
+    mx
 end
 #TOEPHP
 function gmat!(mx, θ, ct::TOEPHP_)
@@ -223,7 +212,7 @@ function gmat!(mx, θ, ct::TOEPHP_)
     @inbounds @simd for m = 1:s
         mx[m, m] = mx[m, m] * mx[m, m]
     end
-    nothing
+    mx
 end
 #UN
 function gmat!(mx, θ, ::UN_)
@@ -242,7 +231,7 @@ function gmat!(mx, θ, ::UN_)
         v = mx[m, m]
         mx[m, m] = v * v
     end
-    nothing
+    mx
 end
 
 function tpnum(m, n, s)
@@ -252,3 +241,44 @@ function tpnum(m, n, s)
     end
     b -= s - n
 end
+
+################################################################################
+
+# Grads
+#=
+function sim_gmat(s, θ::AbstractVector{T}, ct::AbstractCovarianceType) where T
+    mxt = zeros(T, s, s)
+    Symmetric(gmat!(mxt, θ, ct))
+end
+
+function grad_gmat!(s, θ, ct::AbstractCovarianceType)
+    mxg = zeros(s, s, length(θ))
+    gm  = ForwardDiff.jacobian!(mxg, x-> sim_gmat(s, x, ct), θ)
+    mxg
+end
+
+function grad_gmatvec(θ::Vector{T}, covstr) where T
+    gt = [Matrix{T}(undef, covstr.q[r], covstr.q[r], length(θ)) for r in 1:covstr.rn] #<: fix to vec
+    for r = 1:covstr.rn
+        if covstr.random[r].covtype.z
+            gt[r] = grad_gmat!(covstr.q[r], view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
+        end
+    end
+    gt
+end
+
+function grad_zgz_base_inc!(mx::AbstractArray, G, θ::AbstractArray{T}, covstr, block, sblock) where T
+    if covstr.random[1].covtype.z
+        gl = size(mx, 3)
+        for r = 1:covstr.rn
+            zblock    = view(covstr.z, block, covstr.zrndur[r])
+            @inbounds for i = 1:length(sblock[r])
+                for gn = 1:gl
+                    mulαβαtinc!(view(mx, sblock[r][i], sblock[r][i], gn), view(zblock, sblock[r][i], :), view(G[r], :, :, gn))
+                end
+            end
+        end
+    end
+    mx
+end
+=#
diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl
@@ -19,7 +19,7 @@ end
 
 Change θ (only upper triangle). B is symmetric.
 """
-function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix) where T
+function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     axb  = axes(B, 1)
     sa   = size(A, 1)
     for m ∈ 1:sa
@@ -38,7 +38,7 @@ end
 
 Change θ (only upper triangle). B is symmetric.
 """
-function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix, B::AbstractMatrix, alpha) where T
+function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, alpha)
     if  !(size(B, 1) == size(B, 2) == size(A, 2)) || !(size(A, 1) == size(θ, 1) == size(θ, 2)) throw(ArgumentError("Wrong dimentions!")) end
     axb  = axes(B, 1)
     sa   = size(A, 1)
@@ -58,7 +58,7 @@ end
 
 Change θ (only upper triangle). B is symmetric.
 """
-function mulαβαtinc!(θ::AbstractVector{T}, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) where T
+function mulαβαtinc!(θ::AbstractVector, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha)
     if !(size(B, 2) == length(a) == length(b)) || size(B, 1) != size(A, 2) || size(A, 1) != length(θ) throw(ArgumentError("Wrong dimentions.")) end
     axb  = axes(B, 1)
     sa   = size(A, 1)

diff --git a/src/lmm.jl b/src/lmm.jl
@@ -287,3 +287,14 @@ Return fitting log.
 function getlog(lmm::LMM)
     lmm.log
 end
+
+################################################################################
+
+function Base.getproperty(x::LMM, s::Symbol)
+    if s == :θ
+        return x.result.theta
+    elseif s == :β
+        return x.result.beta
+    end
+    getfield(x, s)
+end
diff --git a/src/reml.jl b/src/reml.jl
@@ -33,12 +33,7 @@ end
 ################################################################################
 #                     REML without provided β
 ################################################################################
-#=
-function reml_sweep_β(lmm, θ::Vector{T}; syrkblas::Bool = false) where T <: Number
-    data = LMMDataViews(lmm)
-    reml_sweep_β(lmm, data, θ; syrkblas = syrkblas)
-end
-=#
+
 function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T
     n             = length(lmm.covstr.vcovblock)
     N             = length(lmm.data.yv)
@@ -110,29 +105,23 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; syrkblas::Bool = false) where T
         if issuccess(cθs₂)
             # β calculation
             mul!(β, inv(cθs₂), βm)
-        # final θ₂
+            # final θ₂
             logdetθ₂ = logdet(cθs₂)
+            # θ₃
             @inbounds @simd for i = 1:n
                 θ₃ += mulθ₃(data.yv[i], data.xv[i], β, V⁻¹[i])
             end
         else
             β       .= NaN
             return   Inf, β, Inf, Inf, false
         end
-        # θ₃
-
 
     return   θ₁ + logdetθ₂ + θ₃ + c, β, θs₂, θ₃, noerror #REML, β, iC, θ₃, errors
 end
 ################################################################################
 #                     REML with provided β
 ################################################################################
-#=
-function reml_sweep_β(lmm, θ::Vector{T}, β::Vector) where T <: Number
-    data = LMMDataViews(lmm)
-    reml_sweep_β(lmm, data, θ, β)
-end
-=#
+
 function core_sweep_β(lmm, data, θ::Vector{T}, β, n) where T
     ncore     = min(num_cores(), n, 16)
     accθ₁     = zeros(T, ncore)