An alternative amalgamate strategy

[dmbates]

• indmat extractor for ReMat returns a Bool indicator matrix of potential nonzeros in lambda
• use BlockDiagonals package to create indmat of amalgamated ReMats from original indmat
• BlockDiagonals creates some conflict with BlockArrays but BlockArrays should be removed anyway
• some rewrites according to the github.com/invenia/BlueStyle guide
• remove DataFramesMeta as a test dependency

[codecov-io]

Codecov Report

Merging #216 into master will decrease coverage by 0.16%.
The diff coverage is 94.87%.

@@            Coverage Diff             @@
##           master     #216      +/-   ##
==========================================
- Coverage   93.25%   93.09%   -0.17%     
==========================================
  Files          18       18              
  Lines        1290     1288       -2     
==========================================
- Hits         1203     1199       -4     
- Misses         87       89       +2

Impacted Files	Coverage Δ
src/utilities.jl	84% <ø> (ø)	⬆️
src/MixedModels.jl	100% <ø> (ø)	⬆️
src/varcorr.jl	100% <100%> (ø)	⬆️
src/linearmixedmodel.jl	96.35% <100%> (-0.05%)	⬇️
src/remat.jl	96.04% <93.1%> (-0.78%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ef49e5e...c7a199c. Read the comment docs.

[palday]

Looks good, only a minor style suggestion.

[palday]

How much work will it take to actually remove / write internal replacements for BlockArrays?

[Nosferican]

There aren't that many people familiar with BlockArrays. @dlfivefifty would be the person best positioned to answer whether heterogenous types of blocks could be supported or how to best have an internal workaround (e.g., other than just copying BlockArrays 0.7 as internal module...). The incompatibility is currently affecting a few packages such as Pumas/Bioequivalence which are meant to be working together, but currently can't coexist in the same environment due to the issue.
@dmbates, would it be possible to promote the blocks to the same type?

[palday]

I'm guessing no because some of them are sparse and some of them are not, and the sparse ones may have different substructures of sparsity.

[Nosferican]

Would it be feasible to have a copy of BlockArrays 0.7 version code bundled as an internal module in case there isn't a good solution in the meantime?

[palday]

I haven't looked at how the indexing is actually done in the linear algebra, but assuming that the blocks are used and everything else is ignored (as zeroing out), it might be possible to replace BlockArrays with a dictionary indexed by Block / ordered pair. You lose the ability to index between blocks (at least without adding a wrapper around that / using an indexer that returns a default zero value there) and the sense of overall indexing, but otherwise that should really easy to implement.

[Nosferican]

I couldn't find it, but some package had a simple wrapper view for extending indexing to block-like matrices (a non BlockArrays). At this point I think a temporary solution even if not ideal should be the priority.

[dlfivefifty]

I don't really understand the issue, there shouldn't be any reason to stay on BlockArrays v0.7, and I don't see why BlockDiagonals.jl and BlockArrays.jl should be incompatible.

[Nosferican]

The current implementation used BlockArrays where some blocks are sparse and others dense. Is this supported (or could be) in current versions of BlockArrays?

[dlfivefifty]

Yes:

julia> B = Matrix{AbstractMatrix{Float64}}(undef,2,1)

julia> B[1,1] = fill(1.0,2,2)

julia> B[2,1] = sparse(I,2,2)

julia> A = mortar(B)
2×1-blocked 4×2 BlockArray{Float64,2,Array{AbstractArray{Float64,2},2},BlockArrays.BlockSizes{2,Tuple{Array{Int64,1},Array{Int64,1}}}}:
 1.0  1.0
 1.0  1.0
 ────────
 1.0  0.0
 0.0  1.0

julia> A[Block(1,1)]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  1.0

julia> A[Block(2,1)]
2×2 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 1]  =  1.0
  [2, 2]  =  1.0

[palday]

I think the "incompatibility" @dmbates was referring to was the collision in imported method names and thus the necessity of using fully qualified names for nblocks and blocksize.

[dlfivefifty]

The "right" way to resolve this is to move out the core features (blocksize, nblocks, Block, etc.) to a new package called BlockArraysBase.jl and make a PR to BlockDiagonals.jl to use that.

[dmbates]

I was mistaken about the need to freeze the version of BlockArrays at 0.7. I believe @palday is working on removing that restriction.

The reason I said that the dependence on BlockArrays could be removed is because the full capabilities of the package are not used here. BlockArrays is just used to keep track of the blocks in the .A and .L fields in a LinearMixedModel. Using BlockArrays allows for cleaner output in a few examples showing the structure of those arrays but, other than that, it should be possible to use something like Matrix{AbstractMatrix{T}} instead.

[dmbates]

I do think, as @palday mentioned, that using BlockDiagonal in the evaluation of the ind field for amalgamated ReMats makes the structure clearer. And the BlockDiagonals package is pretty lightweight. So we could try to resolve the issues of using both BlockArrays and BlockDiagonals in the same module or just go with BlockDiagonals and replace the uses of BlockArrays.

[dlfivefifty]

I think it’s a better use of time, and more “community oriented” take make the two packages work well together. Others will run into this problem and so it would be good to solve it properly...

[dkarrasch]

Hi there, I came across this issue by pure chance, from BlockDiagonals.jl. I understood so far that you have use cases in which you have blocked arrays where the blocks have different structure: dense or sparse, and if sparse, then with different sparsity structure. In case you are "only" interested in multiplication with vectors, you may wish to consider LinearMaps.jl for that. It has support for block maps (and soon for block diagonal maps). Blocking them is all lazy, and the types of the blocks are stored in a parametric type of the BlockMap object. In terms of multiplication performance, I guess it is pretty much unbeatable, because every block (whose type is inferred) gets specialized multiplication treatment. On the "downside", LinearMap objects are not (and won't be) subtypes of AbstractArrays, so there is no getindex and setindex!, in case this should be of interest.

[dmbates]

Hi there, I came across this issue by pure chance, from BlockDiagonals.jl. I understood so far that you have use cases in which you have blocked arrays where the blocks have different structure: dense or sparse, and if sparse, then with different sparsity structure. In case you are "only" interested in multiplication with vectors, you may wish to consider LinearMaps.jl for that. It has support for block maps (and soon for block diagonal maps). Blocking them is all lazy, and the types of the blocks are stored in a parametric type of the BlockMap object. In terms of multiplication performance, I guess it is pretty much unbeatable, because every block (whose type is inferred) gets specialized multiplication treatment. On the "downside", LinearMap objects are not (and won't be) subtypes of AbstractArrays, so there is no getindex and setindex!, in case this should be of interest.

Thank you for pointing this out. The use case here is a blocked Cholesky factorization, which may or may not be a good fit for LinearMaps.jl. The matrix being factored is constructed from a fixed, blocked matrix and a variable matrix that depends upon a, generally small, vector of parameters. In the benchmarking that I have done, it is the Cholesky factorization that is taking the majority of the time, rather than the construction of the matrix to be factored. It may be possible to phrase this using LinearMaps.jl. I'll take a look and see what I can come up with.

[dmbates]

As both @dlfivefifty and @nickrobinson251 have made contributions to https://github.com/JuliaArrays/BlockArrays.jl, perhaps one of them could create a BlockArraysBase package under JuliaArrays?