Additional utilities for identity functions, `+`, and `-` #23

MilesCranmer · 2023-06-12T03:51:47Z

This fixes the behavior of one, zero, oneunit as pointed out by @mcabbott in #6.

This also extends + and - for ::Quantity to ::Number, given the quantity is dimensionless.

Could you take a look @mcabbott?

github-actions · 2023-06-12T03:57:35Z

Benchmark Results

	main	`0257746`...	t[main]/t[`0257746`...]
creation/Quantity(x)	2.8 ± 0.1 ns	3.2 ± 0.1 ns	0.875
creation/Quantity(x, length=y)	3.3 ± 0.1 ns	3.2 ± 0.1 ns	1.03
time_to_load	0.0842 ± 0.001 s	0.0989 ± 0.00073 s	0.851
with_numbers/*real	3.2 ± 0.1 ns	2.9 ± 0.1 ns	1.1
with_numbers/^int	0.034 ± 0.0027 μs	0.034 ± 0.0031 μs	1
with_numbers/^int * real	0.0348 ± 0.0021 μs	0.0334 ± 0.0024 μs	1.04
with_quantity/+y	5.1 ± 0 ns	5.8 ± 0.2 ns	0.879
with_quantity//y	1.6 ± 0.1 ns	3.3 ± 0.1 ns	0.485
with_self/dimension	1.6 ± 0.1 ns	1.6 ± 0.1 ns	1
with_self/inv	1.6 ± 0.1 ns	3.3 ± 0.1 ns	0.485
with_self/ustrip	1.6 ± 0.1 ns	1.6 ± 0.1 ns	1

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

src/utils.jl

j-fu · 2023-06-12T10:19:48Z

Essentially, this PR seems to hit in the same direction as #18, but just to remark: this does not define an universal zero(DynamicQuantity(T)), and without this, DynamicQuantity is unable to act like a "normal" number.

BTW I found an intersting paper here:

https://link.springer.com/article/10.1023/B:JOMC.0000034933.47311.93

And two other papers citing it concerned with programmable formalization:
https://scholar.google.de/scholar?cites=5148028135146171185&as_sdt=2005&sciodt=0,5&hl=de

Up to what I could understand, these don't work with "universal zero" as well.

MilesCranmer · 2023-06-12T10:35:58Z

Thanks, I see. Do you want to take a look at this too?

I don't know what you mean by zero(DynamicQuantity(T)), could you elaborate? The zero definition is given here:

DynamicQuantities.jl/src/utils.jl

Lines 58 to 61 in 38ef94a

    
           Base.zero(q::Quantity) = Quantity(zero(ustrip(q)), dimension(q)) 
        
           Base.zero(::Dimensions) = error("There is no such thing as an additive identity for a `Dimensions` object, as + is only defined for `Quantity`.") 
        
           Base.zero(::Type{<:Quantity}) = error("Cannot create an additive identity for a `Quantity` type, as the dimensions are unknown. Please use `zero(::Quantity)` instead.") 
        
           Base.zero(::Type{<:Dimensions}) = error("There is no such thing as an additive identity for a `Dimensions` type, as + is only defined for `Quantity`.")

As is stated, it can only be defined for ::Quantity. Other definitions lack the dimensional information (when given as a type), or are defined on a Dimension themselves, and are therefore ill-posed.

(Only in Unitful.jl is the dimension information accessible from the type itself; whereas here it is stored in the value - and therefore constrains the zero() definition to be like this.)

MilesCranmer · 2023-06-12T10:41:48Z

Btw, note that with this PR, you can now add numbers to dimensionless quantities:

julia> x = Quantity(0.5)
0.5 

julia> x + 1.0
1.5 

julia> typeof(x + 1.0)
Quantity{Float64, DynamicQuantities.FixedRational{Int32, 25200}}

src/utils.jl

mcabbott · 2023-06-12T04:17:00Z

src/utils.jl

+Base.one(d::Dimensions) = one(typeof(d))
+
+# Additive identities:
+Base.zero(q::Q) where {T,Q<:Quantity{T}} = Quantity(zero(ustrip(q)), dimension(q))


Maybe a style thing, but you could shorten many lines by deleing unused type parameters:

Suggested change

Base.zero(q::Q) where {T,Q<:Quantity{T}} = Quantity(zero(ustrip(q)), dimension(q))

Base.zero(q::Quantity) = Quantity(zero(ustrip(q)), dimension(q))

src/utils.jl

mcabbott · 2023-06-12T12:47:55Z

src/math.jl

+Base.:*(l::Quantity, r) = Quantity(l.value * r, l.dimensions)
+Base.:*(l, r::Quantity) = Quantity(l * r.value, r.dimensions)
+Base.:*(l::Dimensions, r) = Quantity(r, l)
+Base.:*(l, r::Dimensions) = Quantity(l, r)


It seems that both Dimensions and Quantity are given many number-like methods. Might it be simpler to only do this for Quantity, and use Quantity(true, Dimension(...)) in place of a naked Dimension?

Or some special One if true is too narrow. But I'm not so sure how much it's intended to support Quantity([1,2,3], time=1). It would seem simplest to me to restrict Quantity to real numbers, but maybe you have other plans.

Actually, when would methods like *(::Number, ::Dimension) be used at all? #22 defines const m = Quantity(1.0, length=1). It could instead use this method and define m::Dimension, but instead seems to focus on numbers.

It seems that both Dimensions and Quantity are given many number-like methods. Might it be simpler to only do this for Quantity, and use Quantity(true, Dimension(...)) in place of a naked Dimension?

I'm not sure if this is possible because Quantity.value and Dimensions are different number systems: .value is a regular number, while Dimensions.(field) are powers. So, for example, :*(::Dimensions, ::Dimensions) actually adds fields, rather than multiplies them.

Actually, when would methods like *(::Number, ::Dimension) be used at all? #22 defines const m = Quantity(1.0, length=1). It could instead use this method and define m::Dimension, but instead seems to focus on numbers.

This is just a convenience thing that felt more intuitive than using Quantity() directly. For example:

x = Quantity(1.0, mass=3, length=1) y = 0.5 * dimension(x)

which puts y in the units of x.

You could do Quantity(0.5, dimension(x)), but I think after #22 merges, the Quantity() constructor will be used less often, so this multiplication strategy seems useful. Wdyt?

(Sorry, just realized all my messages were listed as pending because they are submitted as a batch when I hit "Review changes"...)

j-fu · 2023-06-12T14:53:30Z

Ok, now found a way to make an MWE for what I am asking for:
Basically, this is about linear algebra with units.
This one has no mixed units yet, but highlights some interesting points.

julia> using DynamicQuantities, LinearAlgebra
julia> A=rand(10,10).*Quantity(1,length=1);
julia> b=rand(10).*Quantity(1,time=1);

In order to make A\b work , we need:

julia> Base.oneunit(::Type{Quantity{T,R}}) where {T,R} = Quantity(one(T))
julia> Base.zero(::Type{Quantity{T,R}}) where {T,R} = Quantity(zero(T))
julia> Base.adjoint(a::Quantity{T,R}) where {T,R} = Quantity(adjoint(a.value),a.dimensions)

If not for RowNonZero pivoting (>=1.9) , we also would need Base.abs stripping the dimesions, or Base.isless. RowNonZero was made for having linear algebra working with finite fields.

julia> LinearAlgebra.lupivottype(::Type{Quantity{T,R}}) where {T,R} = LinearAlgebra.RowNonZero()

So with these settings we get

julia> A\b
10-element Vector{Quantity{Float64, DynamicQuantities.FixedRational{Int32, 25200}}}:
 -0.7353181472910584 m⁻¹ s
 3.697912215371167 m⁻¹ s
 4.3285995603007335 m⁻¹ s
 5.1254898127729 m⁻¹ s
 -0.7304999810589047 m⁻¹ s
 -6.27314096327331 m⁻¹ s
 0.2909458298595246 m⁻¹ s
 -2.954423433552737 m⁻¹ s
 -2.833125030348257 m⁻¹ s
 -2.6987760225010198 m⁻¹ s

If I find time I als will try to make an example with mixed quantities in the matrix and the RHS.

Besides of defining these methods in this package, with AbstractQuantities, it would be possible to define a DirtyQuantity (possibly in a different package) which would work the way I described.

MilesCranmer · 2023-06-12T15:26:38Z

In order to make A\b work , we need:

julia> Base.oneunit(::Type{Quantity{T,R}}) where {T,R} = Quantity(one(T))
julia> Base.zero(::Type{Quantity{T,R}}) where {T,R} = Quantity(zero(T))
julia> Base.adjoint(a::Quantity{T,R}) where {T,R} = Quantity(adjoint(a.value),a.dimensions)

Base.adjoint seems like something we can definitely implement.

I'm not sure we should add Base.zero(::Type{Quantity{T,R}}) or Base.oneunit(::Type{Quantity{T,R}}). See the docstring for zero:

help?> zero
search: zero zeros iszero set_zero_subnormals get_zero_subnormals count_zeros RoundToZero leading_zeros trailing_zeros RoundFromZero finalizer

  zero(x)
  zero(::Type)


  Get the additive identity element for the type of x (x can also specify the type itself).

  See also iszero, one, oneunit, oftype.

Meaning that zero(x) needs to be able to be added to x. But zero(::Type{Quantity}) violates this as it can produce dimension errors. So I think only zero(::Quantity) makes sense to implement here.

Here is the docstring for oneunit:

help?> oneunit
search: oneunit DimensionlessQuantity

  oneunit(x::T)
  oneunit(T::Type)


  Return T(one(x)), where T is either the type of the argument or (if a type is passed) the argument. This differs from one for dimensionful quantities:
  one is dimensionless (a multiplicative identity) while oneunit is dimensionful (of the same type as x, or of type T).

so it needs to know about dimensions - hence it requires ::Quantity as input (since dimensions are stored in value), rather than ::Type{Quantity}.

we also would need Base.abs stripping the dimesions

We should be very hesitant to use a different algebra than standard unit packages unless there is a very good reason to do so. Both astropy and Unitful propagate units through abs, so DynamicQuantities should probably do so as well:

Astropy:

[ins] In [2]: abs(1.5*u.km/u.s)
Out[2]: <Quantity 1.5 km / s>

Unitful:

julia> abs(1.5u"km/s")
1.5 km s⁻¹

Base.isless

We can definitely add this. And I guess should raise DimensionError when the dimensions are different.

Besides of defining these methods in this package, with AbstractQuantities, it would be possible to define a DirtyQuantity (possibly in a different package) which would work the way I described.

See #24 for a first attempt at implementing AbstractQuantity, which might make this possible.

MilesCranmer · 2023-06-12T15:44:10Z

For the specific example, it's probably better to just do:

julia> Quantity(ustrip.(A) \ ustrip.(b), dimension(first(b)) / dimension(first(A)))

I think by overloading abs with ustrip, and oneunit with a dimensionless number, one is avoiding the dimensionful calculations regardless.

Edit: also note that this doesn't work for Unitful.jl as well:

julia> A \ b
ERROR: MethodError: no method matching (Quantity{Float64})(::Float64)

mcabbott · 2023-06-12T17:07:12Z

In order to make A\b work [...] Vector{Quantity{Float64,

This is one path, arrays of Quantity. In Unitful.jl, arrays of uniform unit are in principle free, and can be re-interpreted, like PainterQubits/Unitful.jl#437 . But here, they involve storing many Dimensions objects, and can't be re-interpreted to feed to BLAS.

The other possible path is Quantity{Vector{Float64}} etc. To wrap a dense array, and define linear algebra by un-wrapping, calling the usual (possibly BLAS) routines, and re-wrapping the result.

Edit: Quantity(ustrip.(A) \ ustrip.(b), ... is I think suggesting that methods for \(::Vector{Quantity{Float64,... could similarly un-wrap and re-wrap. But more copies are involved.

Both astropy and Unitful propagate units through abs,

This seems essential, abs(x) surely is still in the same units as x, and will show up in user code not just (as here, I think) in the internals of \.

j-fu · 2023-06-12T20:50:33Z

Well, arrays of quantity are certainly interesting for heterogeneous physical systems. I am interested in "multiphysics" PDE systems in this context.

With the switch to AbstractQuantity it is quite straightforward to define a type with working semantics for linear algebra.

With

Base.zero(::Type{XQuantity{T,R}}) where {T,R} = XQuantity(zero(T))
Base.zero(::XQuantity{T,R}) where {T,R} = XQuantity(zero(T))
Base.one(::Type{XQuantity{T,R}}) where {T,R} = XQuantity(one(T))
Base.one(::XQuantity{T,R}) where {T,R} = XQuantity(one(T))
function Base.:+(l::XQuantity, r::XQuantity)
    if iszero(l.value)
        r
    elseif iszero(r.value)
        l
    else
        dimension(l) == dimension(r) ? new_quantity(typeof(l), ustrip(l) + ustrip(r), dimension(l)) : throw(DimensionError(l, r))
    end
end

adjoint(q::XQuantity{T,R}) where {T,R} = XQuantity(adjoint(q.value),q.dimensions)
Base.isless(q,r) = q.dimensions==r.dimensions && q.value<r.value

this seems to get me through dense A\b and through sparspaklu(sparse(A))\b .

Now adjoint and isless seem to be undisputed. I agree that both zero and one are mathematially inconsistent, and that + is even worse. OTOH having pluggable arithmetics of this kind is really intriguing to me.

src/utils.jl

MilesCranmer · 2023-06-12T14:44:51Z

src/math.jl

+Base.:*(l::Quantity, r) = Quantity(l.value * r, l.dimensions)
+Base.:*(l, r::Quantity) = Quantity(l * r.value, r.dimensions)
+Base.:*(l::Dimensions, r) = Quantity(r, l)
+Base.:*(l, r::Dimensions) = Quantity(l, r)


It seems that both Dimensions and Quantity are given many number-like methods. Might it be simpler to only do this for Quantity, and use Quantity(true, Dimension(...)) in place of a naked Dimension?

I'm not sure if this is possible because Quantity.value and Dimensions are different number systems: .value is a regular number, while Dimensions.(field) are powers. So, for example, :*(::Dimensions, ::Dimensions) actually adds fields, rather than multiplies them.

MilesCranmer · 2023-06-12T14:48:24Z

src/math.jl

+Base.:*(l::Quantity, r) = Quantity(l.value * r, l.dimensions)
+Base.:*(l, r::Quantity) = Quantity(l * r.value, r.dimensions)
+Base.:*(l::Dimensions, r) = Quantity(r, l)
+Base.:*(l, r::Dimensions) = Quantity(l, r)


Actually, when would methods like *(::Number, ::Dimension) be used at all? #22 defines const m = Quantity(1.0, length=1). It could instead use this method and define m::Dimension, but instead seems to focus on numbers.

This is just a convenience thing that felt more intuitive than using Quantity() directly. For example:

x = Quantity(1.0, mass=3, length=1) y = 0.5 * dimension(x)

which puts y in the units of x.

You could do Quantity(0.5, dimension(x)), but I think after #22 merges, the Quantity() constructor will be used less often, so this multiplication strategy seems useful. Wdyt?

src/utils.jl

MilesCranmer · 2023-06-14T18:19:21Z

src/math.jl

+Base.:*(l::Quantity, r) = Quantity(l.value * r, l.dimensions)
+Base.:*(l, r::Quantity) = Quantity(l * r.value, r.dimensions)
+Base.:*(l::Dimensions, r) = Quantity(r, l)
+Base.:*(l, r::Dimensions) = Quantity(l, r)


(Sorry, just realized all my messages were listed as pending because they are submitted as a batch when I hit "Review changes"...)

MilesCranmer · 2023-06-14T18:20:35Z

Going to merge this now so we can have the fixed one/zero/oneunit included in 0.4.0, but we can refactor these in the future if needed.

MilesCranmer commented Jun 12, 2023

View reviewed changes

src/utils.jl Outdated Show resolved Hide resolved

MilesCranmer added 5 commits June 12, 2023 07:10

Add proper one, zero, and oneunit calls

957ec74

Extend + and - to regular numbers

ca3f6f9

Test new behavior of zero, one, oneunit, +, -

060d147

Remove ::Number restrictions on math

5bcce75

Clean up zero and oneunit

0257746

MilesCranmer force-pushed the additional-utils branch from 38ef94a to 0257746 Compare June 12, 2023 11:10

MilesCranmer changed the base branch from main to units June 12, 2023 11:10

mcabbott reviewed Jun 12, 2023

View reviewed changes

MilesCranmer added 2 commits June 12, 2023 10:34

Fix behavior of one(::Quantity)

3300a00

Ensure dimensionless when using one(::Quantity)

6538750

MilesCranmer added 2 commits June 12, 2023 17:58

Add isless and == for numbers

bada5cf

Merge branch 'units' into additional-utils

2412d55

Base automatically changed from units to main June 14, 2023 18:16

MilesCranmer commented Jun 14, 2023

View reviewed changes

MilesCranmer merged commit 2044168 into main Jun 14, 2023

MilesCranmer deleted the additional-utils branch June 14, 2023 18:20

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Additional utilities for identity functions, `+`, and `-` #23

Additional utilities for identity functions, `+`, and `-` #23

MilesCranmer commented Jun 12, 2023

github-actions bot commented Jun 12, 2023 •

edited

Loading

j-fu commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023 •

edited

Loading

MilesCranmer commented Jun 12, 2023 •

edited

Loading

mcabbott Jun 12, 2023

mcabbott Jun 12, 2023

MilesCranmer Jun 12, 2023

MilesCranmer Jun 12, 2023

MilesCranmer Jun 14, 2023

j-fu commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023 •

edited

Loading

mcabbott commented Jun 12, 2023 •

edited

Loading

j-fu commented Jun 12, 2023 •

edited

Loading

MilesCranmer Jun 12, 2023

MilesCranmer Jun 12, 2023

MilesCranmer Jun 14, 2023

MilesCranmer commented Jun 14, 2023

	Base.zero(q::Q) where {T,Q<:Quantity{T}} = Quantity(zero(ustrip(q)), dimension(q))
	Base.zero(q::Quantity) = Quantity(zero(ustrip(q)), dimension(q))

Additional utilities for identity functions, +, and - #23

Additional utilities for identity functions, +, and - #23

Conversation

MilesCranmer commented Jun 12, 2023

github-actions bot commented Jun 12, 2023 • edited Loading

Benchmark Results

Benchmark Plots

j-fu commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023 • edited Loading

MilesCranmer commented Jun 12, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

j-fu commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023

MilesCranmer commented Jun 12, 2023 • edited Loading

mcabbott commented Jun 12, 2023 • edited Loading

j-fu commented Jun 12, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MilesCranmer commented Jun 14, 2023

Additional utilities for identity functions, `+`, and `-` #23

Additional utilities for identity functions, `+`, and `-` #23

github-actions bot commented Jun 12, 2023 •

edited

Loading

MilesCranmer commented Jun 12, 2023 •

edited

Loading

MilesCranmer commented Jun 12, 2023 •

edited

Loading

MilesCranmer commented Jun 12, 2023 •

edited

Loading

mcabbott commented Jun 12, 2023 •

edited

Loading

j-fu commented Jun 12, 2023 •

edited

Loading