DynamicQuantities.jl

[odow]

Unitful has a lot of quirks.

This new package seems a lot simpler: https://discourse.julialang.org/t/ann-dynamicquantities-jl-type-stable-physical-quantities/99963

The author's other work is very solid, so it's a new package, but worth exploring.

[trulsf]

I just had a look at it, but having just dimensions and not units, you seem to loose a bit of flexibility with regard to scaling and choice of units. As far as I understood it you have to have your quantities in the base SI units. You would need to involve Unitful in some manner anyway to get the parsing/conversion between different units.

[odow]

Yeah you'd need a normalization pass that converted everything just before adding to the JuMP model. So I guess something at runtime instead of what Unitful does at compile time.

[trulsf]

I did a quick test of using DynamicalQuantities as an extension of JuMP. The code can be found here: https://github.com/trulsf/JuMP.jl/tree/tf/dyn_quant

Currently, only rudimentary support for LP/MIP, but it seems to run ok and the code is much simpler (mainly due to the fact that we do not need a wrapper for variables and constraints, but can include them directly in a DynamicQuantities.Quantity.

The following example show how it can be used:

using Unitful
using DynamicQuantities
using JuMP
using HiGHS

const DQ = DynamicQuantities

m = Model()

@variable(m, x >= 0, Dimensions(length = 1, time = -1))
@variable(m, y >= 0, convert(Dimensions, Unitful.dimension(u"m/s")))
@variable(m, z, Bin, Dimensions())

max_speed = DQ.Quantity(4, length = 1, time = -1)
@constraint(m, x + y  <= max_speed * z)

obj = @objective(m, Max, x + 3y)

set_optimizer(m, HiGHS.Optimizer)

optimize!(m)

println("obj = $(value(obj))")

println("x = $(value(x))")
println("y = $(value(y))")
println("z = $(value(z))")

If you think this is something to explore further, I can open a PR on JuMP for discussion or I can move it to a separate package if that is more suitable.

[odow]

The lack of scaling seems an issue, because that's the main source of bugs for JuMP users. And we don't want people using the base units if it's going to cause numerical issues for the solver.

Perhaps it just needs a scale dimension? for 10^x? Then km = Dimensions(length = 1, scale = 3).

I don't know if we need to add to JuMP immediately. I think it's probably better if we wait a few months for other packages to adopt (or not) it first to see how things pan out.

[trulsf]

I did another test with the recent version of DynamicQuantities (v0.4) and used a Quantity to represent the units of a variable instead of Dimensions. The syntax even without Unitful is now easier. The relative scaling of different variables are then handled, but all constraints are still in SI base units. If there was an option in the package for setting the scale of the base units to use, I guess one are getting close to something acceptable.

The problem currently is that values for the variables are returned in the base units, regardless of the scale used for them in the model. I guess this will require a version with a scaled quantity (as you have been discussing).

The following example shows the results of running the latest version

julia> using DynamicQuantities, JuMP, HiGHS

julia> m = Model(HiGHS.Optimizer);

julia> set_silent(m)

julia> @variable(m, x >= 0, u"m/s")
x [1.0 m s⁻¹]

julia> @variable(m, y >= 0, u"km/h")
y [0.2777777777777778 m s⁻¹]

julia> @variable(m, z, Bin, Dimensions())
z [1.0 ]

julia> max_speed = 4u"km/s"
4000.0 m s⁻¹

julia> @constraint(m, x + y  <= max_speed * z)
x + 0.2777777777777778 y - 4000 z <= 0 [m s⁻¹]

julia> obj = @objective(m, Max, x + 3y)
0.8333333333333334 y + x [m s⁻¹]

julia> optimize!(m)

julia> println("obj = $(value(obj))")
obj = 12000.0 m s⁻¹

julia> println("x = $(value(x))")
x = 0.0 m s⁻¹

julia> println("y = $(value(y))")
y = 4000.0 m s⁻¹

julia> println("z = $(value(z))")
z = 1.0

[MilesCranmer]

Two comments:

1. Not sure if this helps at all, but one thing I have found useful in making the to-be-implemented QuantityArray type is the following pattern:

f(v, units::Quantity) = (v .* ustrip(units), dimension(units))

Which gives you the following:

julia> f(ones(3), u"km/s")
([1000.0, 1000.0, 1000.0], m s⁻¹)

In other words, rather than storing the scale in the units, you just multiply any scaling outside of the SI base units into the expression.

So in e.g.,

julia> @variable(m, y >= 0, u"km/h")

it could be automatically converted into

julia> @variable(m, (y/0.28) >= 0, u"m/s")

I'm still not sure I fully understand the desired implementation, but perhaps with the new AbstractQuantity interface in https://github.com/SymbolicML/DynamicQuantities.jl/pull/24, you can set up type-based scaling, like:

struct ScaledQuantity{T,D<:AbstractDimensions,S} <: AbstractQuantity{T,D}
    value::T
    dimensions::D
end
ScaledQuantity(v::T, d::D, ::Val{S}=Val(1.0)) where {T,D<:AbstractDimensions,S} = ScaledQuantity{T,D,S}(v, d)
ScaledQuantity(q::Quantity{T}) where {T} = ScaledQuantity(one(T), dimension(q), Val(ustrip(q)))
Base.show(io::IO, q::ScaledQuantity{T,D,S}) where {T,D,S} = print(io, ustrip(q), " [", S, " ", dimension(q), "]")

where the S is some sort of scaling inherent to the type. You also need to modify the constructor function to propagate S:

function DynamicQuantities.constructor_of(::Type{SQ}) where {T,D,S,SQ<:ScaledQuantity{T,D,S}}
    return (args...) -> ScaledQuantity(args..., Val(S))
end

Then you could create quantities with, e.g.,

julia> ScaledQuantity(u"km/h")
1.0 [0.2777777777777778 m s⁻¹]

julia> typeof(ScaledQuantity(u"km/h"))
ScaledQuantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}, 0.2777777777777778}

julia> ScaledQuantity(u"km/h") * 10.3
10.3 [0.2777777777777778 m s⁻¹]

Though you would need to overload some functions to throw an error if S is different (if that is your desired use case).

[trulsf]

Thanks for the feedback and the tips - as well as the great work with the DynamicQuantities package.

The current implementation is based on an approach similar to your first pattern, i.e. we use a JuMP-variable as the value of a quantity by stripping off units and adding the dimension

function JuMP.add_variable(
    model::JuMP.Model,
    x::_DQVariable,
    name::String,
) 
    variable = JuMP.add_variable(model, x.variable, name)
    return DQ.Quantity(DQ.ustrip(x.unit) * variable, DQ.dimension(x.unit)) 
end

As stated above this works quite smoothly, but what I would like is to have the solution value in the same units as specified for the variable. If the variable is given as

@variable(m, x, u"kW*h") 

I want the solution value to be given with the same units - not as some value times m² kg s⁻² - which becomes large and hard to relate to. I will look a bit further into the option of using a specialized quantity as in your second example.

[MilesCranmer]

I want the solution value to be given with the same units - not as some value times m² kg s⁻² - which becomes large and hard to relate to. I will look a bit further into the option of using a specialized quantity as in your second example.

This is indeed tricky, as DynamicQuantities.jl converts everything into the same base units (defined by the particular AbstractDimensions). This is done for type stability, whereas Unitful.jl has a different type for any combination of units.

I see a few options for solving this:

1. You could perhaps try to define an PreferredDimensions <: AbstractDimensions based off of the variables the user adds to a particular model. But if the user passes, e.g., u"W" for one variable, and u"kW" for the other, you would have to make some sort promotion (or just throw an error). This solution is a bit trickier though, especially if you want to have the variables displayed as derived SI units rather than base units.

2. You could record the string the user passes, and store (1) the raw string, and (2) the Quantity. Then, when printing the solution, you could factor in the Quantity, divide out the units, and print the raw string.

e.g.,

@variable(m, x, "kW*h")

would turn into

raw_string = "kW*h"
units = uparse(raw_string)
value = units * x

...

solution = ...

println(solution / units, raw_string)

3. Maybe we could just add (2) to DynamicQuantities.jl directly, so that you could specify the preferred representation. In keeping with the type stability theme, perhaps this could be stored as a value rather than a type.

[trulsf]

I have done some more testing and experimenting based upon your ideas and possible options. Since I want to allow the user to combine variables with different scale of the same dimension, I have not included the scale as part of the dimension, but introduced it as a separate field of an adjusted quantity

struct QuantityScale{T,D} <: DQ.AbstractQuantity{T,D}
    value::T
    dimensions::D
    scale::Float64
end

QuantityScale(value, dim) = QuantityScale(value, dim, 1.0)

function Base.show(io::IO, qex::QuantityScale)
    print(io, "$(DQ.ustrip(qex) / qex.scale) [$(qex.scale) $(DQ.dimension(qex))]")
end

When creating new optimization variable, we scale the variable and keep the original scaling

function JuMP.add_variable(
    model::JuMP.Model,
    x::_ScaledVariable,
    name::String,
) 
    variable = JuMP.add_variable(model, x.variable, name)
    scale = DQ.ustrip(x.unit)
    return QuantityScale(scale * variable, DQ.dimension(x.unit), scale) 
end

When variables are combined into expressions and constraints, the scaling will always be kept at 1.0.

julia> @variable(m, x >= 0, u"m/s")
x [m s⁻¹]

julia> x.scale
1.0

julia> @variable(m, y >= 0, u"km/h")
y [0.2777777777777778 m s⁻¹]

julia> y.scale
0.2777777777777778

julia> z = x + y
x + 0.2777777777777778 y [m s⁻¹]

julia> z.scale
1.0

After optimization the optimal value is returned as a QuantityScale according to its original scaling

julia> println("y = $(value(y))")
y = 14400.0 [0.2777777777777778 m s⁻¹]

It is difficult to store the preferred unit (e.g. kW*h), as we only get the Quantity when creating the variable. I have experimented a little bit on simplifying scale/units when outputting (i.e. reverse engineering the original units).

If you want to have a look at the full implementation, it is available here https://github.com/trulsf/JuMP.jl/blob/tf/dq_scaled/ext/units_new.jl

[MilesCranmer]

I drafted some new functionality that I think will be particularly useful for your desired feature in the verbatim-units branch.

With the new AbstractDimensions struct, I declared a UnitDimensions that essentially has one dimension for every single unit symbol. This means the regular pretty-printing will work. It also lets you lazily convert to SI units.

For example:

julia> using DynamicQuantities

julia> uparse_symb = DynamicQuantities.Units.UnitSymbols.uparse
uparse (generic function with 1 method)

julia> q = uparse_symb("nm * kW * h")
1.0 nm kW h

julia> q^2
1.0 nm² kW² h²

Then we can convert to SI at the end:

julia> convert(Quantity{Float64,Dimensions}, q^2)
1.296e-5 m⁶ kg² s⁻⁴

The downside is that when you are performing computation with this dimension space, it has to do 10x more computations, as there are 10x more fields to propagate around. It also increases startup time quite a bit due to my choice to generate the struct from all the declared units. Fixed by lazily creating these symbols.

But, it is good if you want to display the user's entered units which seems like what you are looking for.

Thoughts?

[trulsf]

Thanks for following up on this! As you say, my primary use for this is to display variable values in the same units as declared by the user. With these additions I can store the original unit and use it for display purposes

struct ScaledQuantity{T,D} <: DQ.AbstractQuantity{T,D}
    value::T
    dimensions::D
    scale::Float64
    unit::Union{Nothing, DQ.Quantity}
end

ScaledQuantity(value, dim) = ScaledQuantity(value, dim, 1.0, nothing)

function Base.show(io::IO, q::ScaledQuantity)
    expr = DQ.ustrip(q)
    if isnothing(q.unit)
        print(io, "$(expr / q.scale) [$(DQ.dimension(q))]") 
        return
    end   
    print(io, "$(expr / q.scale) [$(DQ.dimension(q.unit))]")
end

And then convert to base SI units when calculating scale for the variable

function JuMP.add_variable(
    model::JuMP.Model,
    x::_ScaledVariable,
    name::String,
) 
    variable = JuMP.add_variable(model, x.variable, name)
    q = convert(DQ.Quantity{Float64,DQ.Dimensions}, x.unit)
    scale = DQ.ustrip(q)
    return ScaledQuantity(scale * variable, DQ.dimension(q), scale, x.unit) 
end

Running this produces the desired output

julia> usym = DynamicQuantities.Units.UnitSymbols.uparse
uparse (generic function with 1 method)

julia> m = Model(HiGHS.Optimizer);

julia> @variable(m, x >= 0, usym("m/s"))
x [m s⁻¹]

julia> @variable(m, y >= 0, usym("km/h"))
y [km h⁻¹]

julia> @variable(m, z, Bin, Dimensions())
z []

julia> max_speed = 4u"km/s"
4000.0 m s⁻¹

julia> @constraint(m, x + y  <= max_speed * z)
x + 0.2777777777777778 y - 4000 z <= 0 [1.0 m s⁻¹]

julia> obj = @objective(m, Max, x + 3y)
0.8333333333333334 y + x [m s⁻¹]

julia> optimize!(m)

julia> println("obj = $(value(obj))")
obj = 12000.0 [m s⁻¹]

julia> println("x = $(value(x))")
x = 0.0 [m s⁻¹]

julia> println("y = $(value(y))")
y = 14400.0 [km h⁻¹]

I think it is not ideal for the user that they need to differentiate between @variable(m, x, u"km/h") and @variable(m, x, usym("km/h")). Should both be allowed, or only the second.

I am off for vacation in a few days, but I will follow up on this. I'll join the JuMPDev in Boston and will try to bring it up as a discussion topic.

[MilesCranmer]

Awesome, I am glad it works!

To give you another update, I just refactored UnitDimensions to use a sparse vector for storage. This means that it is significantly more memory efficient to work with, and startup time is drastically reduced.

julia> q = sym_uparse("5.2*nm * kW * h")
5.2 nm kW h

julia> @btime q^2
  85.065 ns (5 allocations: 240 bytes)
27.040000000000003 nm² kW² h²

So it is actually feasible to perform calculations directly with UnitDimensions. However, note there is no automatic conversion available. So trying to do usym("nm") + usym("m") would fail as it thinks they are different dimensions. (But this would of course work if you were using u"nm" + u"m", as it converts to base units first.)

Regarding

I think it is not ideal for the user that they need to differentiate between @variable(m, x, u"km/h") and @variable(m, x, usym("km/h")). Should both be allowed, or only the second.

If I were a user, I think I wouldn't see this choice as much of an issue. One automatically converts to SI and the other leaves it, it seems pretty reasonable to have both options.

For consistency, I just exported a@us_str macro, and now export sym_uparse. To convert, I also export a function expand_units which converts to Quantity{<:Any,Dimensions} automatically. So you can do, e.g.,

julia> using DynamicQuantities

julia> q = 125us"kW*h/day"
125.0 kW h day⁻¹

julia> q2 = q^2 / 5.2
3004.8076923076924 kW² h² day⁻²

julia> expand_units(q2)
5.216680021367522e6 m⁴ kg² s⁻⁶