JuliaGeometry / CoordinateTransformations.jl

A fresh approach to coordinate transformations...
Other
179 stars 25 forks source link

Unitful compatibility #60

Closed SebastianM-C closed 4 years ago

SebastianM-C commented 4 years ago

I did not consider <:Number because I didn't know how to do it without getting a stack overflow and also because Unitful subtypes Number and promoting to a common type is exactly the opposite of what this PR is trying to do.

I tried to cover the most common cases with the tests. Let me know what you think of this. This closes #58

andyferris commented 4 years ago

The Inetger vs AbstractFloat thing will be very helpful. I supposing there will be corner cases with Dual and so-on, but if Unitful is using Number (which I take as at least being a semi-ring!) then there's not much we can do...

SebastianM-C commented 4 years ago

There might be one or two places where you need some type (for zero(T)) but looking at that, that logic might need cleaning up in any case (these days I tend to use zero(x) where x is a value).

Do you mean to say that I should add zero for Polar and the others?

Regarding the Inetger with AbstractFloat combinations, the old implementation had some restrictions for things like

julia> x=SVector{3}(1,0,0);

julia> CylindricalFromCartesian()(x)
ERROR: MethodError: no method matching Cylindrical(::Float64, ::Float64, ::Int64)
Closest candidates are:
  Cylindrical(::T, ::T, ::T) where T at C:\Users\sebastian\.julia\packages\CoordinateTransformations\HeCSs\src\coordinatesystems.jl:84
Stacktrace:
 [1] (::CylindricalFromCartesian)(::SArray{Tuple{3},Int64,1,3}) at C:\Users\sebastian\.julia\packages\CoordinateTransformations\HeCSs\src\coordinatesystems.jl:153
 [2] top-level scope at none:0

with this PR

julia> x=SVector{3}(1,0,0);

julia> CylindricalFromCartesian()(x)
Cylindrical(r=1.0, θ=0.0 rad, z=0.0)

Regarding the conversions, there is a slight discrepancy between the unitless and unitful cases for things like Polar(1.0,2) vs Polar(1.0u"m",2).

julia> Polar(1.0,2)
Polar(r=1.0, θ=2.0 rad)

julia> Polar(1.0u"m",2)
Polar(r=1.0 m, θ=2 rad)

This is because the number type is wrapped by Quantity (more about this at https://painterqubits.github.io/Unitful.jl/stable/types/) and I was not sure how to implement this. There is something that addresses this in RecursiveArrayTools.jl, more precisely recursive_unitless_eltype which uses one in the Unitful case

https://github.com/SciML/RecursiveArrayTools.jl/blob/84f3051388de0c0ea984f34ff94a1db454bf46a8/src/utils.jl#L94

I was thinking that I could do something similar, but first I wanted to ask whether the logic for changing the type should be in the constructor or in the CartesianFromX transformations. I also wanted to ask why is important to have a uniform type inside Polar (and friends).

Also I'm not sure what happens when you combine Dual numbers and Unitful or other types, but if you have ideas for some tests I can include them.

andyferris commented 4 years ago

The constructor is the right place to deal with promotion.

Note that these packages all define promote fine.

julia> using Unitful, DualNumbers

julia> x = 1u"m"
1 m

julia> y = 2.0
2.0

julia> d = dual(0x01, 0x02)
1 + 2ɛ

julia> promote(x, y)
(1.0 m, 2.0)

julia> promote(x, d)
(1 + 0ɛ m, 1 + 2ɛ)

julia> promote(y, d)
(2.0 + 0.0ɛ, 1.0 + 2.0ɛ)

For Polar I would do this non-recursive outer constructor:

function Polar(r, theta)
    (r2, theta2) = promote(r, theta)
    return Polar{typeof(r2), typeof(theta)}(r2, theta)
end

and let people use the inner constructor to override.

SebastianM-C commented 4 years ago

Thanks for the suggestion, this is much better than my previous implementation. I implemented promotion as you said, but there is a small problem with testing. It looks like isapprox is broken in the case of dimensionless numbers (<: DimensionlessQuantity) due to a missing dispatch. Should I wait for the upstream fix or should I just test the individual components (and use ustrip) to avoid the problem?

andyferris commented 4 years ago

I think @test_broken and waiting for upstream is quite sensible.