Improve `AffineDecomposition`

[pbrehmer]

The type restriction on the coefficient_map should be lifted in order to represent affine decompositions with constant coefficients. Also, as it turns out, this seems necessary to allow for a working custom serialization of AffineDecompositions.

@mfherbst Do you think we should implement a compress(::AffineDecomposition{<:AbstractArray,<:AbstractArray}, ...) for affine decompositions with constant coefficients that returns the compressed matrix (as we shortly discussed with regards to m_reduced in the docs examples)?

It could be that this makes things more confusing, since compress would then return 3 objects (the matrix, the compressed AffineDecomposition and a AffineDecomposition with only the matrix elements, which is needed for later truncation).

• closes #40

[mfherbst]

@mfherbst Do you think we should implement a compress(::AffineDecomposition{<:AbstractArray,<:AbstractArray}, ...) for affine decompositions with constant coefficients that returns the compressed matrix (as we shortly discussed with regards to m_reduced in the docs examples)?

It could be that this makes things more confusing, since compress would then return 3 objects (the matrix, the compressed AffineDecomposition and a AffineDecomposition with only the matrix elements, which is needed for later truncation).

Generally I would indeed discourage this. On the other hand if you use named tuples this is not a problem, since then one can either choose to use this feature or not. In my opinion the deciding factor is whether this either improves convenience (in the sense that one saves a lot of repetitive code) or it improves performance. I don't think it is the latter, but it might be the former. What do you think?

[pbrehmer]

@mfherbst Do you think we should implement a compress(::AffineDecomposition{<:AbstractArray,<:AbstractArray}, ...) for affine decompositions with constant coefficients that returns the compressed matrix (as we shortly discussed with regards to m_reduced in the docs examples)?

It could be that this makes things more confusing, since compress would then return 3 objects (the matrix, the compressed AffineDecomposition and a AffineDecomposition with only the matrix elements, which is needed for later truncation).

Generally I would indeed discourage this. On the other hand if you use named tuples this is not a problem, since then one can either choose to use this feature or not. In my opinion the deciding factor is whether this either improves convenience (in the sense that one saves a lot of repetitive code) or it improves performance. I don't think it is the latter, but it might be the former. What do you think?

My Feeling is that including this option does not avoid that much code repetition. Also, I think this would be confusing to describe in the Docs examples, even more so than explaining that some AffineDecompositions have constant coefficients such that the explicit sum can be obtained by calling just (::AffineDecomposition)(). So for now, I would prefer to exclude the matrix-return option.

[mfherbst]

Ok, then let's don't. Just rebase the code and the we merge.

[pbrehmer]

Okay, I will adjust the alignment and change the type symbol to C, and resolve the conflicts later.

I would have one more thing we might add in this PR on the topic of improving AffineDecompositions. So in practice one needs a parallelized version of compress, e.g. when computing the compressed terms for a structure factor — they can be evaluated all in parallel since there is no data dependency between the terms.

Now in my private RB applications repository, I have already implemented such a method for affine decompositions with matrix-like terms:

function compress(ad::AffineDecomposition{<:Matrix}, basis::RBasis; symmetric_terms=false)
    matrixel = similar(ad.terms, Matrix{Number})
    is_nonredundant(a, b) = (size(ad.terms, 1) > size(ad.terms, 2)) ? ≥(a, b) : ≤(a, b)
    Threads.@threads for idx in CartesianIndices(matrixel)
        if is_nonredundant(first(idx.I), last(idx.I))  # Redundant upper/lower triangle
            matrixel[idx] = compress(ad.terms[idx], basis.snapshots)
        end
    end  

    Threads.@threads for idx in findall(x -> !is_nonredundant(first(x.I), last(x.I)), matrixel)
        if symmetric_terms
            matrixel[idx] = matrixel[last(idx.I), first(idx.I)]
        else
            matrixel[idx] = compress(ad.terms[idx], basis.snapshots)
        end
    end  # Use "symmetry" to fill redundant elements
    matrixel = promote_type.(matrixel)  # Promote to common floating-point type
    rbterms = map(m -> basis.vectors' * m * basis.vectors, matrixel)

    (; rb=AffineDecomposition(rbterms, ad.coefficients),
     raw=AffineDecomposition(matrixel, ad.coefficients))
end

Of course this could be generalized to any type of terms, but this particular method is especially needed at the moment. Unfortunately, this still uses the promote_type but I don't know how to avoid this at the moment.

Note that this also touches the issue #35, which we would have to rethink anyway in order to implement a parallelized compress.

[pbrehmer]

@mfherbst I think I'll leave this compress feature for now and just use the one I implemented in my (still private) applications repository. I'll just resolve the conflicts and then we can merge.