Add `binomial` and `factorial` as potential dimensionless functions

[TorkelE]

Fixes https://github.com/SciML/ModelingToolkit.jl/issues/2554 https://github.com/SymbolicML/DynamicQuantities.jl/issues/124

[github-actions[bot]]

Benchmark Results

	main	bf03a593892734...	main/bf03a593892734...
Quantity/creation/Quantity(x)	3.41 ± 0.01 ns	3.71 ± 0.92 ns	0.919
Quantity/creation/Quantity(x, length=y)	3.11 ± 0.001 ns	3.42 ± 0.01 ns	0.909
Quantity/with_numbers/*real	3.11 ± 0.01 ns	3.11 ± 0.01 ns	1
Quantity/with_numbers/^int	8.05 ± 2.2 ns	8.37 ± 2.2 ns	0.963
Quantity/with_numbers/^int * real	8.67 ± 2.2 ns	8.05 ± 1.8 ns	1.08
Quantity/with_quantity/+y	4.04 ± 0.001 ns	4.04 ± 0.001 ns	1
Quantity/with_quantity//y	3.42 ± 0.01 ns	3.11 ± 0.001 ns	1.1
Quantity/with_self/dimension	3.1 ± 0.01 ns	2.79 ± 0.009 ns	1.11
Quantity/with_self/inv	3.11 ± 0.001 ns	3.11 ± 0.001 ns	1
Quantity/with_self/ustrip	2.79 ± 0.01 ns	2.79 ± 0 ns	1
QuantityArray/broadcasting/multi_array_of_quantities	0.143 ± 0.00076 ms	0.147 ± 0.0015 ms	0.972
QuantityArray/broadcasting/multi_normal_array	0.0529 ± 0.00019 ms	0.0498 ± 0.00019 ms	1.06
QuantityArray/broadcasting/multi_quantity_array	0.155 ± 0.00065 ms	0.155 ± 0.0007 ms	1
QuantityArray/broadcasting/x^2_array_of_quantities	24.1 ± 2.3 μs	24.5 ± 1.8 μs	0.982
QuantityArray/broadcasting/x^2_normal_array	4.84 ± 0.97 μs	5.14 ± 1.2 μs	0.941
QuantityArray/broadcasting/x^2_quantity_array	7.05 ± 0.28 μs	6.94 ± 0.24 μs	1.02
QuantityArray/broadcasting/x^4_array_of_quantities	0.0787 ± 0.00058 ms	0.0786 ± 0.00048 ms	1
QuantityArray/broadcasting/x^4_normal_array	0.0467 ± 0.00018 ms	0.0498 ± 0.00019 ms	0.94
QuantityArray/broadcasting/x^4_quantity_array	0.0501 ± 0.00019 ms	0.0499 ± 0.00019 ms	1
time_to_load	0.128 ± 0.0014 s	0.129 ± 0.0011 s	0.994

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).

[MilesCranmer]

I think only :factorial goes there. But :binomial needs to be elsewhere as it has two arguments rather than 1.

Also can you add it to the unittests please?

[TorkelE]

Thanks.

I've updated. I removed binomial from that set of functions, not sure where to put it back though.

[MilesCranmer]

For binomial I would just write the function definition explicitly.

[MilesCranmer]

Like

function Base.binomial(q::UnionAbstractQuantity, n::Integer)
    iszero(dimension(q)) || throw(DimensionError(q))
    return binomial(ustrip(q), n)
end

Or maybe you need to allow quantities for n as well(?) in which case you would add more methods for that (also checking for dimensionality, and stripping).

[TorkelE]

Currently, these are my tests:

    # Tests factorial (single integer input).
    for Q in (Quantity{Int64}(2u"1")), D in (Dimensions, SymbolicDimensions)
        for x in rand(3, 1:10)
            qx_dimensionless = Q(x, D)
            qx_dimensions = convert(with_type_parameters(Q, Float64, D), Q(x, dimension(u"m/s")))
            @eval @test $f($qx_dimensionless) == $f($x)
            @eval @test_throws DimensionError $f($qx_dimensions)
        end
    end

    # Tests binomial (two integer inputs).
    for Q in (Quantity{Int64}(2u"1")), D in (Dimensions, SymbolicDimensions)
        for x in rand(3, 1:10), y in rand(3, 1:10)
            qx_dimensionless = Q(x, D)
            qx_dimensions = convert(with_type_parameters(Q, Float64, D), Q(x, dimension(u"m/s")))
            qy_dimensionless = Q(y, D)
            qy_dimensions = convert(with_type_parameters(Q, Float64, D), Q(y, dimension(u"m/s")))
            @eval @test $f($y, $qx_dimensionless) == $f($y, $x)
            @eval @test $f($qy_dimensionless, $x) == $f($y, $x)
            @eval @test $f($qy_dimensionless, $qx_dimensionless) == $f($y, $x)
            @eval @test $f($qy_dimensions, $qx_dimensions) == $f($y, $x)
            @eval @test_throws DimensionError $f($qy_dimensions, $x)
            @eval @test_throws DimensionError $f($y, $qx_dimensions)
        end
    end

and this is my binomial definition

# Explicit declaration of `binomial` function.
function Base.binomial(q1::UnionAbstractQuantity, q2::UnionAbstractQuantity)
    iszero(dimension(q1)) || throw(DimensionError(q1))
    iszero(dimension(q2)) || throw(DimensionError(q2))
    return binomial(ustrip(q), ustrip(q))
end

[MilesCranmer]

Tests are failing: https://github.com/SymbolicML/DynamicQuantities.jl/actions/runs/8345114243/job/22841107243?pr=125

You can check the tests pass locally with:

julia --project=/path/to/repository -e 'using Pkg; Pkg.test()'

[TorkelE]

I am trying to get the tests to work, but having an issue here:

    # Tests factorial (single integer input).
    for Q in (Quantity{Int64},), D in (Dimensions, SymbolicDimensions)
        for x in rand(1:10, 3)
            qx_dimensionless = Q(x, D)
            qx_dimensions = convert(with_type_parameters(Q, Int64, D), Q(x, dimension(u"m/s")))
            @eval @test $factorial($qx_dimensionless) == $ffactorial($x)
            @eval @test_throws DimensionError $factorial($qx_dimensions)
        end
    end

Specifically, this is the minimal error:

julia> Quantity{Int64}(5, dimension(u"m/s"))
ERROR: MethodError: no method matching (Quantity{Int64})(::Int64, ::Dimensions{DynamicQuantities.FixedRational{Int32, 25200}})

Closest candidates are:
  (::Type{Q})(::T; kws...) where {T<:Number, T2, Q<:(AbstractQuantity{T2})}
   @ DynamicQuantities ~/.julia/dev/DynamicQuantities/src/types.jl:219
  (::Type{Q})(::Number; kws...) where Q<:AbstractQuantity
   @ DynamicQuantities ~/.julia/dev/DynamicQuantities/src/types.jl:220
  (::Type{Q})(::T, ::Type{D}; kws...) where {D<:AbstractDimensions, T<:Number, T2, Q<:(AbstractQuantity{T2})}
   @ DynamicQuantities ~/.julia/dev/DynamicQuantities/src/types.jl:217
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[35]:1

which comes from this line:

qx_dimensions = convert(with_type_parameters(Q, Int64, D), Q(x, dimension(u"m/s")))

[MilesCranmer]

I think that constructor just doesn't exist. Maybe do Q(x*u"m/s") instead?

[TorkelE]

That worked, thanks a lot! The tests passed locally now, I will push them and see what happens.

[TorkelE]

Ok, this seems good now?

[TorkelE]

good catch, sorry for missing that one!

[TorkelE]

@MilesCranmer Does this look good for starting CI again? I looked through the changes one more time and think it looks good. I also rerun CI locally and confirmed it works there.

[MilesCranmer]

Thanks!

[TorkelE]

Thanks a lot for all the help with this :)