Skip to content

Commit

Permalink
test cover, minor opt
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Aug 22, 2024
1 parent 8e7229b commit 6085659
Show file tree
Hide file tree
Showing 5 changed files with 92 additions and 35 deletions.
25 changes: 22 additions & 3 deletions src/reml.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,14 @@ function reml_sweep_β(lmm, data, θ::Vector{T}; maxthreads::Int = 4) where T #
#-----------------------------------------------------------------------
end
θ₁ = sum(accθ₁)
θ₂ = sum(accθ₂)
#θ₂ = 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]
Expand Down Expand Up @@ -194,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

Check warning on line 208 in src/reml.jl

View check run for this annotation

Codecov / codecov/patch

src/reml.jl#L206-L208 

Added lines #L206 - L208 were not covered by tests
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

Check warning on line 215 in src/reml.jl

View check run for this annotation

Codecov / codecov/patch

src/reml.jl#L213-L215 

Added lines #L213 - L215 were not covered by tests
end
βtc = accβm[1]
# Beta calculation
copyto!(θ₂, θ₂tc)
ldθ₂, info = LinearAlgebra.LAPACK.potrf!('U', θ₂tc)
Expand Down
2 changes: 1 addition & 1 deletion src/rmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ Base.@propagate_inbounds function rmat!(mx, θ, ::AbstractMatrix, ::SI_, ::Int)
@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)
Expand Down
4 changes: 2 additions & 2 deletions src/sweep.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ 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)
Expand Down Expand Up @@ -44,7 +44,7 @@ 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))
Expand Down
65 changes: 37 additions & 28 deletions test/devtest.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,71 +27,80 @@ b = Vector{Any}(undef, 4)
# Metida
################################################################################
# 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),
)
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),
)
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),
)
b[3] = @benchmark Metida.fit!($lmm, hes = false; maxthreads = 16) seconds = 15
# MODEL 4
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

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

#=
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%
MODEL 1
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%
█
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▄ ▆
11.7 ms Histogram: log(frequency) by time 1.42 s <
█
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▁▄▃ ▆
8.19 ms Histogram: log(frequency) by time 533 ms <
Memory estimate: 19.38 MiB, allocs estimate: 31265.
MODEL 2
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%
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%
█ ▁
▆▆▆▁▁▆█▆█▁▁▁▆▁▁▁▁▁▁▁▁▁▁▆▆▆▁▁▁▆▆▆▁▁▁▁▆▁▆▆▁▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆ ▁
657 ms Histogram: frequency by time 719 ms <
▁ ▁▁ ▁ ██ ▁ ▁ ▁▁ ▁ ▁ ▁ ▁ █ ▁ █ ▁ ▁
█▁▁▁▁▁██▁█▁██▁▁█▁█▁▁▁██▁▁▁█▁▁█▁▁█▁▁▁█▁▁▁█▁▁▁▁█▁█▁▁▁█▁▁▁▁▁▁▁▁█ ▁
658 ms Histogram: frequency by time 704 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%
Memory estimate: 132.83 MiB, allocs estimate: 5085.
MODEL 3
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%
█
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▅▅ ▅
17.4 ms Histogram: log(frequency) by time 1.19 s <
█
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▄▄▄▁▅▄▄▄▄▁▁▁▁▁▁▄ ▆
13.7 ms Histogram: log(frequency) by time 510 ms <
Memory estimate: 42.28 MiB, allocs estimate: 23492.
Memory estimate: 42.28 MiB, allocs estimate: 23491.
MODEL 4
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%
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.7 s Histogram: frequency by time 4.99 s <
█ █ █ █
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
4.64 s Histogram: frequency by time 4.88 s <
Memory estimate: 2.56 GiB, allocs estimate: 40772.
=#
31 changes: 30 additions & 1 deletion test/test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -148,7 +150,15 @@ include("testdata.jl")
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;
Expand Down Expand Up @@ -234,6 +244,7 @@ include("testdata.jl")
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;
Expand Down Expand Up @@ -808,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;
Expand Down

0 comments on commit 6085659

Please sign in to comment.