Supporting Unitful LinearMaps

[JeffFessler]

This is my WIP.
Here's the main issue:
if A is a LM with elements that are integers with units seconds,
and B is LM with elements that are Float32 with units meters,
then the appropriate eltype for the composite C = A * B cannot be computed by promote_type (AFAIK) but rather is something likeeltype(oneunit(eltype(A)) * oneunit(eltype(B))
which in this case would be Float32 with units meter-seconds.
In my other packages I have managed to use oneunit like this without needing to add Unitful as a package dependency. I hope for the same approach here.

In the above example, Float64 with units meters-seconds would also be a legit eltype for C and I think the promote_type check can handle that.

The existing code forces the caller to determine the eltype of the product.
Why not provide a convenience constructor?

The changes here pass all existing tests except for one type inference.
Before I chase that down I wanted to get your general reaction.
If you are open to something along these lines then I will add Unitful tests.
If you have ideas about the type inference issue, comments welcome!

[JeffFessler]

Could eventually address #195

[dkarrasch]

Thanks @JeffFessler for starting this. As you may have noticed, I started a similar effort in LinearAlgebra, so this is very welcome. It's perhaps a matter of taste, but I like the current behaviour where we determine the new eltype at the "implicit" constructor *. It has a bit of code repetition smell, but it's very clear and has the other advantage that only two (el)types are involved. For type inference, this might be easier than going down the reduction path with prod. I'm curious how the multiplication code will work, because it currently tries to reuse temporary memory. But when all linear maps have units, then each intermediate vector will have a different unit. Once we have worked out/agreed on a design pattern, then this pattern should probably be transferred to the Kronecker stuff (and even further to the KhatriRao maps). I'm very excited about this!

[JeffFessler]

As you may have noticed, I started a similar effort in LinearAlgebra

Actually I had not noticed, but I am very glad to hear it!

and has the other advantage that only two (el)types are involved

Yeah, I kind of thought you might say that and I'm fine with retaining the small repetition to help with type inference. In some of my other code I had to tell the compiler what the return type would be to get type 'inference' to work and I hope to avoid that here
https://github.com/JuliaImageRecon/ImagePhantoms.jl/blob/0c995ab8db32b67a734f31edf72821c89400b64e/src/object3.jl#L182

curious how the multiplication code will work

I suspect we may have to use reinterpret as suggested here for fft.

Let me fix up the constructors first and then take a stab at multiply.

[JeffFessler]

@dkarrasch I think I have adopted all your suggestions and the existing tests pass fine with Julia 1.8.3 on my computer, but one of the type inferences fails on 1.6 in the CI process in a linearcombination test. Is it possible that the Julia 1.6 compiler is simply not as good at type inference? If so, personally I would just abandon 1.6 because I don't want to put time into supporting old versions. But I suspect you want to keep it. Would you mind taking a look at the type inference that failed and see if you can identify the issue? To my eyes it seems unrelated to this PR, but probably I'm missing something :)

Update: now 1.6 seems to be passing too, so perhaps your merge of master did some magic - thanks for that!

[JeffFessler]

Update: I've made progress on the products, but there is one function that will be a sticking point.
For most cases the source,dest recycling in the composite mul! works fine with reinterpret,
but I could not get it to work for one test case where one of the maps is complex and earlier ones are not.
The problem is that reinterpret from, say, complex vector to a real vector doubles the length.
So I just modified it to allocate a new dest every time, which is what _resize does anyway in cases where the sizes don't match.

I think I am partly responsible for encouraging the temporary vectors source,dest that are recycled, recall:
#90
One option here would be to keep the existing method for Real,Complex (and quaternion?) cases and use the new (always allocating) method for general Number types (which would include Unitful).
All of the tests pass on 1.8.3 on my Mac. So it's ready for more review from you.

[codecov]

Codecov Report

Patch coverage: 100.00% and project coverage change: -0.63 ⚠️

Comparison is base (65eb422) 99.68% compared to head (2ccfb0c) 99.05%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #196      +/-   ##
==========================================
- Coverage   99.68%   99.05%   -0.63%     
==========================================
  Files          22       22              
  Lines        1591     1592       +1     
==========================================
- Hits         1586     1577       -9     
- Misses          5       15      +10     

Impacted Files	Coverage Δ
src/LinearMaps.jl	100.00% <100.00%> (ø)
src/composition.jl	92.42% <100.00%> (-7.58%)	⬇️
src/scaledmap.jl	100.00% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[dkarrasch]

This is nice. I think we need unitful test cases where we would want memory reuse in compositions to work. For instance, I noticed a special pattern here in the _compositemul! loop:

dest = similar([], _multype(A.maps[n], source), (size(A.maps[n], 1), size(source)[2:end]...))
_unsafe_mul!(dest, A.maps[n], source)
# dest, source = source, dest # alternate dest and source
dest, source = [], dest

The last line should be simply source = dest, because dest is allocated in the next loop, so we better don't allocate a Vector{Any} and confuse type inference completely. Moreover, the lines above can be equivalently reduced to source = A.maps[n] * source, which means we don't reuse anything and such code is better handled by foldr(*, reverse(A.maps), init=x) as we already have for two-arg CompositeMaps. But for classic cases I think we do want to keep the memory reuse and the source-dest-swapping. So we need to draw the line between the classic cases and those that we would be fine if each factor produced its own return array. One such classic case could be where the overall eltype coincides with each single map eltype (which excludes units).

[JeffFessler]

I hope to get back to this in a couple weeks...