From f6801f069626e8bdfa111e16732f0e5cce9b2026 Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Fri, 27 Oct 2017 12:17:43 -0500
Subject: [PATCH 1/4] Add Ac_mul_B method needed in example from
 marketingattribution.com

---
 src/modelterms.jl | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/src/modelterms.jl b/src/modelterms.jl
index a867f7b7c..3f936ff7a 100644
--- a/src/modelterms.jl
+++ b/src/modelterms.jl
@@ -464,6 +464,24 @@ function Base.Ac_mul_B(A::ScalarFactorReTerm{T,V,R},
     end
 end
 
+function Base.Ac_mul_B(A::VectorFactorReTerm, B::ScalarFactorReTerm)
+    nzeros = copy(A.wtz)
+    k, n = size(nzeros)
+    rowind = Matrix{Int32}(k, n)
+    refs = A.f.refs
+    bwtz = B.wtz
+    for j in 1:n
+        bwtzj = bwtz[j]
+        offset = refs[j] - one(eltype(refs))
+        for i in 1:k
+            rowind[i, j] = i + offset
+            nzeros[i, j] += bwtzj
+        end
+    end
+    sparse(vec(rowind), Vector{Int32}(repeat(B.f.refs, inner=k)), vec(nzeros),
+           k * nlevs(A), nlevs(B))
+end
+
 function Base.Ac_mul_B(A::VectorFactorReTerm{T}, B::VectorFactorReTerm{T}) where T
     if A === B
         l = vsize(A)

From 38e5dee3302deac0678c0da97524c072885ac351 Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Fri, 27 Oct 2017 13:07:27 -0500
Subject: [PATCH 2/4] Change vectorized call to broadcast

---
 src/mixedmodel.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mixedmodel.jl b/src/mixedmodel.jl
index fde83ec00..dc686c41b 100644
--- a/src/mixedmodel.jl
+++ b/src/mixedmodel.jl
@@ -36,7 +36,7 @@ function StatsBase.coeftable(m::MixedModel)
     fe = fixef(m)
     se = stderr(m)
     z = fe ./ se
-    pvalue = ccdf(Chisq(1), abs2.(z))
+    pvalue = ccdf.(Chisq(1), abs2.(z))
     CoefTable(hcat(fe, se, z, pvalue), ["Estimate", "Std.Error", "z value", "P(>|z|)"],
         fixefnames(m), 4)
 end

From 531ca2db4854542e5f86b95841e839e8badb202b Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Fri, 27 Oct 2017 13:07:52 -0500
Subject: [PATCH 3/4] Fix recently added method.  Add test.

---
 src/modelterms.jl    | 4 ++--
 test/FactorReTerm.jl | 4 ++++
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/modelterms.jl b/src/modelterms.jl
index 3f936ff7a..d68bdcbfc 100644
--- a/src/modelterms.jl
+++ b/src/modelterms.jl
@@ -472,10 +472,10 @@ function Base.Ac_mul_B(A::VectorFactorReTerm, B::ScalarFactorReTerm)
     bwtz = B.wtz
     for j in 1:n
         bwtzj = bwtz[j]
-        offset = refs[j] - one(eltype(refs))
+        offset = k * (refs[j] - 1)
         for i in 1:k
             rowind[i, j] = i + offset
-            nzeros[i, j] += bwtzj
+            nzeros[i, j] *= bwtzj
         end
     end
     sparse(vec(rowind), Vector{Int32}(repeat(B.f.refs, inner=k)), vec(nzeros),
diff --git a/test/FactorReTerm.jl b/test/FactorReTerm.jl
index f7d5f3618..020728909 100644
--- a/test/FactorReTerm.jl
+++ b/test/FactorReTerm.jl
@@ -116,6 +116,10 @@ end
     @test isa(vrp, UniformBlockDiagonal{Float64})
     @test size(vrp) == (36, 36)
 
+    scl = ScalarFactorReTerm(slp[:G], Array(slp[:U]), Array(slp[:U]), :G, ["U"], 1.0)
+
+    @test sparse(corr)'sparse(scl) == corr'scl
+
     @test MixedModels.Λ_mul_B!(Vector{Float64}(36), corr, ones(36)) == repeat([0.5, 1.0], outer=18)
     
     @testset "reweight!" begin

From 17c4d17e71b81548dcb3320d6c89c15ccd19f95a Mon Sep 17 00:00:00 2001
From: Douglas Bates <dmbates@gmail.com>
Date: Fri, 27 Oct 2017 13:59:14 -0500
Subject: [PATCH 4/4] Methods for special case of very sparse nested factors.

---
 src/linalg/rankUpdate.jl | 8 +++++---
 src/pls.jl               | 4 +++-
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/src/linalg/rankUpdate.jl b/src/linalg/rankUpdate.jl
index a75a8590a..03bc1bbb9 100644
--- a/src/linalg/rankUpdate.jl
+++ b/src/linalg/rankUpdate.jl
@@ -60,9 +60,11 @@ function rankUpdate!(α::T, A::SparseMatrixCSC{T}, C::Diagonal{T}) where T <: Nu
     rv = rowvals(A)
     for j in 1:n
         nzr = nzrange(A, j)
-        length(nzr) == 1 || throw(ArgumentError("A*A' has off-diagonal elements"))
-        k = nzr[1]
-        @inbounds dd[rv[k]] += α * abs2(nz[k])
+        if !isempty(nzr)
+            length(nzr) == 1 || throw(ArgumentError("A*A' has off-diagonal elements"))
+            k = nzr[1]
+            @inbounds dd[rv[k]] += α * abs2(nz[k])
+        end
     end
     C
 end
diff --git a/src/pls.jl b/src/pls.jl
index e8e1f6e30..3cf8a39fc 100644
--- a/src/pls.jl
+++ b/src/pls.jl
@@ -5,10 +5,12 @@ Convert sparse `S` to `Diagonal` if `S` is diagonal or to `full(S)` if
 the proportion of nonzeros exceeds `threshold`.
 """
 function densify(S::SparseMatrixCSC, threshold::Real = 0.3)
+    dropzeros!(S)
     m, n = size(S)
     if m == n && isdiag(S)  # convert diagonal sparse to Diagonal
         Diagonal(diag(S))
-    elseif nnz(S)/(*(size(S)...)) ≤ threshold # very sparse matrices left as is
+    elseif nnz(S)/(*(size(S)...)) ≤ threshold ||   # very sparse matrices left as is
+        all(d -> iszero(d) || d == 1, diff(S.colptr))
         S
     else
         full(S)