rankUpdate! method for Diagonal + Dense (#446)

* add oxide dataset from nlme

* add oxide to model cache

* new rankupdate method + test

* version bump

Author: palday

Artifacts.toml

@@ -5,21 +5,21 @@
 [TestData]
 # compute this using
 # using Tar, Inflate, SHA
-# filename = "download?version=4" # I just used wget for the URL below and this is how it saved it
+# filename = "download?version=5" # I just used wget for the URL below and this is how it saved it
 # println("sha256: ", bytes2hex(open(sha256, filename)))
 # println("git-tree-sha1: ", Tar.tree_hash(IOBuffer(inflate_gzip(filename))))
 # from https://julialang.github.io/Pkg.jl/dev/artifacts/
-git-tree-sha1 = "75260d131b693f26e5834adf855f4c35d627348d"
+git-tree-sha1 = "91132469677f725c2e4097493ae8b1d566f90a3f"
 lazy = false
 
     [[TestData.download]]
-    # this is the SHA from https://osf.io/djaqb/download?version=4
-    sha256 = "d179cadfb27fc638fd8a06b0d58f5f2916f786c5719ab22de93fedb6129726f4"
+    # this is the SHA from https://osf.io/djaqb/download?version=5
+    sha256 = "8040933246179d4b46cf37e4dd1076752102dfcb4bd937d5f12c0f724775e119"
     # when updating this, make sure to change to change the version number,
     # because if the version number isn't included, it will always point to the
     # latest version, which means it will break existing users when we update
     # between releases.
-    url = "https://osf.io/djaqb/download?version=4"
+    url = "https://osf.io/djaqb/download?version=5"
 
     # for future work on using xz-compressed data:
     # Julia invokes wget without using HTTP metadata, so we need the link

Project.toml

@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <[email protected]>", "Douglas Bates <[email protected]>", "Jose Bayoan Santiago Calderon <[email protected]>"]
-version = "3.1.1"
+version = "3.1.2"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"

src/linalg/rankUpdate.jl

@@ -49,6 +49,18 @@ function rankUpdate!(C::HermOrSym, A::BlockedSparse, α, β)
     rankUpdate!(C, sparse(A), α, β)
 end
 
+function rankUpdate!(C::HermOrSym{T,Diagonal{T,Vector{T}}}, A::StridedMatrix{T}, α, β) where {T,S}
+    Cdiag = C.data.diag
+    @. Cdiag = β * Cdiag
+
+    for i in 1:length(Cdiag)
+        Arow = view(A, i, :)
+        Cdiag[i] = Cdiag[i] + α * Arow'Arow
+    end
+
+    C
+end
+
 function rankUpdate!(
     C::HermOrSym{T,Diagonal{T,Vector{T}}},
     A::SparseMatrixCSC{T},

test/modelcache.jl

@@ -9,6 +9,7 @@ using MixedModels: dataset
 )
 
 @isdefined(fms) || const global fms = Dict(
+    :oxide => [@formula(Thickness ~ 1 + (1|Lot/Wafer))],
     :dyestuff => [@formula(yield ~ 1 + (1|batch))],
     :dyestuff2 => [@formula(yield ~ 1 + (1|batch))],
     :d3 => [@formula(y ~ 1 + u + (1+u|g) + (1+u|h) + (1+u|i))],

test/pls.jl

@@ -407,6 +407,12 @@ end
     @test "-0.89" in tokens
 end
 
+@testset "oxide" begin
+    # this model has an interesting structure with two diagonal blocks
+    m = only(models(:oxide))
+    @test all(isapprox.(m.θ, [1.689182746, 2.98504262]; atol=1e-4))
+end
+
 @testset "Rank deficient" begin
     rng = MersenneTwister(0);
     x = rand(rng, 100);