micro opt

Project.toml

@@ -6,7 +6,6 @@ authors = ["Vladimir Arnautov <[email protected]>"]
 version = "0.15.1"
 
 [deps]
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -20,7 +19,6 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
-Compat = "3.46, 4"
 DiffResults = "1"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 ForwardDiff = "0.10"

src/fit.jl

@@ -63,7 +63,7 @@ Fit LMM model.
 * `refitinit` - true/false - if `true` - use last values for initial condition  (`false` by default)
 * `optmethod` - Optimization method. Look at Optim.jl documentation. (Newton by default)
 * `singtol` - singular tolerance = 1e-8
-* `maxthreads` - maximum threads = num_cores()
+* `maxthreads` - maximum threads = min(num_cores(), Threads.nthreads())
 
 """
 function fit!(lmm::LMM{T}; kwargs...) where T
@@ -87,7 +87,7 @@ function fit!(lmm::LMM{T}; kwargs...) where T
     :refitinit ∈ kwkeys ? refitinit = kwargs[:refitinit] : refitinit = false
     :optmethod ∈ kwkeys ? optmethod = kwargs[:optmethod] : optmethod = :default
     :singtol ∈ kwkeys ? singtol = kwargs[:singtol] : singtol = 1e-8
-    :maxthreads ∈ kwkeys ? maxthreads = kwargs[:maxthreads] : maxthreads = num_cores()
+    :maxthreads ∈ kwkeys ? maxthreads = kwargs[:maxthreads] : maxthreads = min(num_cores(), Threads.nthreads())
 
     # If model was fitted, previous results can be used if `refitinit` == true
     # Before fitting clear log

src/gmat.jl

@@ -148,7 +148,7 @@ function gmat!(mx, θ, ::TOEP_)
     if s > 1
         for n = 2:s
             @inbounds @simd for m = 1:n-1
-                mx[m, n] = de * θ[n-m+1]
+                mx[m, n] = de * θ[n - m + 1]
             end
         end
     end
@@ -177,8 +177,9 @@ function gmat!(mx, θ, ::TOEPH_)
     end
     if s > 1
         for n = 2:s
+            @inbounds mxnn = mx[n, n]
             @inbounds @simd for m = 1:n-1
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[n-m+s]
+                mx[m, n] = mx[m, m] * mxnn * θ[n - m + s]
             end
         end
     end
@@ -195,8 +196,9 @@ function gmat!(mx, θ, ct::TOEPHP_)
     end
     if s > 1 && ct.p > 1
         for m = 1:s - 1
+            @inbounds mxmm = mx[m, m]
             for n = m + 1:(m + ct.p - 1 > s ? s : m + ct.p - 1)
-                @inbounds mx[m, n] = mx[m, m] * mx[n, n] * θ[n - m + s]
+                @inbounds mx[m, n] = mxmm * mx[n, n] * θ[n - m + s]
             end
         end
     end
@@ -213,8 +215,9 @@ function gmat!(mx, θ, ::UN_)
     end
     if s > 1
         for n = 2:s
+            @inbounds mxnn = mx[n, n]
             @inbounds @simd for m = 1:n - 1
-                mx[m, n] = mx[m, m] * mx[n, n] * θ[s + tpnum(m, n, s)]
+                mx[m, n] = mx[m, m] * mxnn * θ[s + tpnum(m, n, s)]
             end
         end
     end

src/reml.jl

@@ -33,7 +33,7 @@ end
 ################################################################################
 #                     REML without provided β
 ################################################################################
-function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 16) where T # Main optimization way - make gradient / hessian analytical / semi-analytical functions
+function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 4) where T # Main optimization way - make gradient / hessian analytical / semi-analytical functions
     n             = length(lmm.covstr.vcovblock)
     N             = length(lmm.data.yv)
     c             = (N - lmm.rankx)*log(2π)

src/rmat.jl

@@ -309,8 +309,9 @@ function unrmat(θ::AbstractVector{T}, rz) where T
     end
     if rm > 1
         for m = 1:rm - 1
+            @inbounds mxmm = mx[m, m] 
             @inbounds @simd for n = m + 1:rm
-                mx[m, n] += mx[m, m] * mx[n, n] * θ[rm+tpnum(m, n, rm)]
+                mx[m, n] += mxmm * mx[n, n] * θ[rm + tpnum(m, n, rm)]
             end
         end
     end

src/utils.jl

@@ -204,8 +204,8 @@ end
 function logreml(lmm)
     -m2logreml(lmm)/2
 end
-function m2logreml(lmm, theta)
-    reml_sweep_β(lmm, LMMDataViews(lmm), theta)[1]
+function m2logreml(lmm, theta; maxthreads::Int = num_cores())
+    reml_sweep_β(lmm, LMMDataViews(lmm), theta; maxthreads = maxthreads)[1]
 end
 ################################################################################