diff --git a/Project.toml b/Project.toml index 4b044b1..86d8102 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ authors = ["Vladimir Arnautov "] 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" diff --git a/src/fit.jl b/src/fit.jl index ec3b85b..a6d0788 100644 --- a/src/fit.jl +++ b/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 diff --git a/src/gmat.jl b/src/gmat.jl index c6346c9..452cf9f 100644 --- a/src/gmat.jl +++ b/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 diff --git a/src/reml.jl b/src/reml.jl index 5b8500c..10c29e5 100644 --- a/src/reml.jl +++ b/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π) diff --git a/src/rmat.jl b/src/rmat.jl index 04fcaf7..7681c7e 100644 --- a/src/rmat.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 2bab8fa..3a8efad 100644 --- a/src/utils.jl +++ b/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 ################################################################################ diff --git a/test/profile.pb.gz b/test/profile.pb.gz new file mode 100644 index 0000000..fd252a9 Binary files /dev/null and b/test/profile.pb.gz differ