Merge pull request #32 from PharmCat/dev14

optim, minor fix, update

.github/workflows/Documenter.yml

@@ -23,7 +23,7 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 30
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: "1"

.github/workflows/Tier1.yml

@@ -40,7 +40,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -49,6 +49,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v2
+      - uses: codecov/codecov-action@v3
         with:
           files: lcov.info

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.14.4"
+version = "0.14.5"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

docs/src/api.md

@@ -141,14 +141,24 @@ Metida.hessian
 Metida.lcontrast
 ```
 
+### Metida.nblocks
+```@docs
+Metida.nblocks
+```
+
 ### Metida.rand
 ```@docs
 Metida.rand
 ```
 
 ### Metida.rand!
 ```@docs
-Metida.rand
+Metida.rand!
+```
+
+### Metida.raneff
+```@docs
+Metida.raneff
 ```
 
 ### Metida.rankx
@@ -215,7 +225,7 @@ Metida.confint
 
 ### Metida.crossmodelmatrix
 ```@docs
-StatsBase.crossmodelmatrix
+Metida.crossmodelmatrix
 ```
 
 ### Metida.dof
@@ -253,21 +263,11 @@ Metida.loglikelihood
 Metida.modelmatrix
 ```
 
-### Metida.nblocks
-```@docs
-Metida.nblocks
-```
-
 ### Metida.nobs
 ```@docs
 Metida.nobs
 ```
 
-### Metida.raneff
-```@docs
-Metida.raneff
-```
-
 ### Metida.response
 ```@docs
 Metida.response

docs/src/details.md

@@ -33,7 +33,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}| \\

docs/src/nlopt.md

@@ -19,3 +19,9 @@ random   = VarEffect(@covstr(formulation|subject), CSH),
 repeated = VarEffect(@covstr(formulation|subject), VC))
 fit!(lmm; solver = :nlopt)
 ```
+
+
+NLopt is a free/open-source library for nonlinear optimization, providing a common interface for a number of different free optimization routines available online as well as original implementations of various other algorithms.
+
+
+Optimization with NLopt.jl using gradient-free algirithms is less stable, that why two-step optimization schema used. Results can be slightly different for differens OS and Julia versions. Always look into logs. 

src/Metida.jl

@@ -8,7 +8,7 @@ using  ProgressMeter, LinearAlgebra, ForwardDiff, DiffResults, Random, Optim, Li
 
 import MetidaBase: StatsBase, StatsModels, CategoricalArrays, Distributions
 
-import MetidaBase.CategoricalArrays: CategoricalArray
+import MetidaBase.CategoricalArrays: CategoricalArray, AbstractCategoricalVector
 import MetidaBase.Distributions: Normal, TDist, FDist, Chisq, MvNormal, FullNormal, ccdf, cdf, quantile
 
 import MetidaBase: Tables, MetidaModel, AbstractCovarianceStructure, AbstractCovmatMethod, AbstractCovarianceType, AbstractLMMDataBlocks, MetidaTable, metida_table, PrettyTables, indsdict!

src/gmat.jl

@@ -57,9 +57,10 @@ function gmat!(mx, θ, ::AR_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = de * θ[2] ^ (n - m)
+        @inbounds θ2 = θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = de * θ2 ^ (n - m)
             end
         end
     end
@@ -72,14 +73,16 @@ function gmat!(mx, θ, ::ARH_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[end] ^ (n - m)
+        θe = last(θ)
+        for n = 2:s
+            mxnn = mx[n, n]
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = mx[m, m] * mxnn * θe ^ (n - m)
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m] 
     end
     mx
 end
@@ -88,9 +91,10 @@ function gmat!(mx, θ, ::CS_)
     s = size(mx, 1)
     mx .= θ[1]^2
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = mx[m, m] * θ[2]
+        mxθ2 = θ[1]^2 * θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = mxθ2
             end
         end
     end
@@ -103,14 +107,15 @@ function gmat!(mx, θ, ::CSH_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[end]
+        for n = 2:s
+            @inbounds mxnθe = mx[n, n] * last(θ)
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = mx[m, m] * mxnθe
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -123,9 +128,10 @@ function gmat!(mx, θ, ::ARMA_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
-                mx[m, n] = de * θ[2] * θ[3] ^ (n - m - 1)
+        deθ2 = de * θ[2]
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
+                mx[m, n] = deθ2 * θ[3] ^ (n - m - 1)
             end
         end
     end
@@ -139,8 +145,8 @@ function gmat!(mx, θ, ::TOEP_)
         mx[i, i] = de
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = de * θ[n-m+1]
             end
         end
@@ -169,14 +175,14 @@ function gmat!(mx, θ, ::TOEPH_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s]
             end
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -194,7 +200,7 @@ function gmat!(mx, θ, ct::TOEPHP_)
         end
     end
     @inbounds @simd for m = 1:s
-        mx[m, m] = mx[m, m] * mx[m, m]
+        mx[m, m] *= mx[m, m]
     end
     mx
 end
@@ -205,15 +211,14 @@ function gmat!(mx, θ, ::UN_)
         mx[m, m] = θ[m]
     end
     if s > 1
-        for m = 1:s - 1
-            @inbounds @simd for n = m + 1:s
+        for n = 2:s
+            @inbounds @simd for m = 1:n-1
                 mx[m, n] = mx[m, m] * mx[n, n] * θ[s+tpnum(m, n, s)]
             end
         end
     end
     @inbounds @simd for m = 1:s
-        v = mx[m, m]
-        mx[m, m] = v * v
+        mx[m, m] *= mx[m, m]
     end
     mx
 end

src/linearalgebra.jl

@@ -2,17 +2,19 @@
 
 # Fine
 """
+    mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
+
 θ + A * B * A'
 
 Change θ (only upper triangle). B is symmetric.
 """
 @noinline function mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
     axb  = axes(B, 1)
     sa   = size(A, 1)
-    @simd  for j ∈ axb
-        @simd for i ∈ axb
+    for j ∈ axb
+        for i ∈ axb
             @inbounds Bij = B[i, j]
-            @simd  for n ∈ 1:sa
+           for n ∈ 1:sa
                 @inbounds Anj = A[n, j]
                 BijAnj = Bij * Anj
                 @simd for m ∈ 1:n
@@ -31,6 +33,8 @@ function mulαβαtinc!(θ::AbstractMatrix{T}, A::AbstractMatrix{T}, B::Abstract
 end
 =#
 """
+    mulαβαtinc!(θ::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix, alpha)
+
 θ + A * B * A' * alpha
 
 Change θ (only upper triangle). B is symmetric.
@@ -39,10 +43,10 @@ Change θ (only upper triangle). B is symmetric.
     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)
-    @simd for j ∈ axb
-        @simd for i ∈ axb
+   for j ∈ axb
+        for i ∈ axb
             @inbounds Bij = B[i, j]
-            @simd  for n ∈ 1:sa 
+            for n ∈ 1:sa 
                 @inbounds Anj = A[n, j]
                 BijAnjalpha = Bij * Anj * alpha
                 @simd for m ∈ 1:n
@@ -54,6 +58,8 @@ Change θ (only upper triangle). B is symmetric.
     θ
 end
 """
+    mulαβαtinc!(θ::AbstractVector{T}, A::AbstractMatrix, B::AbstractMatrix, a::AbstractVector, b::AbstractVector, alpha) where T
+
 θ + A * B * (a - b) * alpha
 
 Change θ (only upper triangle). B is symmetric.
@@ -62,9 +68,9 @@ Change θ (only upper triangle). B is symmetric.
     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)
-    @simd for i ∈ axb
+    for i ∈ axb
         @inbounds abi = a[i] - b[i]
-        @simd for j ∈ axb
+       for j ∈ axb
             @inbounds Bji = B[j, i]
             Bjiabialpha = Bji * abi * alpha
             @simd for m ∈ 1:sa
@@ -95,13 +101,13 @@ use only upper triangle of V
         return -V[1, 1] * (y[1] - cs)^2
     end
     c = zeros(T, q)
-    @simd for m = 1:p
+    for m = 1:p
         @inbounds βm = β[m]
         @simd for n = 1:q
             @inbounds c[n] += X[n, m] * βm
         end
     end
-    @simd for m = 2:q
+    for m = 2:q
         @inbounds ycm2 = (y[m] - c[m])*2
         @simd for n = 1:m-1
             @inbounds θ -= V[n, m] * (y[n] - c[n]) * ycm2
@@ -114,6 +120,8 @@ use only upper triangle of V
 end
 
 """
+    mulαtβinc!(θ::AbstractVector{T}, A::AbstractMatrix, b::AbstractVector) where T
+
 θ + A' * b
 
 Change θ.
@@ -124,7 +132,7 @@ Change θ.
     p = size(A, 2)
     for n in 1:p
         θn = zero(T)
-        for m in 1:q
+        @simd for m in 1:q
             @inbounds θn += b[m] * A[m, n]
         end
         @inbounds θ[n] += θn
@@ -136,7 +144,7 @@ end
     vec = zeros(T, size(rz, 1))
     for i ∈ axes(rz, 2)
         @inbounds θi = θ[i]
-        for r ∈ axes(rz, 1)
+        @simd for r ∈ axes(rz, 1)
             @inbounds vec[r] += rz[r, i] * θi
         end
     end
@@ -146,8 +154,8 @@ end
 @inline function diag!(f, v, m)
     l = checksquare(m)
     l == length(v) || error("Length not equal")
-    for i = 1:l
-        v[i] = f(m[i, i])
+    @simd for i = 1:l
+        @inbounds v[i] = f(m[i, i])
     end
     v
 end

src/lmm.jl

@@ -238,7 +238,7 @@ function Base.show(io::IO, lmm::LMM)
         end
         pretty_table(io, mx; show_header = false, alignment=:l, tf = tf_borderless)
     else
-        if !any(x-> x.type == :WARN, lmm.log)
+        if !any(x-> x.type == :ERROR, lmm.log)
             println(io, "Not fitted.")
         else
             printstyled(io, "Not fitted. See error(s) in log.\n"; color = :red)