|
| 1 | +module UnitfulLinearAlgebra |
| 2 | + |
| 3 | +using Unitful:@u_str, FreeUnits, unit |
| 4 | +using LinearAlgebra: dot |
| 5 | + |
| 6 | +# We want TagsAlongAxis to act like an AbstactVector |
| 7 | +# TagsAlongAxis{T} <: AbstractVector{Unitful.FreeUnits} means that *instances* of TagsAlongAxis act like an AbstractVector. Not what we want. |
| 8 | +struct TagsAlongAxis{T} |
| 9 | +end |
| 10 | + |
| 11 | +#TODO also define `lastindex` |
| 12 | +#TODO restrict T to be tuple? |
| 13 | +Base.getindex(taa::Type{TagsAlongAxis{T}}, i)where {T}= T[i] |
| 14 | +Base.iterate(taa::Type{TagsAlongAxis{T}}, args...)where {T}=iterate(T, args...) |
| 15 | +Base.length(taa::Type{TagsAlongAxis{T}})where {T}=length(T) |
| 16 | + |
| 17 | +import LinearAlgebra |
| 18 | + |
| 19 | +# generic fallback won't work unless we define `iterate` on `typeof(u"m")`. |
| 20 | +# u"m" is an object which supports only `:*`. `u"m" + u"cm"` has no meaning. |
| 21 | +# to see if two `TagsAlongAxis` are compatible for contraction, give them units and check that. |
| 22 | +function LinearAlgebra.dot(v1::Type{TagsAlongAxis{UNITS1}}, v2::Type{TagsAlongAxis{UNITS2}})where {UNITS1, UNITS2} |
| 23 | +@assertlength(UNITS1)==length(UNITS2) |
| 24 | + u1_one=map(u->true* u, UNITS1) |
| 25 | + u2_one=map(u->true* u, UNITS2) |
| 26 | + |
| 27 | + s= u1_one[1]* u2_one[1] |
| 28 | +for i=2:length(UNITS1) |
| 29 | + s+= u1_one[i]* u2_one[i] |
| 30 | +end |
| 31 | +returnunit(s) |
| 32 | +end |
| 33 | + |
| 34 | +#TODO consider using (Holy-) Traits |
| 35 | +Multipliable= Union{Number, FreeUnits} |
| 36 | + |
| 37 | +function Base.:*(a::Multipliable, v::Type{TagsAlongAxis{UNITS}})where {UNITS} |
| 38 | +return TagsAlongAxis{map(u-> a* u, UNITS)} |
| 39 | +end |
| 40 | + |
| 41 | + |
| 42 | +function Base.:*(v::Type{TagsAlongAxis{UNITS}}, a::Multipliable)where {UNITS} |
| 43 | +# field multiplication should be commutative, but maintain multipliation order anyway. |
| 44 | +return TagsAlongAxis{map(u-> u* a, UNITS)} |
| 45 | +end |
| 46 | + |
| 47 | +struct TagsOuterProduct{T} |
| 48 | +end |
| 49 | + |
| 50 | +#TODO types vs singletons. Can dispatch on types. |
| 51 | +function Base.getindex(top::Type{TagsOuterProduct{Tuple{TAA}}}, ind)where {TAA} |
| 52 | +returngetindex(TAA, ind) |
| 53 | +end |
| 54 | + |
| 55 | + |
| 56 | +#TODO define recursively, dispatch on Tuple{TAA1, VarArgs} |
| 57 | +function Base.getindex(top::Type{TagsOuterProduct{Tuple{TAA1, TAA2}}}, ind1, ind2)where {TAA1, TAA2} |
| 58 | +returngetindex(TAA1, ind1)*getindex(TAA2, ind2) |
| 59 | +end |
| 60 | + |
| 61 | + |
| 62 | +#TODO restrict TT to TT<:TagsOuterProduct |
| 63 | +struct TaggedTensor{T, N, A<:AbstractArray, TT}<:AbstractArray{T, N} |
| 64 | + x::A |
| 65 | +end |
| 66 | + |
| 67 | +#TODO put in inner constructor |
| 68 | +functionTTensor(tag, x) |
| 69 | +TaggedTensor{eltype(x), ndims(x), typeof(x), tag}(x) |
| 70 | +end |
| 71 | + |
| 72 | +get_tag(::TaggedTensor{T, N, A, TAG})where{T, N, A, TAG}= TAG |
| 73 | + |
| 74 | +Base.size(t::TaggedTensor)=size(t.x) |
| 75 | +Base.getindex(t::TaggedTensor, ind...)= t.x[ind...]*get_tag(t)[ind...] |
| 76 | +#get_freeunits(quantity::Unitful.Quantity{T, Dim, Unit}) = Unit |
| 77 | + |
| 78 | +# matrix-vector product |
| 79 | +function Base.:*(m::Type{TagsOuterProduct{Tuple{mTAA1, mTAA2}}}, v::Type{TagsOuterProduct{Tuple{vTAA}}})where {mTAA1, mTAA2, vTAA} |
| 80 | + taa= mTAA1*dot(mTAA2, vTAA) |
| 81 | +return TagsOuterProduct{Tuple{taa}} |
| 82 | +end |
| 83 | + |
| 84 | +# matrix-matrix product |
| 85 | +function Base.:*(m::Type{TagsOuterProduct{Tuple{m1TAA1, m1TAA2}}}, v::Type{TagsOuterProduct{Tuple{m2TAA1, m2TAA2}}})where {m1TAA1, m1TAA2, m2TAA1, m2TAA2} |
| 86 | +return m1TAA1*dot(m1TAA2, m2TAA1)* m2TAA2 |
| 87 | +end |
| 88 | + |
| 89 | +function Base.:*(t1::TaggedTensor, t2::TaggedTensor) |
| 90 | + TAG1=get_tag(t1) |
| 91 | + TAG2=get_tag(t2) |
| 92 | + TAG= TAG1* TAG2 |
| 93 | + t= t1.x* t2.x |
| 94 | + |
| 95 | +returnTTensor(TAG, t) |
| 96 | +end |
| 97 | + |
| 98 | +end# module |