Inner product between bases on different grids, basis transforms

Say A isa Basis and B isa Basis, where axes(A,1) != axes(B,1), then we get the following error:

julia> A'B
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match

However, if I have a function $f$ that is expanded over the basis A, $|f\rangle = \sum_i c_i|a_i\rangle$, and I want to project it onto the basis B, $|f\rangle = \sum_j d_i|b_j\rangle$, then the above basis overlap is precisely the one I need, since $d_j = \sum_i c_i \langle b_j|a_i\rangle$, or equivalently $\mathbf{d}=A^HB\mathbf{c}$.

To me it is unimportant that the axes are different, i.e. I think it is up to the user to make sure that the end result makes sense (e.g. projecting onto a "larger" basis, or maybe I am interested in the subspace, etc).

I can of course override this on a case-by-case basis, but it would be convenient if this just worked, and the only thing you'd need to provide is the calculation of basis function overlaps $\langle b_j|a_i\rangle$. Then it would be very easy to transform between e.g. B-splines and finite-elements or -differences.

What do you think?

EDIT: It seems it will be difficult to work around this on a case-by-case basis, since I cannot circumvent the dimension check easily:

julia> ApplyQuasiArray(*, A', B)
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match
Stacktrace:
 [1] _check_mul_axes(A::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, B::QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:91
 [2] check_mul_axes
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:107 [inlined]
 [3] check_applied_axes
   @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:22 [inlined]
 [4] instantiate
   @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:63 [inlined]
 [5] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:53 [inlined]
 [6] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:54 [inlined]
 [7] ApplyQuasiArray(M::LazyArrays.Applied{LazyArrays.MulStyle, typeof(*), Tuple{QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false}}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:55
 [8] ApplyQuasiArray(::Function, ::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, ::Vararg{Any})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:59

Start at the multiplication entry points shown in the trace: ApplyQuasiArray and ArrayLayouts' check_mul_axes. Review how A'B validates axes and how basis overlaps are represented. Done means a tested path for computing overlaps and transforming between bases with different grids.