Units

[Tim Holy]

Hi folks,

I just pushed a small set of utilities for dealing with variables that have
units. I have to say, Julia's design really shines for this kind of stuff...

Here are some examples of what you can do:

julia> load("units.jl")

julia> v = parse_sivalue("3.5 us")
3.5 us

julia> pshow(v)
3.5 μs
julia> T = SIUnit(Second)
SIUnit{,s}

julia> convert(T, v)
3.5e-6 s

julia> T = SIUnit(Inch)
SIUnit{,in}

julia> convert(T, v)
Not convertable
in convert at /home/tim/src/julia/extras/units.jl:217

julia> v = SIValue(Meter, 1.0)
1.0 m

julia> convert(T, v)
39.37007874015748 in

julia> fshow(v)
1.0 meter
julia> convert(SIUnit(Zetta,Meter), v)
1.0e-21 Zm

julia> fshow(convert(SIUnit(Zetta,Meter), v))
1.0e-21 zettameter


While I was comprehensive about the SI prefixes, there are many simple units
that are not in there. Add as you see fit.

And I haven't even thought about compound units, like joule = kg m^2/s^2.
Really what prompted this was parsing measurements from header files (though
obviously I got a bit carried away...), and it so happens that all of the
measurements in our text files are simple.

--Tim

[Matthias Schabel]

Tim,

(Small plug of my own work) You might want to take a look at the Boost.Units library. It is a C++ TMP library that achieves fully generic dimensional analysis at compile time with zero runtime overhead. Would be very interesting to see if the same thing can be accomplished in Julia without the enormous amount of template ugliness that it took to achieve in C++...

Matthias

> --
>
>
>

[Tim Holy]

On Thursday, July 19, 2012 04:21:11 PM Matthias Schabel wrote:
> (Small plug of my own work) You might want to take a look at the Boost.Units
> library. It is a C++ TMP library that achieves fully generic dimensional
> analysis at compile time with zero runtime overhead. Would be very
> interesting to see if the same thing can be accomplished in Julia without
> the enormous amount of template ugliness that it took to achieve in C++...

I looked, and it's very impressive. I can only imagine how hard this must have
been to achieve in C++ (although of course C++'s facility for template
metaprogamming is very useful). This Julia file is currently much simpler, with
the most significant limitation being the lack of compound unit types. To be a
little more complete, I added quite a few more elementary types, although I
didn't go so far as to add "furlong" and "league" :-).

In this Julia implementation, conversion among types should also have pretty
close to no runtime overhead: all the unit information is is encoded using
type parameters, which (at least in principle) the compiler should be able to
optimize away. The only overhead might be for the numerical factors when the
value is not known at compile time. Even that seems like it should be able to
be optimized the minimum of operations (one multiplication, except for cases
like temperature which are more complex). String parsing (for I/O), of course,
will always have some runtime overhead, but I used Dicts so at least it's
pretty efficient.

Currently the whole thing in Julia weighs in at just 326 lines. It makes heavy
use of Julia's ability to generate code programatically (so, of course, there
is some _load time_ overhead). I think in C++ terms you'd say this is one big
exercise in traits. If you're interested, the code should be pretty readable:
basically it sets up a bunch of tables with all the information you need, and
then generates functions (presumably all inline-able) that perform the
operations.

It will be interesting to see how much uglier it gets when somebody implements
support for compound types. Hopefully, that too can be implementable in a
reasonably clean fashion, but I am sure it will increase the line count
noticeably.

Oh, and I should say that I decided my original type names stank. Compared to
my first email, SIUnit is now Unit, and SIValue is now Quantity. Many of the
units were not in the SI system, so the name didn't really make sense.

--Tim

[Tim Holy]

On Monday, August 06, 2012 11:53:14 PM Gunnar Farnebäck wrote:
> Compound units can be implemented using types parameterized by integers,
> encoding the exponent of each base unit, like in this incomplete proof of
> concept:

Very cool! Presumably we'd want the SI type to be extended to all 7 of the
basic SI units, but your proof-of-concept is more readable because you focused
on the big three.

If the negative-integer dispatch issue can be fixed (on that I cannot comment),
I'd definitely merge a pull request containing this functionality.

> type SI{m, s, kg}
> value::Float
> end

I wonder if we should base this on g rather than kg? It will presumably make
the prefix extraction more straightforward. If you want to use kg in displayed
output, we could modify the show method to default to this. But I'm happy to
accept counterarguments.

The rest looks great. I'm especially impressed by this example:
> julia> 70kg*9.81m/s^2
> 686.7 N

--Tim

[Aron Ahmadia]

> The rest looks great. I'm especially impressed by this example:
>> julia> 70kg*9.81m/s^2
>> 686.7 N

+1

[Jeffrey Sarnoff]

using _s, _m, _kg for seconds, meters, kilograms .. etc makes it more readable  

[Jeffrey Sarnoff]

with numbers, that is:  75kg 75_kg

[Patrick O'Leary]

On Tuesday, August 7, 2012 4:29:14 AM UTC-5, Tim wrote:

I wonder if we should base this on g rather than kg? It will presumably make
the prefix extraction more straightforward. If you want to use kg in displayed
output, we could modify the show method to default to this. But I'm happy to
accept counterarguments.

The main couterargument is that the kilogram, not the gram, is the SI base unit, so derived units are already expressed in those terms.

This is very cool!

[Stefan Karpinski]

This is ingenious. Slow clap.

On Tue, Aug 7, 2012 at 2:53 AM, Gunnar Farnebäck <gun...@lysator.liu.se> wrote:

And I haven't even thought about compound units, like joule = kg m^2/s^2.

Compound units can be implemented using types parameterized by integers, encoding the exponent of each base unit, like in this incomplete proof of concept:

type SI{m, s, kg}
  value::Float
end

*{m1, m2, s1, s2, kg1, kg2}(x::SI{m1, s1, kg1}, y::SI{m2, s2, kg2}) = SI{m1 + m2, s1 + s2, kg1 + kg2}(x.value * y.value)
*{m, s, kg}(x::Number, y::SI{m, s, kg}) = SI{m, s, kg}(x * y.value)
*{m, s, kg}(x::SI{m, s, kg}, y::Number) = SI{m, s, kg}(x.value * y)

^{m, s, kg}(x::SI{m, s, kg}, n::Int) = SI{n * m, n * s, n * kg}(x.value^n)

/{m1, m2, s1, s2, kg1, kg2}(x::SI{m1, s1, kg1}, y::SI{m2, s2, kg2}) = SI{m1 - m2, s1 - s2, kg1 - kg2}(x.value / y.value)
/{m, s, kg}(x::Number, y::SI{m, s, kg}) = SI{-m, -s, -kg}(x / y.value)
/{m, s, kg}(x::SI{m, s, kg}, y::Number) = SI{m, s, kg}(x.value / y)

+{m, s, kg}(x::SI{m, s, kg}, y::SI{m, s, kg}) = SI{m, s, kg}(x.value + y.value)
+{m1, m2, s1, s2, kg1, kg2}(x::SI{m1, s1, kg1}, y::SI{m2, s2, kg2}) = throw("addition of incompatible units")

-{m, s, kg}(x::SI{m, s, kg}, y::SI{m, s, kg}) = SI{m, s, kg}(x.value - y.value)
-{m1, m2, s1, s2, kg1, kg2}(x::SI{m1, s1, kg1}, y::SI{m2, s2, kg2}) = throw("subtraction of incompatible units")

typealias Meter    SI{ 1,  0,  0}
typealias Second   SI{ 0,  1,  0}
typealias Kilogram SI{ 0,  0,  1}
typealias Hertz    SI{ 0, -1,  0}
typealias Newton   SI{ 1, -2,  1}
typealias Pascal   SI{-1, -2,  1}
typealias Joule    SI{ 2, -2,  1}
typealias Watt     SI{ 2, -3,  1}

m  = Meter(1.0)
s  = Second(1.0)
kg = Kilogram(1.0)
Hz = Hertz(1.0)
N  = Newton(1.0)
Pa = Pascal(1.0)
J  = Joule(1.0)
W  = Watt(1.0)

show(io, x::Meter)    = print(io, "$(x.value) m")
show(io, x::Second)   = print(io, "$(x.value) s")
show(io, x::Kilogram) = print(io, "$(x.value) kg")
show(io, x::Hertz)    = print(io, "$(x.value) Hz")
show(io, x::Newton)   = print(io, "$(x.value) N")
show(io, x::Pascal)   = print(io, "$(x.value) Pa")
show(io, x::Joule)    = print(io, "$(x.value) J")
show(io, x::Watt)     = print(io, "$(x.value) W")

function unit_string(s, n)
  if n == 0
    return ""
  elseif n == 1
    return s
  else
    return "$s^$n"
  end
end

function show{m, s, kg}(io, x::SI{m, s, kg})
  units = [unit_string("m", m), unit_string("s", s), unit_string("kg", kg)]
  units = filter(str -> str != "", units)
  print(io, "$(x.value) " * join(units, "⋅"))
end

With this you can do things like
julia> 5kg*3s*7m^2
105.0 m^2⋅s⋅kg

julia> 70kg*9.81m/s^2
686.7 N

However, dispatch on negative integer parameters doesn't seem to work properly at the moment.
julia> 70kg*(9.81m/s^2)
no method *(SI{0,0,1},SI{1,-2,0})
 in method_missing at base.jl:70

--
 
 
 

[Jeffrey Sarnoff]

Patrick is correct -- it is important to make the kilogram primary and define the gram in terms of the kilogram.

The definition of the kilogram and a few friends are likely to change -- the definitions are to be revisited in 2014.

( e.g.http://www.nist.gov/pml/newsletter/siredef.cfm http://www.sciencedaily.com/releases/2011/09/110929074207.htm )

[Tim Holy]

Oh, I'm well aware that kg is the official base SI unit. However, if you look at
our existing extras/units.jl implementation, which is based on gram rather
than kilogram, you'll see that it's super-easy to autogenerate everything
because the prefix is separable from the base unit. However, if you standardize
on kg, then suddenly for mass milli = 10^-6, whereas for everything else
milli=10^-3.

Naturally, with a little bit of extra code anything can be done. I was just
wondering which way would be easier.

Best,
--Tim

[Jeffrey Sarnoff]

Ahh, so it has been in extra/units.jl.   Really, we need to change that --
for a technical/scientific computing platform, it sends the wrong message.
And it may gum up your stuff. CGS units shouldn't be rough after the SI.

Maybe this will ease the pain (:

type SI{m, s, kg}

  value::Float64
  SI(f::Float64) = new(f)
  SI(f::Float32) = new(float64(f))
  SI(i::Int32) = new(float64(i*1.0))
  SI(i::Int64) = new(float64(i*1.0))
end

typealias Meter    SI{ 1,  0,  0}
typealias Second   SI{ 0,  1,  0}
typealias Kilogram SI{ 0,  0,  1}
typealias Hertz    SI{ 0, -1,  0}
typealias Newton   SI{ 1, -2,  1}
typealias Pascal   SI{-1, -2,  1}
typealias Joule    SI{ 2, -2,  1}
typealias Watt     SI{ 2, -3,  1}

_m  = Meter(1)
_s  = Second(1)
_kg = Kilogram(1)
_Hz = Hertz(1)
_N  = Newton(1)
_Pa = Pascal(1)
_J  = Joule(1)
_W  = Watt(1)


function (^)(x::SI,y::Integer)
  val = x.value
  si  = map( a->(a*y), typeof(x).parameters )
  SI{si[1],si[2],si[3]}(val)
end


function (*)(x::SI,y::SI)
  val = x.value * y.value)
  val = x.value
  si  = map( (a,b)->(a+b), typeof(x).parameters, typeof(y).parameters )
  SI{si[1],si[2],si[3]}(val)
end

function (/)(x::SI,y::SI)
  val = x.value * y.value
  si  = map( (a,b)->(a-b), typeof(x).parameters, typeof(y).parameters )
  SI{si[1],si[2],si[3]}(val)
end

function (*)(x::Number,y::SI)
    val = x * y.value)
  val = x.value
    typeof(y)(val)
end

function (/)(x::SI,y::Number)
    val = x.value / y
    typeof(x)(val)
end

show(io, x::Meter)    = print(io, "$(x.value) m")
show(io, x::Second)   = print(io, "$(x.value) s")
show(io, x::Kilogram) = print(io, "$(x.value) kg")
show(io, x::Hertz)    = print(io, "$(x.value) Hz")
show(io, x::Newton)   = print(io, "$(x.value) N")
show(io, x::Pascal)   = print(io, "$(x.value) Pa")
show(io, x::Joule)    = print(io, "$(x.value) J")
show(io, x::Watt)     = print(io, "$(x.value) W")

function show(io, x::SI)
    fns(a,b) = ((a>0) ? "$(b)^$(a) " : (a<0) ? "$(b)^($(a)) " : "")
)
  val = x.value
    m_pow,t_pow,kg_pow = typeof(x).parameters
    m_pow = fns(m_pow,"m")
    t_pow = fns(t_pow,"t")
    kg_pow = fns(kg_pow,"kg")
    units = "$(m_pow)$(t_pow)$(kg_pow)"
    if (length(units) > 0) units = "[$(units[1:end-1])]" end
 
    print(io, "$(x.value)$(units)")
end


julia> 70_kg*9.81_m/_s^2
686.7 N

julia> 70_kg*(9.81_m/_s^2)  # negative integers no problem
686.7 N

[Tim Holy]

On Tuesday, August 07, 2012 03:42:37 PM Jeffrey Sarnoff wrote:
> Ahh, so it has been in extra/units.jl. Really, we need to change that --
> for a technical/scientific computing platform, it sends the wrong message.
> And it may gum up your stuff. CGS units shouldn't be rough after the SI.

My only point was that the issue of internal representation could be a
different matter from what you show the world, _if_ it makes your life easier.
I was thinking mostly of the string parsing and type generation, which is the
bulk of what's in extras/units.jl---the main aim of that work was to be able
to parse text files that have tables of values that may or may not be
associated with strings that denote units. Turned out to be a super-easy
problem, given Julia's wonderful ability to generate code inside for loops
:-), but of course autogeneration is a little easier if you can count on the
separability (in the mathematical sense) of the prefix and the base unit.

But as you point out, there are potential problems with that approach---having
internal and external representations differ can be confusing. I think someone
sitting down and actually implementing it (the full version, not just the
proof-of-principle) will make the best path clear. The bottom line is that I
don't care how it gets done, I'll be happy to merge whatever works and people
are happy with.

Somebody add this functionality to units.jl & generate a pull request! It will
be a very nice addition indeed.

--Tim