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 b95c833..e8950eb 100644 --- a/src/Metida.jl +++ b/src/Metida.jl @@ -20,6 +20,7 @@ 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/reml.jl b/src/reml.jl index 10c29e5..bc30de5 100644 --- a/src/reml.jl +++ b/src/reml.jl @@ -88,7 +88,12 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 4) where T # end θ₁ = sum(accθ₁) θ₂ = sum(accθ₂) - βm = sum(accβm) + 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(θ₂) diff --git a/src/utils.jl b/src/utils.jl index f2de44c..a142c87 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -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 diff --git a/test/devtest.jl b/test/devtest.jl index d70c114..e1e52d8 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,86 @@ 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 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 - 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 - 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 -################################################################################ - 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() - -#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. - - -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. - - -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. - - -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. - - -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. - -=# +for i = 1:4 + display(b[i]) + push!(results, (now(), "Model $i", minimum(b[i]).time, minimum(b[i]).memory, minimum(b[i]).allocs)) +end #= -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 +BenchmarkTools.Trial: 411 samples with 1 evaluation. + Range (min … max): 11.681 ms … 1.511 s ┊ GC (min … max): 0.00% … 99.24% + Time (median): 14.456 ms ┊ GC (median): 0.00% + Time (mean ± σ): 36.240 ms ± 170.372 ms ┊ GC (mean ± σ): 57.08% ± 11.90% + + █ + █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▄ ▆ + 11.7 ms Histogram: log(frequency) by time 1.42 s < + + Memory estimate: 19.38 MiB, allocs estimate: 31265. +BenchmarkTools.Trial: 23 samples with 1 evaluation. + Range (min … max): 657.256 ms … 719.396 ms ┊ GC (min … max): 0.00% … 4.31% + Time (median): 670.301 ms ┊ GC (median): 0.00% + Time (mean ± σ): 677.402 ms ± 16.735 ms ┊ GC (mean ± σ): 1.46% ± 2.20% + + █ ▁ + ▆▆▆▁▁▆█▆█▁▁▁▆▁▁▁▁▁▁▁▁▁▁▆▆▆▁▁▁▆▆▆▁▁▁▁▆▁▆▆▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁ + 657 ms Histogram: frequency by time 719 ms < + + Memory estimate: 132.83 MiB, allocs estimate: 5086. +BenchmarkTools.Trial: 297 samples with 1 evaluation. + Range (min … max): 17.368 ms … 1.552 s ┊ GC (min … max): 0.00% … 98.55% + Time (median): 20.409 ms ┊ GC (median): 0.00% + Time (mean ± σ): 54.268 ms ± 197.857 ms ┊ GC (mean ± σ): 60.26% ± 15.96% + + █ + █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▅▅ ▅ + 17.4 ms Histogram: log(frequency) by time 1.19 s < + + Memory estimate: 42.28 MiB, allocs estimate: 23492. +BenchmarkTools.Trial: 4 samples with 1 evaluation. + Range (min … max): 4.696 s … 4.993 s ┊ GC (min … max): 2.54% … 8.32% + Time (median): 4.727 s ┊ GC (median): 2.74% + Time (mean ± σ): 4.786 s ± 138.981 ms ┊ GC (mean ± σ): 4.11% ± 2.86% + + █ ██ █ + █▁▁▁▁██▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁ + 4.7 s Histogram: frequency by time 4.99 s < + + Memory estimate: 2.56 GiB, allocs estimate: 40772. =# \ No newline at end of file diff --git a/test/test.jl b/test/test.jl index fff428a..bd6ecc0 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), ) @@ -130,21 +130,26 @@ 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) + # Repeated effect only lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0; repeated = Metida.VarEffect(Metida.@covstr(formulation|nosubj)), @@ -157,14 +162,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 +188,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 +200,12 @@ include("testdata.jl") # Int dependent variable, function Term in random part df0.varint = Int.(ceil.(df0.var2)) - @test_warn "Response variable not <: AbstractFloat" 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,7 +219,6 @@ 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)) @@ -227,7 +236,6 @@ include("testdata.jl") @test_nowarn fit!(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)