Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

dispatch on units of angle #38

Closed
mweastwood opened this issue Oct 18, 2016 · 4 comments
Closed

dispatch on units of angle #38

mweastwood opened this issue Oct 18, 2016 · 4 comments

Comments

@mweastwood
Copy link
Contributor

I am interested in adding Unitful as a dependency for a package I maintain, but it is important that I am able to dispatch on units of angle. For example I have a 3-vector type, and it would be nice to provide constructors like (without units all of the arguments are Float64 and you don't get the convenience of having both constructors):

f(x::Unitful.Length, y::Unitful.Length, z::Unitful.Length) = ...
f(r::Unitful.Length, theta::Unitful.Angle, phi::Unitful.Angle) = ...

I am also going to argue that "Angle" should be a fundamental (not derived dimension). The idea that angles are "dimensionless" seems to originate from the idea that it is a length / length (see table 3 http://physics.nist.gov/cuu/Units/units.html), but this is only correct in the small angle approximation. In general you need to use a trig function to get from an angle to something that is truly dimensionless. Or rather you need an inverse trig function to get an angle from something that is dimensionless.

Finally it would also be really cool if sin(1°) dispatched to sind(1).
(see https://github.com/Keno/SIUnits.jl/blob/d25c43483828c954cd63c19cbe3ff209a614bbd1/src/nonsiunits.jl#L15-L22)

@ajkeller34
Copy link
Collaborator

Wow, small world, someone else at Caltech found my package...

First, let me argue that no small angle approximation is required. Consider that the circumference of a circle is (2π radians)*(radius). Surely the circumference has dimensions of length, not dimensions of length*angle. The dimension of 2π radians must then be that of circumference over radius, or length/length, which is dimensionless. So, angle dimensions are kind of problematic for consistency when doing dimensional analysis, which is probably why there is no angle dimension in SI.

To do the dispatch you want, I offer two solutions:

  1. You can specify theta::Number, phi::Number and then do a dimensions check in the function body (dimension(theta) != NoDims && throw(Unitful.DimensionError())). There will be no ambiguity with the other method you defined, and in the past I have found that the compiler can totally elide the dimensions check: even though it looks like a run-time check, you shouldn't see a branch in the generated code when you pass in a valid argument. This isn't something guaranteed by Unitful, rather it just appears to work that way based on some black magic going on inside Julia. I think it works because the dimension can be decided based on the type of theta, rather than its value. In other words, there should be no performance penalty for doing this.
  2. If you really want to be very specific about the allowed types---and remember, there's no performance advantage to being specific with your argument types, you only need to be specific if you want to for dispatch reasons---you can try something like the following, which will only accept reals or dimensionless quantities that are based on reals. For the record I wouldn't do this unless you had very strong reasons, since it would unfairly exclude other potentially valid numeric types.
julia> using Unitful

julia> using Unitful: °, rad, kg

julia> typealias Angle{T} Union{T, Unitful.DimensionlessQuantity{T}}
Union{T,Unitful.Quantity{T,Unitful.Dimensions{()},U}}

julia> g{T<:Real}(::Angle{T}) = 1
g (generic function with 1 method)

julia> g(::Any) = 2
g (generic function with 2 methods)

julia> g(1°)
1

julia> g(1rad)
1

julia> g(1.0)
1

julia> g(1+2im)
2

julia> g((1+2im)rad)
2

julia> g(1kg)
2

Regarding sind, without having studied the implementation carefully it does seem like sind has been optimized in base julia but strangely when I benchmark sin(1.0°) vs. sind(1.0) I don't really see a performance difference, if anything sind(1.0) is slower somehow. Sounds suspicious so maybe I'll look into what's going on there.

@mweastwood
Copy link
Contributor Author

Wow, small world, someone else at Caltech found my package...

Hi! Julia users at Caltech are a rare breed right now.

First, let me argue that no small angle approximation is required.

For some reason I was only thinking about chords (and not arcs) when I made this statement. You are definitely correct here.

Between options 1 and 2 my preference is for option 2 because it is stylistically more Julian (although I agree that the 2 options are functionally equivalent because the compiler will delete the dead branch if you call it with the right types). However there is nothing stopping me from putting the typealias in my own package so if you are somewhat hesitant, this might be the best solution.

I was never asking about performance (I actually don't mind a runtime check) but rather being able to conveniently define functions that act on angles.

For example:

julia> using Unitful

julia> type Vec3
           x::Float64
           y::Float64
           z::Float64
       end

julia> Vec3(::Unitful.Length, ::Unitful.Length, ::Unitful.Length) = ... cartesian constructor ...

julia> Vec3(::Unitful.Length, ::Angle, ::Angle) = ... spherical constructor ...

In this example, my motivation is a convenient interface for the user to pick between the cartesian and spherical constructors. Without attaching units to the arguments I can't discern the user's intention from the types of the arguments. I would instead have to define two functions Vec3_cartesian and Vec3_spherical, which is somewhat more cumbersome.

For this reason the typealias solution is more than sufficient for me. If you don't want to add that definition to Unitful, feel free to close this issue.

Lastly I believe sind is for numerical accuracy, not speed.

julia> sind(360)
0.0

julia> sin(deg2rad(360))
-2.4492935982947064e-16

Also thanks for developing this package! I've always felt like SIUnits needed some TLC, so it's great to see Unitful with some great ideas.

@ajkeller34
Copy link
Collaborator

You're welcome! And thank you for alerting me to the purpose of sind. I will spin off into a separate issue.

Unless I misunderstand what you want, the correct method or constructor should be called even if you leave off the ::Angle or replace with ::Number, I think:

julia> using Unitful

julia> using Unitful: °

julia> type Vec3
    x::Float64
    y::Float64
    z::Float64
end

julia> Vec3(::Unitful.Length, ::Unitful.Length, ::Unitful.Length) = 1 # or ... cartesian constructor ...
Vec3

julia> Vec3(::Unitful.Length, y::Number, z::Number) = 2 # or ... spherical constructor with dim checks ... 
Vec3

julia> Vec3(0.0, 0.0, 0.0)
Vec3(0.0,0.0,0.0)

julia> Vec3(1u"m", 2u"m", 3u"m")
1

julia> Vec3(1u"m", 0°, 2°)
2

I guess I am hesitant to include an Angle definition, at least for the time being. I see the convenience, but since we agree angles have no dimensions, I'd need to come up with a fair restriction on the numeric types an angle is allowed to be, which is hard to decide in general when you interact with other packages. The typealias I provided could prevent you from playing nicely with other package-defined Number types. I don't have an example in mind but for instance in the past I've run into trouble with dual numbers (see ForwardDiff.jl) when I over-specify types.

@cmichelenstrofer
Copy link
Contributor

I know this is a very old thread, but just in case someone else finds it: checkout DimensionfulAngles.jl. It treats angles as a dimension which allows for dispatching.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants