*(::AbstractTriangular, ::DimArray) discards dimensions of triangular matrix

[sethaxen]

An UpperTriangular(::AbstractDimMatrix) has DimUnitRange axes. However, multiplying it by a DimVector results in compatibility of the inner dimension not being validated and of the outer dimension being discarded, e.g.

julia> using DimensionalData, LinearAlgebra, Test

julia> A = UpperTriangular(DimMatrix(randn(5, 5), (X, Y)));

julia> B = DimVector(randn(5), Y);

julia> axes(A)
(X(Base.OneTo(5)), Y(Base.OneTo(5)))

julia> A * B  # should this have dim X?
┌ 5-element DimArray{Float64, 1} ┐
├────────────────────────── dims ┤
  ↓ AnonDim Base.OneTo(1)
└────────────────────────────────┘
   1     -0.192035
 #undef   3.22561
 #undef   0.765636
 #undef   2.12322
 #undef   0.239097

julia> A' * B  # should this error?
┌ 5-element DimArray{Float64, 1} ┐
├────────────────────────── dims ┤
  ↓ AnonDim Base.OneTo(1)
└────────────────────────────────┘
   1     0.00633928
 #undef  0.147744
 #undef  0.0511726
 #undef  1.29501
 #undef  1.95363

I wonder if some special-casing can be performed on the following lines to check if the AbstractMatrix has DimUnitRange axes, and if so, to use those:

https://github.com/rafaqz/DimensionalData.jl/blob/6db30de4b2e1fc7f8611b7e1dc3f89dc02c78598/src/array/matmul.jl#L114-L124

Also, for some reason the multiplication is performed twice here, and the AnonDim is given a lookup with a different length than the dimension length. I wonder if the following is more sensible:

function DimensionalData._rebuildmul(A::AbstractMatrix, B::AbstractDimVector)
    ax1, ax2 = axes(A)
    if ax2 isa Dimensions.DimUnitRange
        isstrict = DimensionalData.strict_matmul()
        Dimensions.comparedims(Dimensions.dims(ax2), first(Dimensions.dims(B)); 
            order=isstrict, val=isstrict, length=false
        )
    end
    newdata = A * parent(B)
    if newdata isa AbstractArray
        out_dim = if ax1 isa Dimensions.DimUnitRange
            Dimensions.dims(ax1)
        else
            Dimensions.AnonDim(ax1)
        end
        DimensionalData.rebuild(B, newdata, (out_dim,))
    else
        newdata
    end
end

With this method, we get

julia> @inferred A * B
┌ 5-element DimArray{Float64, 1} ┐
├────────────────────────── dims ┤
  ↓ X
└────────────────────────────────┘
 -0.192035
  3.22561
  0.765636
  2.12322
  0.239097

julia> A' * B
ERROR: DimensionMismatch: X and Y dims on the same axis.
...

[rafaqz]

That looks sensible to me, I think matmuls are from before we had DimUnitRange

[sethaxen]

Would it make sense to have the following dims fallback?

dims(x::AbstractArray) = dims(axes(x))

With this, the following works:

julia> A = UpperTriangular(DimMatrix(randn(5, 5), (X, Y)));

julia> dims(A)
(↓ X, → Y)

julia> dims(DimArray(fill(1), ()))
()

julia> dims(randn(5, 5))

so its only behavior change seems to be to expand dims support to non-AbstractDimArrays that have dimensional indices. Then e.g. _comparedims_mul could just return true if dims(a) or dims(b) are nothing and otherwise call comparedims on the dims, and then this method is generalized to support non-AbstractDimArrays.

[rafaqz]

Yeah thats a great solution. But should we put it in the breaking branch just in case?