diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
index 5ca7853..21695a8 100644
--- a/.github/workflows/Documenter.yml
+++ b/.github/workflows/Documenter.yml
@@ -23,12 +23,11 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 45
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: "1.8"
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-docdeploy@v1
+      - uses: julia-actions/cache@v1
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy
diff --git a/.github/workflows/Tier1.yml b/.github/workflows/Tier1.yml
index 7639b18..9431170 100644
--- a/.github/workflows/Tier1.yml
+++ b/.github/workflows/Tier1.yml
@@ -39,15 +39,26 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v3
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: julia-actions/cache@v1
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v3
+      - uses: codecov/codecov-action@v4
+        if: ${{ matrix.os == 'ubuntu-latest' && matrix.version == '1' && matrix.arch == 'x64' }}
         with:
-          files: lcov.info
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/Project.toml b/Project.toml
index 86d8102..3f7ad97 100644
--- 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 <mail@pharmcat.net>"]
-version = "0.15.1"
+version = "0.16.0"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
diff --git a/docs/src/api.md b/docs/src/api.md
index b2c824b..ee1300e 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -124,6 +124,12 @@ Metida.ToeplitzParameterized
 Metida.Unstructured
 ```
 
+### Metida.ScaledWeightedCov
+
+```@docs
+Metida.ScaledWeightedCov
+```
+
 ### Methods
 
 ### Metida.caic
diff --git a/docs/src/custom.md b/docs/src/custom.md
index 7789d09..6ad8de7 100644
--- a/docs/src/custom.md
+++ b/docs/src/custom.md
@@ -26,10 +26,10 @@ function Metida.gmat!(mx, θ, ::YourCovarianceStruct)
 end
 ```
 
-Function `rmat!` have 4 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
+Function `rmat!` have 5 arguments and add repeated effect to V': V = V' + R (so V = Z * G * Z' + R), `mx` - V' matrix, `θ` - theta vector for this effect, `rz` - subject effect matrix, `ct` - your covariance type object, `sb` = block number. For example, `rmat!` for Heterogeneous Toeplitz Parameterized structure is specified bellow (`TOEPHP_  <: AbstractCovarianceType`).
 
 ```
-function Metida.rmat!(mx, θ, rz, ct::TOEPHP_)
+function Metida.rmat!(mx, θ, rz, ct::TOEPHP_, ::Int)
     l     = size(rz, 2)
     vec   = rz * (θ[1:l])
     s   = size(mx, 1)
@@ -123,7 +123,7 @@ Metida.fit!(lmm)
 
 # for R matrix
 
-function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure)
+function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure, ::Int)
     vec = Metida.tmul_unsafe(rz, θ)
     rn    = size(mx, 1)
     if rn > 1
diff --git a/docs/src/details.md b/docs/src/details.md
index 4fceee5..a6959f5 100644
--- a/docs/src/details.md
+++ b/docs/src/details.md
@@ -13,6 +13,33 @@ In matrix notation a mixed effect model can be represented as:
 y = X\beta + Zu + \epsilon
 ```
 
+where:
+
+
+```math
+\epsilon \sim N(0, R)
+
+\\
+
+u \sim N(0, G)
+
+\\
+
+y \sim N(X\beta, V)
+
+```
+
+where V depends on covariance sructure and parameters ``\theta``:
+
+```math
+V = CovStruct(\theta)
+
+```
+
+The unknown parameters include the regression parameters in ``\beta`` and covariance parameters in ``\theta``.
+
+Estimation of these model parameters relies on the use of a Newton-Ralphson (by default) algorithm. When we use either algorithm for finding REML solutions, we need to compute ``V^{-1}`` and its derivatives with respect to ``\theta``, which are computationally difficult for large ``n``, therefor SWEEP (see https://github.com/joshday/SweepOperator.jl) algorithm used to meke oprtimization less computationaly expensive.
+
 #### V
 
 ```math
@@ -33,7 +60,7 @@ logREML(\theta,\beta) = -\frac{N-p}{2} - \frac{1}{2}\sum_{i=1}^nlog|V_{\theta, i
 -\frac{1}{2}log|\sum_{i=1}^nX_i'V_{\theta, i}^{-1}X_i|-\frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_{\theta, i}^{-1}(y_i - X_{i}\beta)
 ```
 
-Actually ```L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c``` used for optimization, where:
+Actually `` L(\theta) = -2logREML = L_1(\theta) + L_2(\theta) + L_3(\theta) + c`` used for optimization, where:
 
 ```math
 L_1(\theta) = \frac{1}{2}\sum_{i=1}^nlog|V_{i}| \\
@@ -51,16 +78,36 @@ L_3(\theta) = \frac{1}{2}\sum_{i=1}^n(y_i - X_{i}\beta)'V_i^{-1}(y_i - X_{i}\bet
 \mathcal{H}\mathcal{L}(\theta) =  \mathcal{H}L_1(\theta)  + \mathcal{H}L_2(\theta) +  \mathcal{H} L_3(\theta)
 ```
 
-#### Weights
+#### [Weights](@id weights_header)
 
 If weights defined:
 
 ```math
+
+W_{i} = diag(wts_{i})
+
+\\
+
 V_{i} = Z_{i} G Z_i'+ W^{- \frac{1}{2}}_i R_{i} W^{- \frac{1}{2}}_i
 ```
 
+where ``W`` - diagonal matrix of weights.
+
+
+If wts is matrix then:
+
+
+```math
+
+W_{i} = wts_{i}
+
+\\
+
+V_{i} = Z_{i} G Z_i'+  R_{i} \circ W_{i}
+```
+
+where ``\circ`` - element-wise multiplication.
 
-where ```W``` - diagonal matrix of weights.
 
 
 #### Multiple random and repeated effects
diff --git a/docs/src/index.md b/docs/src/index.md
index a364593..a42bf86 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -38,7 +38,7 @@ Implemented covariance structures:
 
 Actually Metida can fit datasets with wore than  160k observation and 40k subjects levels on PC with 64 GB RAM. This is not "hard-coded" limitation, but depends on your model and data structure. Fitting of big datasets can take a lot of time. Optimal dataset size is less than 100k observations with maximum length of block less than 400.
 
-!!! note
+!!! warning
 
     Julia v1.8 or higher required.
 
diff --git a/src/Metida.jl b/src/Metida.jl
index a86b214..e8950eb 100644
--- a/src/Metida.jl
+++ b/src/Metida.jl
@@ -17,8 +17,10 @@ import StatsBase: fit, fit!, coef, coefnames, confint, nobs, dof_residual, dof,
 import Base: show, rand, ht_keyindex, getproperty
 import Random: default_rng, AbstractRNG, rand!
 
+
 export @formula, @covstr, @lmmformula,
 SI, ScaledIdentity,
+SWC, ScaledWeightedCov,
 DIAG, Diag,
 AR, Autoregressive,
 ARH, HeterogeneousAutoregressive,
diff --git a/src/fvalue.jl b/src/fvalue.jl
index 1a69e69..1d8e6ba 100644
--- a/src/fvalue.jl
+++ b/src/fvalue.jl
@@ -1,8 +1,5 @@
 # fvalue.jl
 
-#=
-Metida.fvalue(lmm, [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0])
-=#
 """
     fvalue(lmm::LMM, l::Matrix)
 
diff --git a/src/gmat.jl b/src/gmat.jl
index 452cf9f..d4b0418 100644
--- a/src/gmat.jl
+++ b/src/gmat.jl
@@ -9,7 +9,7 @@
             gmat!(gt[r].data, view(θ, covstr.tr[r]), covstr.random[r].covtype.s)
         end
     end
-    gt
+    return gt
 end
 # Main
 @noinline function zgz_base_inc!(mx::AbstractArray, G, covstr, bi)
@@ -23,7 +23,7 @@ end
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 ################################################################################
@@ -32,7 +32,7 @@ function gmat!(::Any, ::Any, ::AbstractCovarianceType)
     error("No gmat! method defined for thit structure!")
 end
 function gmat!(mx, ::Any, ::ZERO)
-    mx
+    return mx
 end
 #SI
 Base.@propagate_inbounds function gmat!(mx, θ, ::SI_)
@@ -40,14 +40,14 @@ Base.@propagate_inbounds function gmat!(mx, θ, ::SI_)
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = val
     end
-    mx
+    return mx
 end
 #DIAG
 function gmat!(mx, θ, ::DIAG_)
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] = θ[i] ^ 2
     end
-    mx
+    return mx
 end
 #AR
 function gmat!(mx, θ, ::AR_)
@@ -64,7 +64,7 @@ function gmat!(mx, θ, ::AR_)
             end
         end
     end
-    mx
+    return mx
 end
 #ARH
 function gmat!(mx, θ, ::ARH_)
@@ -84,7 +84,7 @@ function gmat!(mx, θ, ::ARH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m] 
     end
-    mx
+    return mx
 end
 #CS
 function gmat!(mx, θ, ::CS_)
@@ -99,7 +99,7 @@ function gmat!(mx, θ, ::CS_)
             end
         end
     end
-    mx
+    return mx
 end
 #CSH
 function gmat!(mx, θ, ::CSH_)
@@ -118,7 +118,7 @@ function gmat!(mx, θ, ::CSH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 ################################################################################
 #ARMA
@@ -136,7 +136,7 @@ function gmat!(mx, θ, ::ARMA_)
             end
         end
     end
-    mx
+    return mx
 end
 #TOEP
 function gmat!(mx, θ, ::TOEP_)
@@ -152,7 +152,7 @@ function gmat!(mx, θ, ::TOEP_)
             end
         end
     end
-    mx
+    return mx
 end
 function gmat!(mx, θ, ct::TOEPP_)
     de  = θ[1] ^ 2    #diagonal element
@@ -167,7 +167,7 @@ function gmat!(mx, θ, ct::TOEPP_)
             end
         end
     end
-    mx
+    return mx
 end
 #TOEPH
 function gmat!(mx, θ, ::TOEPH_)
@@ -186,7 +186,7 @@ function gmat!(mx, θ, ::TOEPH_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 #TOEPHP
 function gmat!(mx, θ, ct::TOEPHP_)
@@ -205,7 +205,7 @@ function gmat!(mx, θ, ct::TOEPHP_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 #UN
 function gmat!(mx, θ, ::UN_)
@@ -224,7 +224,7 @@ function gmat!(mx, θ, ::UN_)
     @inbounds @simd for m = 1:s
         mx[m, m] *= mx[m, m]
     end
-    mx
+    return mx
 end
 
 function tpnum(m, n, s)
@@ -233,4 +233,5 @@ function tpnum(m, n, s)
         b += s - i
     end
     b -= s - n
+    return b
 end
diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl
index 79d5093..1f5de27 100644
--- a/src/linearalgebra.jl
+++ b/src/linearalgebra.jl
@@ -23,7 +23,7 @@ Change θ (only upper triangle). B is symmetric.
             end
         end
     end
-    θ
+    return θ
 end
 #=
 function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}) where T <: AbstractFloat
@@ -55,7 +55,7 @@ Change θ (only upper triangle). B is symmetric.
             end
         end
     end
-    θ
+    return θ
 end
 """
     mulαβαtinc!(θ::AbstractVector{T}, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) where T
@@ -78,7 +78,7 @@ Change θ (only upper triangle). B is symmetric.
             end
         end
     end
-    θ
+    return θ
 end
 
 """
@@ -108,8 +108,8 @@ use only upper triangle of V
         end
     end
     for m = 2:q
-        @inbounds ycm2 = (y[m] - c[m])*2
-        @simd for n = 1:m-1
+        @inbounds ycm2 = (y[m] - c[m]) * 2
+        @simd for n = 1:m - 1
             @inbounds θ -= V[n, m] * (y[n] - c[n]) * ycm2
         end
     end
@@ -137,7 +137,7 @@ Change θ.
         end
         @inbounds θ[n] += θn
     end
-    θ
+    return θ
 end
 # Diagonal(b) * A * Diagonal(b) - chnage only A upper triangle 
 @noinline function mulβdαβd!(A::AbstractMatrix, b::AbstractVector)
@@ -149,7 +149,7 @@ end
             @inbounds A[m, n] *= b[m] * b[n]
         end
     end
-    A
+    return A
 end
 
 
@@ -162,7 +162,7 @@ end
             @inbounds vec[r] += rz[r, i] * θi
         end
     end
-    vec
+    return vec
 end
 
 @inline function diag!(f, v, m)
@@ -171,5 +171,5 @@ end
     @simd for i = 1:l
         @inbounds v[i] = f(m[i, i])
     end
-    v
+    return v
 end
diff --git a/src/lmm.jl b/src/lmm.jl
index cdeb5c7..9edbcc2 100644
--- a/src/lmm.jl
+++ b/src/lmm.jl
@@ -11,7 +11,7 @@ struct ModelStructure
 end
 
 """
-    LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing, wts::Union{Nothing, AbstractVector, AbstractString, Symbol} = nothing)
+    LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing, wts::Union{Nothing, AbstractVector, AbstractMatrix, AbstractString, Symbol} = nothing)
 
 Make Linear-Mixed Model object.
 
@@ -27,6 +27,10 @@ Make Linear-Mixed Model object.
 
 `wts`: regression weights (residuals).
 
+Weigts can be set as `Symbol` or `String`, in this case weights taken from tabular data.
+If weights is vector then this vector applyed to R-side part of covariance matrix (see [Weights details](@ref weights_header)).
+If weights is matrix then R-side part of covariance matrix multiplied by corresponding part of weight-matrix.
+
 See also: [`@lmmformula`](@ref)
 """
 struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
@@ -57,7 +61,11 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
         log::Vector{LMMLogMsg}) where T where W <: Union{LMMWts, Nothing}
         new{T, W}(model, f, modstr, covstr, data, dv, nfixed, rankx, result, maxvcbl, wts, log)
     end
-    function LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, wts::Union{Nothing, AbstractVector, AbstractString, Symbol} = nothing)
+    function LMM(model, data; 
+        contrasts=Dict{Symbol,Any}(),  
+        random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, 
+        repeated::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, 
+        wts::Union{Nothing, AbstractVector, AbstractMatrix, AbstractString, Symbol} = nothing)
         #need check responce - Float
         if !Tables.istable(data) error("Data not a table!") end
         if repeated === nothing && random === nothing
@@ -122,12 +130,22 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
             if wts isa Symbol
                 wts = Tables.getcolumn(data, wts)
             end
-            if length(lmmdata.yv) == length(wts)
-                if any(x -> x <= zero(x), wts) error("Only cases with positive weights allowed!") end
-                lmmwts = LMMWts(wts, covstr.vcovblock)
-            else
-                @warn "wts count not equal observations count! wts not used."
-                lmmwts = nothing
+            if wts isa AbstractVector
+                if length(lmmdata.yv) == length(wts)
+                    if any(x -> x <= zero(x), wts) error("Only cases with positive weights allowed!") end
+                    lmmwts = LMMWts(wts, covstr.vcovblock)
+                else
+                    @warn "wts count not equal observations count! wts not used."
+                    lmmwts = nothing
+                end
+            elseif wts isa AbstractMatrix
+                if length(lmmdata.yv) == LinearAlgebra.checksquare(wts)
+                    if any(x -> x <= zero(x), wts) error("Only positive values allowed!") end
+                    lmmwts = LMMWts(wts, covstr.vcovblock)
+                else
+                    @warn "wts count not equal observations count! wts not used."
+                    lmmwts = nothing
+                end
             end
         end
 
diff --git a/src/lmmdata.jl b/src/lmmdata.jl
index bf07fa7..98f6805 100644
--- a/src/lmmdata.jl
+++ b/src/lmmdata.jl
@@ -35,16 +35,23 @@ struct LMMDataViews{T<:AbstractFloat} <: AbstractLMMDataBlocks
     end
 end
 
-struct LMMWts{T<:AbstractFloat} 
-    sqrtwts::Vector{Vector{T}}
-    function LMMWts(sqrtwts::Vector{Vector{T}}) where T
+struct LMMWts{T} 
+    sqrtwts::Vector{T}
+    function LMMWts(sqrtwts::Vector{T}) where T
         new{T}(sqrtwts)
     end
-    function LMMWts(wts::Vector{T}, vcovblock) where T
+    function LMMWts(wts::AbstractVector{T}, vcovblock) where T
         sqrtwts = Vector{Vector{T}}(undef, length(vcovblock))
         for i in eachindex(vcovblock)
             sqrtwts[i] = @. inv(sqrt($(view(wts, vcovblock[i]))))
         end
         LMMWts(sqrtwts)
     end
+    function LMMWts(wts::AbstractMatrix{T}, vcovblock) where T
+        sqrtwts = Vector{Matrix{T}}(undef, length(vcovblock))
+        for i in eachindex(vcovblock)
+            sqrtwts[i] = wts[vcovblock[i], vcovblock[i]]
+        end
+        LMMWts(sqrtwts)
+    end
 end
\ No newline at end of file
diff --git a/src/reml.jl b/src/reml.jl
index 10c29e5..84c5fe1 100644
--- a/src/reml.jl
+++ b/src/reml.jl
@@ -7,7 +7,7 @@ function subutri!(a, b)
             @inbounds a[m,n] -= b[m,n]
         end
     end
-    a
+    return a
 end
 
 function fillzeroutri!(a::AbstractArray{T})  where T
@@ -28,7 +28,7 @@ function checkmatrix!(mx::AbstractMatrix{T}) where T
             e = false
         end
     end
-    e
+    return e
 end
 ################################################################################
 #                     REML without provided β
@@ -87,8 +87,20 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 4) where T #
             #-----------------------------------------------------------------------
         end
         θ₁      = sum(accθ₁)
-        θ₂      = sum(accθ₂)
-        βm      = sum(accβm)
+        #θ₂      = sum(accθ₂)
+        if length(accθ₂) > 1
+            for i = 2:length(accθ₂)
+                accθ₂[1] += accθ₂[i]
+            end
+        end
+        θ₂ = accθ₂[1]
+
+        if length(accβm) > 1
+            for i = 2:length(accβm)
+                accβm[1] += accβm[i]
+            end
+        end
+        βm = accβm[1]
         noerror = all(erroracc)
         noerror = noerror * checkmatrix!(θ₂)
         θs₂     = Symmetric(θ₂)
@@ -189,8 +201,20 @@ function reml_sweep_β_nlopt(lmm, data, θ::Vector{T}; maxthreads::Int = 16) whe
             return   Inf, β, θ₂, Inf, false
         end
         θ₁   = sum(accθ₁)
-        θ₂tc = sum(accθ₂)
-        βtc  = sum(accβm)
+        #θ₂tc = sum(accθ₂)
+        if length(accθ₂) > 1
+            for i = 2:length(accθ₂)
+                accθ₂[1] += accθ₂[i]
+            end
+        end
+        θ₂tc = accθ₂[1]
+        #βtc  = sum(accβm)
+        if length(accβm) > 1
+            for i = 2:length(accβm)
+                accβm[1] += accβm[i]
+            end
+        end
+        βtc = accβm[1]
     # Beta calculation
         copyto!(θ₂, θ₂tc)
         ldθ₂, info = LinearAlgebra.LAPACK.potrf!('U', θ₂tc)
@@ -246,7 +270,7 @@ function core_sweep_β(lmm, data, θ::Vector{T}, β, n; maxthreads::Int = 16) wh
             accθ₃[t]  += mulθ₃(data.yv[i], data.xv[i], β, V)
         end
     end
-    sum(accθ₁), sum(accθ₂), sum(accθ₃), all(erroracc)
+    return sum(accθ₁), sum(accθ₂), sum(accθ₃), all(erroracc)
 end
 ###
 function reml_sweep_β(lmm, data, θ::Vector{T}, β; kwargs...) where T
diff --git a/src/rmat.jl b/src/rmat.jl
index 7681c7e..8c4a098 100644
--- a/src/rmat.jl
+++ b/src/rmat.jl
@@ -10,34 +10,49 @@
         zblock    = view(covstr.rz[j], block, :)
         @simd for i = 1:subjn(covstr, en, bi)
             sb = getsubj(covstr, en, bi, i)
-            rmat!(view(mx, sb, sb), view(θ, rθ[j]), view(zblock, sb, :), covstr.repeated[j].covtype.s)
+            rmat!(view(mx, sb, sb), view(θ, rθ[j]), view(zblock, sb, :), covstr.repeated[j].covtype.s, bi)
         end
     end
-    mx
+    return mx
 end
 ################################################################################
-function rmat!(::Any, ::Any, ::Any,  ::AbstractCovarianceType)
+function rmat!(::Any, ::Any, ::Any,  ::AbstractCovarianceType, ::Int)
     error("No rmat! method defined for thit structure!")
 end
 #SI
-Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_)
+Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_, ::Int)
     val = θ[1]^2
     @inbounds @simd for i ∈ axes(mx, 1)
             mx[i, i] += val
     end
-    mx
+    return mx
+end
+#SWC
+function rmat!(mx, θ, ::AbstractMatrix, ct::SWC_, sbj::Int)
+    s  = size(mx, 1)
+    de  = θ[1] ^ 2
+    if s > 1
+        for n = 1:s
+            @inbounds @simd for m = 1:n
+                mx[m, n] += de * ct.wtsb[sbj][m, n]
+            end
+        end
+    else
+        @inbounds mx[1, 1] += de * ct.wtsb[sbj][1, 1]
+    end
+    return mx
 end
 #DIAG
-function rmat!(mx, θ, rz, ::DIAG_)
+function rmat!(mx, θ, rz, ::DIAG_, ::Int)
     for i ∈ axes(mx, 1)
         @inbounds @simd for c ∈ axes(θ, 1)
             mx[i, i] += rz[i, c] * θ[c] ^ 2
         end
     end
-    mx
+    return mx
 end
 #AR
-function rmat!(mx, θ, ::AbstractMatrix, ::AR_)
+function rmat!(mx, θ, ::AbstractMatrix, ::AR_, ::Int)
     s  = size(mx, 1)
     de  = θ[1] ^ 2
     @inbounds @simd for m = 1:s
@@ -50,10 +65,10 @@ function rmat!(mx, θ, ::AbstractMatrix, ::AR_)
             end
         end
     end
-    mx
+    return mx
 end
 #ARH
-function rmat!(mx, θ, rz, ::ARH_)
+function rmat!(mx, θ, rz, ::ARH_, ::Int)
     vec = tmul_unsafe(rz, θ)
     s    = size(mx, 1)
     if s > 1
@@ -66,13 +81,13 @@ function rmat!(mx, θ, rz, ::ARH_)
     @inbounds  for m ∈ axes(mx, 1)
         mx[m, m] += vec[m] * vec[m]
     end
-    mx
+    return mx
 end
 #CS
-function rmat!(mx, θ, ::AbstractMatrix,  ::CS_)
+function rmat!(mx, θ, ::AbstractMatrix,  ::CS_, ::Int)
     s    = size(mx, 1)
-    θsq   =  θ[1]*θ[1]
-    θsqp  =  θsq*θ[2]
+    θsq   =  θ[1] * θ[1]
+    θsqp  =  θsq  * θ[2]
     @inbounds @simd for i = 1:size(mx, 1)
         mx[i, i] += θsq
     end
@@ -83,17 +98,17 @@ function rmat!(mx, θ, ::AbstractMatrix,  ::CS_)
             end
         end
     end
-    mx
+    return mx
 end
 #CSH
-function rmat!(mx, θ, rz, ::CSH_)
+function rmat!(mx, θ, rz, ::CSH_, ::Int)
     vec = tmul_unsafe(rz, θ)
     s    = size(mx, 1)
     if s > 1
         θend = last(θ)
         for n = 2:s
             @inbounds vecnθend = vec[n] * θend
-            @inbounds @simd for m = 1:n-1
+            @inbounds @simd for m = 1:n - 1
                 mx[m, n] += vec[m] * vecnθend
             end
         end
@@ -101,11 +116,11 @@ function rmat!(mx, θ, rz, ::CSH_)
     @inbounds  for m ∈ axes(mx, 1)
         mx[m, m] += vec[m] * vec[m]
     end
-    mx
+    return mx
 end
 ################################################################################
 #ARMA
-function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_)
+function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_, ::Int)
     s  = size(mx, 1)
     de  = θ[1] ^ 2
     @inbounds @simd for m = 1:s
@@ -119,11 +134,11 @@ function rmat!(mx, θ, ::AbstractMatrix,  ::ARMA_)
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 #TOEPP
-function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_)
+function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_, ::Int)
     de  = θ[1] ^ 2    #diagonal element
     s   = size(mx, 1) #size
     @inbounds @simd for i = 1:s
@@ -136,11 +151,11 @@ function rmat!(mx, θ, ::AbstractMatrix, ct::TOEPP_)
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 #TOEPHP
-function rmat!(mx, θ, rz, ct::TOEPHP_)
+function rmat!(mx, θ, rz, ct::TOEPHP_, ::Int)
     l     = size(rz, 2)
     vec   = rz * (θ[1:l])
     s   = size(mx, 1) #size
@@ -154,7 +169,7 @@ function rmat!(mx, θ, rz, ct::TOEPHP_)
     @inbounds @simd for m = 1:s
         mx[m, m] += vec[m] * vec[m]
     end
-    mx
+    return mx
 end
 ################################################################################
 #=
@@ -181,7 +196,7 @@ function edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T
 end
 ################################################################################
 #SPEXP
-function rmat!(mx, θ, rz,  ::SPEXP_)
+function rmat!(mx, θ, rz,  ::SPEXP_, ::Int)
     σ²    = θ[1]^2
     #θe    = exp(θ[2])
     θe    = θ[2]
@@ -198,11 +213,11 @@ function rmat!(mx, θ, rz,  ::SPEXP_)
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 #SPPOW
-function rmat!(mx, θ, rz,  ::SPPOW_)
+function rmat!(mx, θ, rz,  ::SPPOW_, ::Int)
     σ²    = θ[1]^2
     ρ     = θ[2]
     s    = size(mx, 1)
@@ -218,11 +233,11 @@ function rmat!(mx, θ, rz,  ::SPPOW_)
             end
         end
     end
-    mx
+    return mx
 end
 
 #SPGAU
-function rmat!(mx, θ, rz,  ::SPGAU_)
+function rmat!(mx, θ, rz,  ::SPGAU_, ::Int)
     σ²    = θ[1]^2
     #θe    = exp(θ[2])
     θe    = θ[2]
@@ -240,11 +255,11 @@ function rmat!(mx, θ, rz,  ::SPGAU_)
             end
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 #SPEXPD cos(pidij)
-function rmat!(mx, θ, rz,  ::SPEXPD_)
+function rmat!(mx, θ, rz,  ::SPEXPD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
@@ -260,10 +275,10 @@ function rmat!(mx, θ, rz,  ::SPEXPD_)
             end
         end
     end
-    mx
+    return mx
 end
 #SPPOWD
-function rmat!(mx, θ, rz,  ::SPPOWD_)
+function rmat!(mx, θ, rz,  ::SPPOWD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     ρ     = θ[3]
@@ -278,10 +293,10 @@ function rmat!(mx, θ, rz,  ::SPPOWD_)
             end
         end
     end
-    mx
+    return mx
 end
 #SPGAUD
-function rmat!(mx, θ, rz,  ::SPGAUD_)
+function rmat!(mx, θ, rz,  ::SPGAUD_, ::Int)
     σ²    = θ[2]^2
     σ²d   = θ[1]^2 + σ²
     θe    = θ[3]
@@ -297,7 +312,7 @@ function rmat!(mx, θ, rz,  ::SPGAUD_)
             end
         end
     end
-    mx
+    return mx
 end
 
 #UN
@@ -318,14 +333,14 @@ function unrmat(θ::AbstractVector{T}, rz) where T
     @inbounds @simd for m = 1:rm
         mx[m, m] *= mx[m, m]
     end
-    Symmetric(mx)
+    return Symmetric(mx)
 end
-function rmat!(mx, θ, rz::AbstractMatrix, ::UN_)
+function rmat!(mx, θ, rz::AbstractMatrix, ::UN_, ::Int)
     vec    = tmul_unsafe(rz, θ)
     rm     = size(mx, 1)
     rcov  = unrmat(θ, rz)
     mulαβαtinc!(mx, rz, rcov)
-    mx
+    return mx
 end
 ###############################################################################
 ###############################################################################
diff --git a/src/statsbase.jl b/src/statsbase.jl
index 888be62..a29b318 100644
--- a/src/statsbase.jl
+++ b/src/statsbase.jl
@@ -75,7 +75,7 @@ Model coefficients (β).
 StatsBase.coef(lmm::LMM) = copy(coef_(lmm))
 
 function coef_(lmm::LMM)
-    lmm.result.beta
+    return lmm.result.beta
 end
 """
     StatsBase.coefnames(lmm::LMM) = StatsBase.coefnames(lmm.mf)
@@ -99,7 +99,7 @@ end
 DOF residuals: N - rank(X), where N - total number of observations.
 """
 function StatsBase.dof_residual(lmm::LMM)
-    nobs(lmm) - lmm.rankx
+    return nobs(lmm) - lmm.rankx
 end
 
 """
@@ -108,7 +108,7 @@ end
 DOF.
 """
 function StatsBase.dof(lmm::LMM)
-    lmm.nfixed + lmm.covstr.tl
+    return lmm.nfixed + lmm.covstr.tl
 end
 
 """
@@ -117,7 +117,7 @@ end
 Return loglikelihood value.
 """
 function StatsBase.loglikelihood(lmm::LMM)
-    -lmm.result.reml/2
+    return -lmm.result.reml/2
 end
 
 """
@@ -128,7 +128,7 @@ Akaike Information Criterion.
 function StatsBase.aic(lmm::LMM)
     l = loglikelihood(lmm)
     d = lmm.covstr.tl
-    -2l + 2d
+    return -2l + 2d
 end
 
 """
@@ -140,7 +140,7 @@ function StatsBase.bic(lmm::LMM)
     l = loglikelihood(lmm)
     d = lmm.covstr.tl
     n = nobs(lmm) - lmm.nfixed
-    -2l + d * log(n)
+    return -2l + d * log(n)
 end
 
 """
@@ -152,7 +152,7 @@ function StatsBase.aicc(lmm::LMM)
     l = loglikelihood(lmm)
     d = lmm.covstr.tl
     n = nobs(lmm) - lmm.nfixed
-    -2l + (2d * n) / (n - d - 1.0)
+    return -2l + (2d * n) / (n - d - 1.0)
 end
 
 """
@@ -164,14 +164,14 @@ function caic(lmm::LMM)
     l = loglikelihood(lmm)
     d = lmm.covstr.tl
     n = nobs(lmm) - lmm.nfixed
-    -2l + d * (log(n) + 1.0)
+    return -2l + d * (log(n) + 1.0)
 end
 
 """
     StatsBase.isfitted(lmm::LMM)
 """
 function StatsBase.isfitted(lmm::LMM)
-    lmm.result.fit
+    return lmm.result.fit
 end
 """
     StatsBase.vcov(lmm::LMM)
@@ -187,7 +187,7 @@ Standard error
 StatsBase.stderror(lmm::LMM) = copy(stderror_(lmm))
 
 function stderror_(lmm::LMM)
-    lmm.result.se
+    return lmm.result.se
 end
 
 stderror!(v, lmm::LMM) = copyto!(v, lmm.result.se)
@@ -239,7 +239,7 @@ end
 Responce varible name.
 """
 function StatsBase.responsename(lmm::LMM)
-    StatsBase.coefnames(lmm.f)[1]
+    return StatsBase.coefnames(lmm.f)[1]
 end
 
 
diff --git a/src/sweep.jl b/src/sweep.jl
index 3610621..744e21a 100644
--- a/src/sweep.jl
+++ b/src/sweep.jl
@@ -10,14 +10,14 @@ function nsyrk!(α, x, A)
             @inbounds A[i, j] += x[i] * xjα
         end
     end
-    A
+    return A
 end
 function nsyrk!(α, x, A::AbstractArray{T}) where T <: AbstractFloat
-    BLAS.syrk!('U', 'N', α, x, one(T), A)
+    return BLAS.syrk!('U', 'N', α, x, one(T), A)
 end
 
 function sweep!(A::AbstractArray{T}, k::Integer, inv::Bool = false) where T
-    sweepb!(Vector{T}(undef, size(A, 2)), A, k, inv)
+    return sweepb!(Vector{T}(undef, size(A, 2)), A, k, inv)
 end
 function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, k::Integer, inv::Bool = false) where T <: Number
     p = checksquare(A)
@@ -26,7 +26,7 @@ function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, k::Integer, i
     @simd for j in 1:k
         @inbounds akk[j] = A[j, k]
     end
-    @simd for j in (k+1):p
+    @simd for j in (k + 1):p
         @inbounds akk[j] = A[k, j]
     end
     # syrk!(uplo, trans, alpha, A, beta, C)
@@ -44,11 +44,11 @@ function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, k::Integer, i
         @inbounds A[k, j] = akk[j]
     end
     @inbounds A[k, k] = -d
-    A
+    return A
 end
 function sweep!(A::AbstractArray{T, 2}, ks::AbstractVector{I}, inv::Bool = false; logdet::Bool = false) where {T <: Number, I <: Integer}
     akk = Vector{T}(undef, size(A,2))
-    sweepb!(akk, A, ks, inv; logdet = logdet)
+    return sweepb!(akk, A, ks, inv; logdet = logdet)
 end
 function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, ks::AbstractVector{I}, inv::Bool = false; logdet::Bool = false) where
         {T <: Number, I<:Integer}
@@ -57,7 +57,7 @@ function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, ks::AbstractV
     if logdet
         ld = 0
         for k in ks
-            @inbounds Akk = A[k,k]
+            @inbounds Akk = A[k, k]
             if Akk > 0
                 ld += log(Akk)
             else
@@ -75,5 +75,5 @@ function sweepb!(akk::AbstractArray{T, 1}, A::AbstractArray{T, 2}, ks::AbstractV
             sweepb!(akk, A, k, inv)
         end
     end
-    A, ld, noerror
+    return A, ld, noerror
 end
diff --git a/src/typeiii.jl b/src/typeiii.jl
index 7128330..c7e2345 100644
--- a/src/typeiii.jl
+++ b/src/typeiii.jl
@@ -84,7 +84,7 @@ function contrast(lmm, l::AbstractMatrix; name::String = "Contrast", ddf = :satt
     else
         pval = NaN
     end
-    ContrastTable([name], [F], [ndf], [df], [pval])
+    return ContrastTable([name], [F], [ndf], [df], [pval])
 end
 
 function Base.show(io::IO, at::ContrastTable)
diff --git a/src/utils.jl b/src/utils.jl
index 3a8efad..b714afa 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -7,7 +7,7 @@ function initvar(y::Vector, X::Matrix{T}) where T
     r    = copy(y)
     LinearAlgebra.BLAS.gemv!('N', one(T), X, β, -one(T), r)
     #r    = y .- X * β
-    dot(r, r)/(length(r) - size(X, 2)), β
+    return dot(r, r)/(length(r) - size(X, 2)), β
 end
 ################################################################################
 function fixedeffn(f::FormulaTerm)
@@ -16,10 +16,7 @@ end
 function fixedeffn(lmm::LMM)
     fixedeffn(lmm.f) 
 end
-#=
-function nterms(mf::ModelFrame)
-    mf.schema.schema.count
-=#
+
 function nterms(rhs::Union{Tuple{Vararg{AbstractTerm}}, Nothing, AbstractTerm})
     if isa(rhs, Term)
         p = 1
@@ -28,7 +25,7 @@ function nterms(rhs::Union{Tuple{Vararg{AbstractTerm}}, Nothing, AbstractTerm})
     else
         p = 0
     end
-    p
+    return p
 end
 """
     Term name.
@@ -71,7 +68,7 @@ function lcontrast(lmm::LMM, i::Int)
             mx[j, inds[j]] = 1
         end
     end
-    mx
+    return mx
 end
 ################################################################################
 #                        VAR LINK
@@ -79,17 +76,17 @@ end
 
 function vlink(σ::T) where T <: Real
     if σ < -21.0 return one(T)*7.582560427911907e-10 end #Experimental
-    exp(σ)
+    return exp(σ)
 end
 function vlinkr(σ::T) where T <: Real
-    log(σ)
+    return log(σ)
 end
 
 function vlinksq(σ::T) where T <: Real
-    σ*σ
+    return σ*σ
 end
 function vlinksqr(σ::T) where T <: Real
-    sqrt(σ)
+    return sqrt(σ)
 end
 
 function rholinkpsigmoid(ρ::T) where T <: Real
@@ -144,7 +141,7 @@ function varlinkvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm)
             end
         end
     end
-    v
+    return v
 end
 function varlinkrvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm)
     @inbounds for i = 1:length(v)
@@ -168,7 +165,7 @@ function varlinkrvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm)
             end
         end
     end
-    v
+    return v
 end
 function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm)
     s = similar(v)
@@ -195,22 +192,22 @@ function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm)
             s[i] = v[i]
         end
     end
-    s
+    return s
 end
 ################################################################################
 function m2logreml(lmm)
-    lmm.result.reml
+    return lmm.result.reml
 end
 function logreml(lmm)
-    -m2logreml(lmm)/2
+    return -m2logreml(lmm)/2
 end
 function m2logreml(lmm, theta; maxthreads::Int = num_cores())
-    reml_sweep_β(lmm, LMMDataViews(lmm), theta; maxthreads = maxthreads)[1]
+    return reml_sweep_β(lmm, LMMDataViews(lmm), theta; maxthreads = maxthreads)[1]
 end
 ################################################################################
 
 function optim_callback(os)
-    false
+    return false
 end
 ################################################################################
 """
@@ -219,7 +216,7 @@ end
 Return true if CovarianceType is ZERO.
 """
 function zeroeff(eff)
-    isa(eff.covtype.s, ZERO)
+    return isa(eff.covtype.s, ZERO)
 end
 """
     raneffn(lmm)
@@ -229,7 +226,7 @@ function raneffn(lmm)
     if zeroeff(lmm.covstr.random[1]) 
         return 0
     end
-    length(lmm.covstr.random)
+    return length(lmm.covstr.random)
 end
 
 """
@@ -240,7 +237,7 @@ function gmatrix(lmm::LMM{T}, r::Int) where T
     if r > length(lmm.covstr.random) error("Invalid random effect number: $(r)!") end
     G = zeros(T, lmm.covstr.q[r], lmm.covstr.q[r])
     gmat!(G, view(lmm.result.theta, lmm.covstr.tr[r]), lmm.covstr.random[r].covtype.s)
-    Symmetric(G)
+    return Symmetric(G)
 end
 
 
@@ -250,7 +247,7 @@ end
 Return true if all variance-covariance matrix (G) of random effect is positive definite.
 """
 function gmatrixipd(lmm::LMM)
-    lmm.result.ipd
+    return lmm.result.ipd
 end
 
 """
@@ -263,7 +260,7 @@ function rmatrix(lmm::LMM{T}, i::Int) where T
     R    = zeros(T, q, q)
     rθ   = lmm.covstr.tr[lmm.covstr.rn + 1:end]
     rmat_base_inc!(R, lmm.result.theta, rθ, lmm.covstr, i)
-    Symmetric(R)
+    return Symmetric(R)
 end
 
 #####################################################################
@@ -272,10 +269,14 @@ function applywts!(::Any, ::Int, ::Nothing)
     nothing
 end
 
-function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts)
-    mulβdαβd!(V, wts.sqrtwts[i])
+function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts{<:Vector})
+    return mulβdαβd!(V, wts.sqrtwts[i])
 end
 
+function applywts!(V::AbstractMatrix, i::Int, wts::LMMWts{<:Matrix})
+    V .*= wts.sqrtwts[i]
+    return V
+end
 
 #####################################################################
 
@@ -290,15 +291,14 @@ function vmatrix!(V, θ, lmm::LMM, i::Int) # pub API
     rθ   = lmm.covstr.tr[lmm.covstr.rn + 1:end]
     rmat_base_inc!(V, θ, rθ, lmm.covstr, i) # Repeated vector
     applywts!(V, i, lmm.wts) 
-    zgz_base_inc!(V, gvec, lmm.covstr, i)
-    
+    return zgz_base_inc!(V, gvec, lmm.covstr, i) 
 end
 
 # !!! Main function REML used
 @noinline function vmatrix!(V, G, θ, rθ, lmm::LMM, i::Int)
     rmat_base_inc!(V, θ, rθ, lmm.covstr, i)  # Repeated vector
     applywts!(V, i, lmm.wts) 
-    zgz_base_inc!(V, G, lmm.covstr, i)
+    return zgz_base_inc!(V, G, lmm.covstr, i)
 end
 
 """
@@ -307,7 +307,7 @@ end
 Return variance-covariance matrix V for i bolock.
 """
 function vmatrix(lmm::LMM, i::Int)
-    vmatrix(lmm.result.theta, lmm, i)
+    return vmatrix(lmm.result.theta, lmm, i)
 end
 
 function vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int) where T
@@ -315,7 +315,7 @@ function vmatrix(θ::AbstractVector{T}, lmm::LMM, i::Int) where T
     gvec = gmatvec(θ, lmm.covstr)
     rθ   = lmm.covstr.tr[lmm.covstr.rn + 1:end]
     vmatrix!(V, gvec, θ, rθ, lmm, i) # Repeated vector
-    Symmetric(V)
+    return Symmetric(V)
 end
 # For Multiple Imputation
 function vmatrix(θ::Vector, covstr::CovStructure, lmmwts, i::Int)
@@ -325,11 +325,11 @@ function vmatrix(θ::Vector, covstr::CovStructure, lmmwts, i::Int)
     rmat_base_inc!(V, θ, rθ, covstr, i)  # Repeated vector
     applywts!(V, i, lmmwts) 
     zgz_base_inc!(V, gvec, covstr, i)
-    Symmetric(V)
+    return Symmetric(V)
 end
 
 function blockgmatrix(lmm::LMM{T}) where T
-    blockgmatrix(lmm, (1, 1))
+    return blockgmatrix(lmm, (1, 1))
 end
 
 function blockgmatrix(lmm::LMM{T}, v) where T
@@ -348,7 +348,7 @@ function blockgmatrix(lmm::LMM{T}, v) where T
             s = e + 1
         end
     end
-    G
+    return G
 end
 
 function blockzmatrix(lmm::LMM{T}, i) where T
@@ -368,7 +368,7 @@ function blockzmatrix(lmm::LMM{T}, i) where T
         s  = e + 1
         mx[j] = smx
     end
-   hcat(mx...)
+    return hcat(mx...)
 end
 """
     raneff(lmm::LMM{T}, i)
@@ -399,8 +399,7 @@ function raneff(lmm::LMM{T}, block) where T
             rvsbj[i][j] =  subjname => rv[s:e]
         end
     end
-    rvsbj
-    
+    return rvsbj  
 end
 """
     raneff(lmm::LMM{T})
@@ -418,7 +417,7 @@ function raneff(lmm::LMM)
             end
         end
     end
-    fb
+    return fb
 end
 
 """
@@ -438,19 +437,19 @@ function hessian(lmm, theta)
     vloptf(x) = reml_sweep_β(lmm, lmm.dv, x, lmm.result.beta)[1]
     chunk  = ForwardDiff.Chunk{min(8, length(theta))}()
     hcfg   = ForwardDiff.HessianConfig(vloptf, theta, chunk)
-    ForwardDiff.hessian(vloptf, theta, hcfg)
+    return ForwardDiff.hessian(vloptf, theta, hcfg)
 end
 function hessian(lmm)
     if !lmm.result.fit error("Model not fitted!") end
-    hessian(lmm, lmm.result.theta)
+    return hessian(lmm, lmm.result.theta)
 end
 ################################################################################
 ################################################################################
 function StatsModels.termvars(ve::VarEffect)
-    termvars(ve.formula)
+    return termvars(ve.formula)
 end
 function StatsModels.termvars(ve::Vector{VarEffect})
-    union(termvars.(ve)...)
+    return union(termvars.(ve)...)
 end
 
 ################################################################################
diff --git a/src/varstruct.jl b/src/varstruct.jl
index 49a426a..b2bb594 100644
--- a/src/varstruct.jl
+++ b/src/varstruct.jl
@@ -16,9 +16,7 @@ function StatsModels.ContrastsMatrix(contrasts::RawCoding, levels::AbstractVecto
                              contrasts)
 end
 function StatsModels.modelcols(t::CategoricalTerm{RawCoding, T, N}, d::NamedTuple) where T where N
-    #v = d[t.sym]
-    #reshape(v, length(v), 1)  
-    d[t.sym]
+    return d[t.sym]
 end
 
 ################################################################################
@@ -42,7 +40,7 @@ end
 function modelparse(term::FunctionTerm{typeof(|)})
     eff, subj = term.args
     if !isa(subj, AbstractTerm) || isa(subj, FunctionTerm{typeof(*), Vector{Term}}) throw(FormulaException("Subject term type not <: AbstractTerm. Use `term` or `interaction term` only. Maybe you are using something like this: `@covstr(factor|term1*term2)` or `@covstr(factor|(term1+term2))`. Use only `@covstr(factor|term)` or `@covstr(factor|term1&term2)`.")) end
-    eff, subj
+    return eff, subj
 end
 function modelparse(term)
     throw(FormulaException("Model term type not <: FunctionTerm{typeof(|)}. Use model like this: `@covstr(factor|subject)`. Maybe you are using something like this: `@covstr(factor|term1+term2)`. Use only `@covstr(factor|term)` or `@covstr(factor|term1&term2)`."))
@@ -140,7 +138,7 @@ function sabjcrossdicts(d1, d2)
             end
         end
     end
-    v
+    return v
 end
 
 tabcols(data, symbs) = Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs)
@@ -170,7 +168,7 @@ function raneflenv(covstr, block)
     for i = 1:l
         v[i] = length(covstr.esb.sblock[block, i])
     end
-    v
+    return v
 end
 """
     Covarince structure.
@@ -255,7 +253,6 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
             end
         end
         # RANDOM EFFECTS
-        #if random[1].covtype.z #IF NOT ZERO
             @inbounds for i = 1:rn
                 if length(random[i].coding) == 0
                     fill_coding_dict!(random[i].model, random[i].coding, data)
@@ -283,7 +280,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
                 fillur!(tr, i, t)
                 symbs       = StatsModels.termvars(random[i].subj)
                 if length(symbs) > 0
-                    cdata     = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data_), x) for x in symbs)
+                    cdata     = tabcols(data, symbs) 
                     dicts[i]  = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
                     indsdict!(dicts[i], cdata)
                 else
@@ -310,11 +307,10 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
             end
             
             schema[rn + i] = apply_schema(repeated[i].model, StatsModels.schema(data_, repeated[i].coding))
-            #rz_[i]       = reduce(hcat, modelcols(schema[rn+i], data))
             rz_[i]       = modelcols(MatrixTerm(schema[rn+i]), data_)
             symbs        = StatsModels.termvars(repeated[i].subj)
             if length(symbs) > 0
-                cdata    = tabcols(data, symbs) # Tuple(Tables.getcolumn(Tables.columns(data), x) for x in symbs)
+                cdata    = tabcols(data, symbs) 
                 dicts[rn + i]  = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
                 indsdict!(dicts[rn + i], cdata)
             else
@@ -336,9 +332,6 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
         # Theta length
         tl  = sum(t)
         ########################################################################
-        #if any(x-> 1 in keys(x), dicts[1:end-1])
-        #    blocks = [first(dicts)[1]]
-        #else
             if random[1].covtype.z  # if first random effect not null
                 subjblockdict = dicts[1]
                 if length(dicts) > 2 # if more than 2 random effects
@@ -363,9 +356,7 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
                 dicts[rn+i] = subjblockdict
             end
 
-
             blocks  = collect(values(subjblockdict))
-        #end
 
         sblock = Matrix{Vector{Tuple{Vector{Int}, Int}}}(undef, length(blocks), alleffl)
         nblock = []
@@ -394,6 +385,13 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
         end
         esb = EffectSubjectBlock(sblock, nblock)
         #######################################################################
+        # Postprocessing
+        # Modify repeated effect covariance type for some types
+        # Maybe it will be removed
+        for r in repeated 
+            applycovschema!(r.covtype.s, blocks)
+        end
+        #######################################################################
         new{eltype(z), T2}(random, repeated, schema, rcnames, blocks, rn, rtn, rpn, z, esb, zrndur, rz, q, t, tr, tl, ct, emap, sn, maxn)
     end
 end
@@ -422,13 +420,13 @@ end
 ################################################################################
 
 function fill_coding_dict!(t::T, d::Dict, data) where T <: Union{ConstantTerm, InterceptTerm, FunctionTerm}
-    d
+    return d
 end
 function fill_coding_dict!(t::T, d::Dict, data) where T <: Term
     if typeof(Tables.getcolumn(data, t.sym)) <: AbstractCategoricalVector || !(typeof(Tables.getcolumn(data, t.sym)) <: AbstractVector{V} where V <: Real)
         d[t.sym] = StatsModels.FullDummyCoding()
     end
-    d
+    return d
 end
 #=
 function fill_coding_dict!(t::T, d::Dict, data) where T <: InteractionTerm
@@ -450,7 +448,7 @@ function fill_coding_dict_ct!(t, d, data)
             fill_coding_dict!(i, d, data)
         end
     end
-    d
+    return d
 end
 #=
 function fill_coding_dict!(t::T, d::Dict, data) where T <: Tuple{Vararg{AbstractTerm}}
diff --git a/src/vartypes.jl b/src/vartypes.jl
index 1b4707b..d8bf870 100644
--- a/src/vartypes.jl
+++ b/src/vartypes.jl
@@ -5,6 +5,10 @@
 ################################################################################
 
 struct SI_ <: AbstractCovarianceType end
+mutable struct SWC_{W<:AbstractMatrix, B<:Vector{<:AbstractMatrix}} <: AbstractCovarianceType 
+    wtsm::W
+    wtsb::B
+end
 struct DIAG_ <: AbstractCovarianceType end
 struct AR_ <: AbstractCovarianceType end
 struct ARH_ <: AbstractCovarianceType end
@@ -64,6 +68,45 @@ function ScaledIdentity()
     CovarianceType(SI_())
 end
 const SI = ScaledIdentity()
+
+# Experimental
+"""
+    ScaledWeightedCov(wtsm::AbstractMatrix{T})
+
+!!! warning
+    Experimental
+
+Scaled weighted covariance matrix, where `wtsm` - `NxN` within block correlation matrix (N - total number of observations). 
+Used only for repeated effect. 
+
+SWC = ScaledWeightedCov
+
+```math
+R = Corr(W) * \\sigma_c^2
+```
+
+where ``Corr(W)`` - diagonal correlation matrix. 
+
+example:
+
+```julia
+matwts = Symmetric(UnitUpperTriangular(rand(size(df0,1), size(df0,1))))
+lmm = LMM(@formula(var~sequence+period+formulation), df0;
+    repeated = VarEffect(@covstr(1|subject), SWC(matwts)))
+fit!(lmm)
+
+```
+!!! note
+
+There is no `wtsm` checks for symmetricity or values.
+
+"""
+function ScaledWeightedCov(wtsm::AbstractMatrix{T}) where T
+    wtsb = Matrix{T}[]
+    CovarianceType(SWC_(wtsm, wtsb))
+end
+const SWC = ScaledWeightedCov
+
 """
     Diag()
 
@@ -362,6 +405,9 @@ end
 function covstrparam(ct::SI_, ::Int)::Tuple{Int, Int}
     return (1, 0)
 end
+function covstrparam(ct::SWC_, ::Int)::Tuple{Int, Int}
+    return (1, 0)
+end
 function covstrparam(ct::DIAG_, t::Int, )::Tuple{Int, Int}
     return (t, 0)
 end
@@ -412,6 +458,9 @@ end
 function rcoefnames(s, t, ct::SI_)
     return ["σ² "]
 end
+function rcoefnames(s, t, ct::SWC_)
+    return ["σ² "]
+end
 function rcoefnames(s, t, ct::DIAG_)
     if isa(coefnames(s), AbstractArray{T,1} where T) l = length(coefnames(s)) else l = 1 end
     return fill!(Vector{String}(undef, l), "σ² ") .* string.(coefnames(s))
@@ -483,7 +532,7 @@ function indfromtn(ind, s)
             break
         end
     end
-    m, s + ind - b
+    return m, s + ind - b
 end
 
 function rcoefnames(s, t, ct::UN_)
@@ -510,6 +559,21 @@ function rcoefnames(s, t, ct::AbstractCovarianceType)
     v .= "Val "
     return v
 end
+################################################################################
+# APPLY COV SCHEMA
+################################################################################
+function applycovschema!(::AbstractCovarianceType, ::Any)
+    nothing
+end
+
+function applycovschema!(ct::SWC_{<:AbstractMatrix{T}}, vcovblock) where T
+    if length(ct.wtsb) == 0
+        for i in eachindex(vcovblock)
+            push!(ct.wtsb, ct.wtsm[vcovblock[i], vcovblock[i]])
+        end
+    end
+    return ct
+end
 
 ################################################################################
 # SHOW
@@ -524,6 +588,9 @@ end
 function Base.show(io::IO, ct::SI_)
     print(io, "SI")
 end
+function Base.show(io::IO, ct::SWC_)
+    print(io, "SWC")
+end
 function Base.show(io::IO, ct::DIAG_)
     print(io, "DIAG")
 end
diff --git a/test/devtest.jl b/test/devtest.jl
index d70c114..db017c1 100644
--- a/test/devtest.jl
+++ b/test/devtest.jl
@@ -1,9 +1,8 @@
 #using NLopt
 using Metida
 
-using DataFrames, CSV, StatsModels, LinearAlgebra, ForwardDiff, ForwardDiff, Optim, Distributions, CategoricalArrays
-#using SnoopCompile
-#using LineSearches
+using DataFrames, CSV, StatsModels, LinearAlgebra, CategoricalArrays, Dates
+
 using BenchmarkTools
 path    = dirname(@__FILE__)
 cd(path)
@@ -13,349 +12,95 @@ ftdf         = CSV.File(path*"/csv/1fptime.csv"; types = [String, String, Float6
 ftdf2        = CSV.File(path*"/csv/1freparma.csv"; types = [String, String, Float64, Float64]) |> DataFrame
 ftdf3        = CSV.File(path*"/csv/ftdf3.csv"; types =
 [String,  Float64, Float64, String, String, String, String, String, Float64]) |> DataFrame
+hdp         = CSV.File("hdp.csv") |> DataFrame
+transform!(hdp, :DID => categorical);
+transform!(hdp, :HID=> categorical);
+transform!(hdp, :Sex=> categorical);
+transform!(hdp, :School=> categorical);
+transform!(hdp, :pain=> categorical);
 pkgversion(m::Module) = Pkg.TOML.parsefile(joinpath(dirname(string(first(methods(m.eval)).file)), "..", "Project.toml"))["version"]
-# MODEL 1
+
+
+results = DataFrame(datetime =[], model = [], mintime =[], memory = [], allocs = [])
+b = Vector{Any}(undef, 4)
 ################################################################################
 # Metida
 ################################################################################
-#nt = LinearAlgebra.BLAS.get_num_threads()
-#LinearAlgebra.BLAS.set_num_threads(16)
-#LinearAlgebra.BLAS.set_num_threads(nt)
-
+# MODEL 1
+println("MODEL 1")
 lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
 random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
 )
-b11 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-
-#=
-lmm = Metida.LMM(@formula(response ~1 + factor2), ftdf3;
-repeated = Metida.VarEffect(Metida.@covstr(p|subject), Metida.CSH),
-)
-@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-=#
-#:LN_BOBYQA :LN_NEWUOA 
-#@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16, solver = :nlopt, optmethod = :LN_NEWUOA) seconds = 15
-#@benchmark Metida.fit!($lmm, optmethod = Metida.LBFGS_OM, hes = false; maxthreads = 16) seconds = 15
-#@benchmark Metida.fit!($lmm, optmethod = Metida.BFGS_OM, hes = false; maxthreads = 16) seconds = 15
-#@benchmark Metida.fit!($lmm, optmethod = Metida.CG_OM, hes = false; maxthreads = 16) seconds = 15
-#@benchmark Metida.fit!($lmm, optmethod = Optim.NelderMead(), hes = false; maxthreads = 16) seconds = 15
-
-
-#@time Metida.fit!(lmm, hes = false)
-
-
-
-################################################################################
-# MetidaNLopt
-################################################################################
-using MetidaNLopt
-b12 = @benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15
-
-################################################################################
-# MetidaCu
-################################################################################
-using MetidaCu
-b13 = @benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15
-
-
+b[1] = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
 # MODEL 2
-
+println("MODEL 2")
 lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
 random = Metida.VarEffect(Metida.@covstr(1 + time|factor), Metida.ARH),
 )
-b21 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-b22 = @benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15
-b23 = @benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15
-
+b[2]  = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
 # MODEL 3
-
+println("MODEL 3")
 lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
 random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
+repeated = Metida.VarEffect(Metida.@covstr(1|subject), Metida.AR),
 )
-
-b31 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-b32 = @benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15
-b33 = @benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15
-
+b[3]  = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
 # MODEL 4
-
-hdp         = CSV.File("hdp.csv") |> DataFrame
-transform!(hdp, :DID => categorical);
-transform!(hdp, :HID=> categorical);
-transform!(hdp, :Sex=> categorical);
-transform!(hdp, :School=> categorical);
-transform!(hdp, :pain=> categorical);
-################################################################################
-# Metida
-################################################################################
-
+println("MODEL 4")
 lmm = Metida.LMM(@formula(tumorsize ~ 1 + CancerStage), hdp;
 random = Metida.VarEffect(Metida.@covstr(1|HID), Metida.DIAG),
 )
+b[4]  = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
 
-b41 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-b42 = @benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15
-b43 = @benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15
-
-# MODEL 5 maximum 1437 observation-per-subject (10 subjects)
-lmm = Metida.LMM(@formula(tumorsize ~ 1 + CancerStage), hdp;
-random = Metida.VarEffect(Metida.@covstr(1|ntumors), Metida.SI),
-)
-
-#b51 = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-b52 = @benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15
-b53 = @benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15
-
-# MODEL 6: maximum 3409 observation-per-subject (4 subjects)
-
-
-println("Metida version: ", pkgversion(Metida))
-println("MetidaNLopt version: ", pkgversion(MetidaNLopt))
-println("MetidaCu version: ", pkgversion(MetidaCu))
-
-println("MODEL 1")
-println("# Metida")
-display(b11)
-println("# MetidaNLopt")
-display(b12)
-println("# MetidaCu")
-display(b13)
-println()
-println()
-
-println("MODEL 2")
-println("# Metida")
-display(b21)
-println("# MetidaNLopt")
-display(b22)
-println("# MetidaCu")
-display(b23)
-println()
-println()
-
-println("MODEL 3")
-println("# Metida")
-display(b31)
-println("# MetidaNLopt")
-display(b32)
-println("# MetidaCu")
-display(b33)
-println()
-println()
-
-println("MODEL 4")
-println("# Metida")
-display(b41)
-println("# MetidaNLopt")
-display(b42)
-println("# MetidaCu")
-display(b43)
-println()
-println()
-
-println("MODEL 5")
-#println("# Metida")
-#display(b51)
-println("# MetidaNLopt")
-display(b52)
-println("# MetidaCu")
-display(b53)
-println()
-println()
+for i = 1:4
+  println("MODEL $i")
+  display(b[i])
+  push!(results, (now(), "Model $i", minimum(b[i]).time, minimum(b[i]).memory, minimum(b[i]).allocs))
+end
 
-#Julia 1.6.3
 #=
-julia>
-Metida version: 0.12.0
-MetidaNLopt version: 0.4.0
-MetidaCu version: 0.4.1
 MODEL 1
-# Metida
-BenchmarkTools.Trial: 1184 samples with 1 evaluation.
- Range (min … max):   5.489 ms … 161.314 ms  ┊ GC (min … max):  0.00% … 94.86%
- Time  (median):      8.535 ms               ┊ GC (median):     0.00%
- Time  (mean ± σ):   12.700 ms ±  22.994 ms  ┊ GC (mean ± σ):  33.38% ± 16.90%
-
-  ▇█▅
-  ███▄▄▄▁▁▁▁▁▁▁▁▁▁▁▄▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▄▄▁▅▄▆▁▄▄▄▄▅▄▆▄ █
-  5.49 ms       Histogram: log(frequency) by time       144 ms <
-
- Memory estimate: 22.63 MiB, allocs estimate: 37225.
-# MetidaNLopt
-BenchmarkTools.Trial: 173 samples with 1 evaluation.
- Range (min … max):  74.294 ms … 186.791 ms  ┊ GC (min … max): 0.00% … 58.09%
- Time  (median):     78.964 ms               ┊ GC (median):    0.00%
- Time  (mean ± σ):   86.794 ms ±  22.516 ms  ┊ GC (mean ± σ):  9.19% ± 14.57%
-
-    █▃
-  ▅▆██▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅█▁▁▁▁▁▁▄▁▁▄▁▄▁▄▁▁▁▁▁▁▁▄▁▁▁▄▁▄▁▄▄▁▁▄▁▄ ▄
-  74.3 ms       Histogram: log(frequency) by time       176 ms <
-
- Memory estimate: 56.61 MiB, allocs estimate: 481209.
-# MetidaCu
-BenchmarkTools.Trial: 2 samples with 1 evaluation.
- Range (min … max):  9.128 s …    9.323 s  ┊ GC (min … max): 0.00% … 0.58%
- Time  (median):     9.226 s               ┊ GC (median):    0.29%
- Time  (mean ± σ):   9.226 s ± 137.689 ms  ┊ GC (mean ± σ):  0.29% ± 0.41%
-
-  █                                                        █
-  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  9.13 s         Histogram: frequency by time         9.32 s <
-
- Memory estimate: 143.46 MiB, allocs estimate: 2524366.
+BenchmarkTools.Trial: 675 samples with 1 evaluation.
+ Range (min … max):   8.194 ms … 798.405 ms  ┊ GC (min … max):  0.00% … 98.64%
+ Time  (median):     10.882 ms               ┊ GC (median):     0.00%
+ Time  (mean ± σ):   22.823 ms ±  81.034 ms  ┊ GC (mean ± σ):  52.70% ± 14.48%
 
+  █  
+  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▄▃ ▆
+  8.19 ms       Histogram: log(frequency) by time       533 ms <
 
+ Memory estimate: 19.38 MiB, allocs estimate: 31265.
 MODEL 2
-# Metida
-BenchmarkTools.Trial: 18 samples with 1 evaluation.
- Range (min … max):  829.008 ms … 850.212 ms  ┊ GC (min … max): 0.00% … 1.35%
- Time  (median):     839.290 ms               ┊ GC (median):    0.58%
- Time  (mean ± σ):   839.688 ms ±   7.197 ms  ┊ GC (mean ± σ):  0.62% ± 0.63%
-
-  ▁   █    ▁   ▁▁  ▁   █               ▁       ▁█▁▁      ▁  ▁ ▁
-  █▁▁▁█▁▁▁▁█▁▁▁██▁▁█▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁████▁▁▁▁▁▁█▁▁█▁█ ▁
-  829 ms           Histogram: frequency by time          850 ms <
-
- Memory estimate: 140.68 MiB, allocs estimate: 7919.
-# MetidaNLopt
-BenchmarkTools.Trial: 135 samples with 1 evaluation.
- Range (min … max):  108.482 ms … 119.274 ms  ┊ GC (min … max): 0.00% … 4.83%
- Time  (median):     110.394 ms               ┊ GC (median):    0.00%
- Time  (mean ± σ):   111.933 ms ±   2.719 ms  ┊ GC (mean ± σ):  1.80% ± 2.26%
-
-        ▂ ▂ █▂▄▃                                 ▂▂
-  █▅▆▃▅▅███▇█████▆▅▁▃▁▃▃▃▁▁▁▁▃▁▁▁▁▁▃▁▁▁▁▃▃▃▅▁▆▃▇▇██▇▇▇▅▅▁▆▁▁▁▃▃ ▃
-  108 ms           Histogram: frequency by time          117 ms <
-
- Memory estimate: 91.48 MiB, allocs estimate: 33646.
-# MetidaCu
-BenchmarkTools.Trial: 29 samples with 1 evaluation.
- Range (min … max):  507.681 ms … 536.903 ms  ┊ GC (min … max): 0.00% … 1.58%
- Time  (median):     518.749 ms               ┊ GC (median):    0.00%
- Time  (mean ± σ):   521.751 ms ±  10.164 ms  ┊ GC (mean ± σ):  0.58% ± 0.74%
-
-           █             ▃                          ▃     ▃   ▃
-  ▇▁▁▁▇▇▁▇▁█▇▇▇▇▁▇▁▁▁▇▇▁▁█▁▇▁▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▇▁█▁▇▇▁▇█▇▁▁█ ▁
-  508 ms           Histogram: frequency by time          537 ms <
-
- Memory estimate: 93.94 MiB, allocs estimate: 95784.
+BenchmarkTools.Trial: 23 samples with 1 evaluation.
+ Range (min … max):  657.833 ms … 704.073 ms  ┊ GC (min … max): 0.00% … 3.23%
+ Time  (median):     675.148 ms               ┊ GC (median):    0.00%        
+ Time  (mean ± σ):   677.822 ms ±  12.920 ms  ┊ GC (mean ± σ):  1.41% ± 1.62%
 
+  ▁     ▁▁ ▁ ██  ▁ ▁   ▁▁   ▁  ▁  ▁   ▁   █    ▁ █   ▁        ▁  
+  █▁▁▁▁▁██▁█▁██▁▁█▁█▁▁▁██▁▁▁█▁▁█▁▁█▁▁▁█▁▁▁█▁▁▁▁█▁█▁▁▁█▁▁▁▁▁▁▁▁█ ▁
+  658 ms           Histogram: frequency by time          704 ms <
 
+ Memory estimate: 132.83 MiB, allocs estimate: 5085.
 MODEL 3
-# Metida
-BenchmarkTools.Trial: 1100 samples with 1 evaluation.
- Range (min … max):   5.554 ms … 213.611 ms  ┊ GC (min … max):  0.00% … 95.95%
- Time  (median):      8.962 ms               ┊ GC (median):     0.00%
- Time  (mean ± σ):   13.684 ms ±  28.997 ms  ┊ GC (mean ± σ):  37.50% ± 16.47%
-
-  ▇█
-  ██▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▆▄▁▆▄▄▁▄▄▄▆ ▇
-  5.55 ms       Histogram: log(frequency) by time       181 ms <
-
- Memory estimate: 22.63 MiB, allocs estimate: 37224.
-# MetidaNLopt
-BenchmarkTools.Trial: 174 samples with 1 evaluation.
- Range (min … max):  75.122 ms … 186.340 ms  ┊ GC (min … max): 0.00% … 57.72%
- Time  (median):     79.034 ms               ┊ GC (median):    0.00%
- Time  (mean ± σ):   86.617 ms ±  22.856 ms  ┊ GC (mean ± σ):  8.85% ± 14.37%
-
-   ▁█
-  ▇██▇▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▇▆▁▄▁▁▁▅▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▅▁▁▄▁▆ ▄
-  75.1 ms       Histogram: log(frequency) by time       175 ms <
-
- Memory estimate: 56.61 MiB, allocs estimate: 481212.
-# MetidaCu
-BenchmarkTools.Trial: 2 samples with 1 evaluation.
- Range (min … max):  9.221 s …   9.255 s  ┊ GC (min … max): 0.00% … 0.29%
- Time  (median):     9.238 s              ┊ GC (median):    0.14%
- Time  (mean ± σ):   9.238 s ± 24.138 ms  ┊ GC (mean ± σ):  0.14% ± 0.20%
-
-  █                                                       █
-  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  9.22 s         Histogram: frequency by time        9.26 s <
-
- Memory estimate: 143.47 MiB, allocs estimate: 2524496.
+BenchmarkTools.Trial: 440 samples with 1 evaluation.
+ Range (min … max):  13.693 ms … 644.360 ms  ┊ GC (min … max):  0.00% … 97.15%
+ Time  (median):     17.967 ms               ┊ GC (median):     0.00%
+ Time  (mean ± σ):   34.118 ms ±  86.080 ms  ┊ GC (mean ± σ):  47.99% ± 18.02%
 
+  █  
+  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▄▄▁▅▄▄▄▄▁▁▁▁▁▁▄ ▆
+  13.7 ms       Histogram: log(frequency) by time       510 ms <
 
+ Memory estimate: 42.28 MiB, allocs estimate: 23491.
 MODEL 4
-# Metida
-BenchmarkTools.Trial: 3 samples with 1 evaluation.
- Range (min … max):  6.738 s …   6.852 s  ┊ GC (min … max): 1.36% … 0.52%
- Time  (median):     6.745 s              ┊ GC (median):    1.09%
- Time  (mean ± σ):   6.779 s ± 63.739 ms  ┊ GC (mean ± σ):  0.99% ± 0.43%
-
-  █  █                                                    █
-  █▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  6.74 s         Histogram: frequency by time        6.85 s <
-
- Memory estimate: 2.33 GiB, allocs estimate: 41657.
-# MetidaNLopt
-BenchmarkTools.Trial: 11 samples with 1 evaluation.
- Range (min … max):  1.317 s …   1.438 s  ┊ GC (min … max): 3.19% … 8.44%
- Time  (median):     1.365 s              ┊ GC (median):    6.01%
- Time  (mean ± σ):   1.365 s ± 35.318 ms  ┊ GC (mean ± σ):  5.52% ± 2.00%
-
-  ▁▁      ▁   ▁         █▁       ▁▁ ▁                     ▁
-  ██▁▁▁▁▁▁█▁▁▁█▁▁▁▁▁▁▁▁▁██▁▁▁▁▁▁▁██▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  1.32 s         Histogram: frequency by time        1.44 s <
-
- Memory estimate: 2.00 GiB, allocs estimate: 138279.
-# MetidaCu
-BenchmarkTools.Trial: 3 samples with 1 evaluation.
- Range (min … max):  7.391 s …   7.432 s  ┊ GC (min … max): 1.10% … 1.43%
- Time  (median):     7.426 s              ┊ GC (median):    1.30%
- Time  (mean ± σ):   7.416 s ± 21.787 ms  ┊ GC (mean ± σ):  1.28% ± 0.17%
-
-  █                                               █       █
-  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█ ▁
-  7.39 s         Histogram: frequency by time        7.43 s <
-
- Memory estimate: 1.87 GiB, allocs estimate: 919345.
+BenchmarkTools.Trial: 4 samples with 1 evaluation.
+ Range (min … max):  4.635 s …    4.880 s  ┊ GC (min … max): 1.77% … 4.69%
+ Time  (median):     4.791 s               ┊ GC (median):    3.57%
+ Time  (mean ± σ):   4.774 s ± 102.530 ms  ┊ GC (mean ± σ):  3.42% ± 1.34%
 
+  █                                █      █                █
+  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
+  4.64 s         Histogram: frequency by time         4.88 s <
 
-MODEL 5
-# Metida
-BenchmarkTools.Trial: 1 sample with 1 evaluation.
- Single result which took 216.945 s (0.04% GC) to evaluate,
- with a memory estimate of 3.91 GiB, over 7229 allocations.
-# MetidaNLopt
-BenchmarkTools.Trial: 2 samples with 1 evaluation.
- Range (min … max):  8.973 s …    9.139 s  ┊ GC (min … max): 3.06% … 2.77%
- Time  (median):     9.056 s               ┊ GC (median):    2.91%
- Time  (mean ± σ):   9.056 s ± 117.413 ms  ┊ GC (mean ± σ):  2.91% ± 0.20%
-
-  █                                                        █
-  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  8.97 s         Histogram: frequency by time         9.14 s <
-
- Memory estimate: 7.86 GiB, allocs estimate: 57558.
-# MetidaCu
-BenchmarkTools.Trial: 2 samples with 1 evaluation.
- Range (min … max):  12.039 s …  12.084 s  ┊ GC (min … max): 1.85% … 1.78%
- Time  (median):     12.062 s              ┊ GC (median):    1.81%
- Time  (mean ± σ):   12.062 s ± 31.774 ms  ┊ GC (mean ± σ):  1.81% ± 0.05%
-
-  █                                                        █
-  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
-  12 s           Histogram: frequency by time         12.1 s <
-
- Memory estimate: 8.31 GiB, allocs estimate: 365031.
-
-=#
-
-#=
-lmm = Metida.LMM(@formula(r2 ~ f), spatdf;
-repeated = Metida.VarEffect(Metida.@covstr(x+y|1), Metida.SPEXP),
-)
-@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
-
-
-spatdf.ci = map(x -> CartesianIndex(x[:x], x[:y]), eachrow(spatdf))
-function Metida.edistance(mx::AbstractMatrix{<:CartesianIndex}, i::Int, j::Int)
-    return sqrt((mx[i, 1][1] - mx[j, 1][1])^2 + (mx[i, 1][2] - mx[j, 1][2])^2)
-end
-lmm = Metida.LMM(@formula(r2 ~ f), spatdf;
-repeated = Metida.VarEffect(Metida.@covstr(ci|1), Metida.SPEXP; coding = Dict(:ci => Metida.RawCoding())),
-)
-@benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
+ Memory estimate: 2.56 GiB, allocs estimate: 40772.
 =#
\ No newline at end of file
diff --git a/test/profile.pb.gz b/test/profile.pb.gz
deleted file mode 100644
index fd252a9..0000000
Binary files a/test/profile.pb.gz and /dev/null differ
diff --git a/test/test.jl b/test/test.jl
index afcf3be..4609279 100644
--- a/test/test.jl
+++ b/test/test.jl
@@ -32,15 +32,15 @@ include("testdata.jl")
     Metida.fit!(lmm)
     @test Metida.m2logreml(lmm) ≈ 25.00077786912235 atol=1E-6
 
-    #Missing
+    # Missing
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0m;
     random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
     )
     Metida.fit!(lmm)
     @test Metida.m2logreml(lmm) ≈ 16.636012616466203 atol=1E-6
 
-    #milmm = Metida.MILMM(lmm, df0m)
-    #Basic, Subject block
+    # milmm = Metida.MILMM(lmm, df0m)
+    # Basic, Subject block
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
     )
@@ -109,6 +109,8 @@ include("testdata.jl")
     @test Metida.confint(lmm, 6)[1] ≈ -0.7630380758015894 atol=1E-4
     @test Metida.confint(lmm; ddf = :residual)[end][1] ≈ -0.6740837049617738 atol=1E-4
     @test Metida.responsename(lmm) == "var"
+    @test Metida.nblocks(lmm) == 5
+    @test Metida.msgnum(lmm.log) == 3
 
     Metida.confint(lmm; ddf = :contain)[end][1] #NOT VALIDATED
     @test size(crossmodelmatrix(lmm), 1) == 6
@@ -130,21 +132,34 @@ include("testdata.jl")
     # AI like algo
     Metida.fit!(lmm; aifirst = true, init = Metida.theta(lmm))
     @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
+
     # Score
     Metida.fit!(lmm; aifirst = :score)
     @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
+
     # AI
     Metida.fit!(lmm; aifirst = :ai)
     @test Metida.m2logreml(lmm) ≈ 16.241112644506067 atol=1E-6
-    #Set user coding
+
+    # Set user coding
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     random = Metida.VarEffect(Metida.@covstr(1 + formulation|subject), Metida.CSH; coding = Dict(:formulation => StatsModels.DummyCoding())),
     )
+
     # Test varlink/rholinkf
     Metida.fit!(lmm; rholinkf = :sqsigm)
     @test Metida.dof_satter(lmm, [0, 0, 0, 0, 0, 1]) ≈ 6.043195705464293 atol=1E-2
     @test Metida.m2logreml(lmm) ≈ 10.314822559210157 atol=1E-6
-    @test_nowarn Metida.fit!(lmm; varlinkf = :sq)
+
+    Metida.fit!(lmm; rholinkf = :atan)
+    @test Metida.m2logreml(lmm) ≈ 10.314837309793571 atol=1E-6
+
+    Metida.fit!(lmm; rholinkf = :psigm)
+    @test Metida.m2logreml(lmm) ≈ 10.86212458333098 atol=1E-6
+
+    Metida.fit!(lmm; varlinkf = :sq)
+    @test Metida.m2logreml(lmm) ≈ 10.314822479530243 atol=1E-6
+
     # Repeated effect only
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     repeated = Metida.VarEffect(Metida.@covstr(formulation|nosubj)),
@@ -157,14 +172,19 @@ include("testdata.jl")
     random = formulation|subject:Metida.DIAG), df0);
     @test  Metida.responsename(lmm) == "log(var)"
 
-    #BE like
+    # BE like
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
     repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
     )
     Metida.fit!(lmm; aifirst = :score)
     @test Metida.m2logreml(lmm) ≈ 10.065238626765524 atol=1E-6
-    #incomplete
+
+    # One thread
+    Metida.fit!(lmm; maxthreads = 1)
+    @test Metida.m2logreml(lmm) ≈ 10.065238626765524 atol=1E-6
+
+    # incomplete
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df1;
     random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.CSH),
     repeated = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
@@ -178,6 +198,7 @@ include("testdata.jl")
     )
     Metida.fit!(lmm)
     @test Metida.m2logreml(lmm, [0.222283, 0.444566]) ≈ Metida.m2logreml(lmm) atol=1E-6
+
     # EXPERIMENTAL
     @test Metida.dof_contain(lmm, 1) == 12
     @test Metida.dof_contain(lmm, 5) == 8
@@ -189,13 +210,12 @@ include("testdata.jl")
 
     # Int dependent variable, function Term in random part
     df0.varint = Int.(ceil.(df0.var2))
-    lmmint = Metida.fit(Metida.LMM, Metida.@lmmformula(varint~formulation,
+    lmmint =  @test_warn "Response variable not <: AbstractFloat" Metida.fit(Metida.LMM, Metida.@lmmformula(varint~formulation,
     random = 1+var^2|subject:Metida.SI), df0)
     Metida.fit!(lmmint)
     @test Metida.m2logreml(lmmint) ≈ 84.23373276096902 atol=1E-6
 
     # Wts
-
     df0.wtsc = fill(0.5, size(df0, 1))
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
@@ -209,8 +229,24 @@ include("testdata.jl")
     fit!(lmm)
     @test Metida.m2logreml(lmm) ≈ 17.823729 atol=1E-6 # TEST WITH SPSS 28
 
+    @test_warn "wts count not equal observations count! wts not used." lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
+    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
+    wts = ones(10))
+
+    # Matrix wts
+    matwts = Symmetric(rand(size(df0,1), size(df0,1)))
+    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
+    random = Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG),
+    wts = matwts)
+    @test_nowarn fit!(lmm)
+
+    # experimental weighted covariance 
+    lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
+    repeated = Metida.VarEffect(Metida.@covstr(1|subject), Metida.SWC(matwts)))
+    @test_nowarn fit!(lmm)
+    @test_nowarn show(io, lmm)
+
     # Repeated vector
-    
     lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
     repeated = [Metida.VarEffect(Metida.@covstr(formulation|subject), Metida.DIAG), Metida.VarEffect(Metida.@covstr(1|subject), Metida.SI)])
     fit!(lmm)
@@ -654,7 +690,7 @@ end
     reml = Metida.m2logreml(lmm)
     @test reml_c ≈ reml
 
-    function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure)
+    function Metida.rmat!(mx, θ, rz, ::CustomCovarianceStructure, ::Int)
         vec = Metida.tmul_unsafe(rz, θ)
         rn    = size(mx, 1)
         if rn > 1
@@ -783,6 +819,24 @@ end
     random = Metida.VarEffect(Metida.@covstr(factor|subject), Metida.DIAG),
     repeated = Metida.VarEffect(Metida.@covstr(1|subject+factor), Metida.ARMA),
     )
+
+    @test_throws ErrorException  Metida.LMM(@formula(var~sequence+period+formulation), df0;)
+    
+
+    @test_throws ErrorException  begin
+        # make cov type
+        struct NewCCS <: Metida.AbstractCovarianceType end
+        function Metida.covstrparam(ct::NewCCS, t::Int)::Tuple{Int, Int}
+            return (t, 1)
+        end
+        # try to apply to repeated effect
+        lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
+        repeated = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CovarianceType(NewCCS())),
+        )
+        # try to get V 
+        Metida.vmatrix([1.0, 1.0, 1.0], lmm, 1) 
+    end
+
     # Error messages
     io = IOBuffer();
     lmm = Metida.LMM(@formula(response ~ 1 + factor*time), ftdf2;