`Number`, `Real`, `AbstractFloat`, subtyping, rules, semantics

[goretkin]

(paraphrased and excerpted from a discussion on discourse)

Over the domain x::AbstractFloat, the rule that rewrites x - x to zero(x) is wrong, since the domain includes NaNs and Infs. However:

julia> y = SymbolicUtils.Sym{AbstractFloat}(:y)
y

julia> y - y
0

Over the "real numbers" (or any field), the rule is correct, and that is where that rule is coming from.

It strikes me as un-Julian design to have that rule defined on Sym{Real} (or Sym{Number}), and then have a specialization that turns it off for AbstractFloat. Specializing behavior is essential in the Julia ecosystem for generic programming, but specializations are not, for lack of a better word, done "willy nilly". They are supposed to be consistent with a generic definition, but perhaps more efficient or producing a more compact (symbolic!) representation. For example:

julia> typeof(1:4)
UnitRange{Int64}

julia> typeof(1:4) <: AbstractArray
true

julia> (1:4) .+ 5
6:9

julia> collect(1:4) .+ 5
4-element Vector{Int64}:
 6
 7
 8
 9

julia> collect((1:4) .+ 5) == collect(1:4) .+ 5  # consistency
true

If broadcasted addition on UnitRange{Int64} were defined in such a way that the “# consistency” property did not hold, that would be bad for generic programming!

Does that badness apply in this CAS setting? What consistency should exist across rules that are defined on types that have a subtype relationship among them? I am not sure. That’s what the question is.

My instinct is that this CAS should introduce a new type, e.g. Symbolics.Real, not Base.Real (though perhaps Symbolics.Real <: Base.Real). Defining that simplification rule on Base.Real seems wrong since I can give you an x::Base.Real for which the rule is wrong. But I should not be able to give you an y::Symbolics.Real for which the rule is wrong.

It's already not the case that Sym{Real}<:Real. What is the benefit of re-using this abstract type?

[ChrisRackauckas]

I actually think this is more Symbolics.jl related. SymbolicUtils.jl doesn't care what you do. The hardcoding to "this thing is an X which has these rules" is where Symbolics.jl comes into play.

[goretkin]

I only noticed you asked me to open it on SymbolicUtils.jl when I went to follow-up on the the discourse thread. Go ahead and close here / open there / etc as you wish.

[TymonKilich]

I think it's SymbolicUtils.jl issue since automatic simplifications happens at creation for <: Sym{Number} type(s).

julia> using SymbolicUtils

julia> @syms x::Number
(x,)

julia> x - x
0

It also autoreorders expressions, which is sometimes confusing:

julia> sin(x)^2 + cos(x)^2
cos(x)^2 + sin(x)^2

julia> r = @rule sin(~x)^2 + cos(~x)^2 => 1
sin(~x) ^ 2 + cos(~x) ^ 2 => 1

julia> r(sin(x)^2 + cos(x)^2) === nothing
true

I also agree that this may lead to some surprising behaviour

julia> x/Inf
0.0

julia> Inf/x
Inf(x^-1)

julia> x*(1/0 - 1/0)
NaNx

julia> (1/x)*x
1

julia> (1/x - 1/x)
0

julia> (1/x - 1/x)*x
0

EDIT: Maybe the solution would be to move auto simplification and reordering to some function (cannonicalize()?) but not applying it automatically at creation of <: Sym{Number} expression?

[goretkin]

I think even if simplification rules are not applied automatically/immediately (which seems like a good idea!), I think there's still a question of what good does using the type Number (e.g. in Sym{Number}) do? Is it to allow tracing through some function defined like f(x::Number) = ...?

[ChrisRackauckas]

The data structure based simplifications are something around a 10,000x to 100,000x performance gain, so it's good to do them whenever you can. I think the default numbers of a symbolic modeling language are not floating point numbers: floating point like behavior should be opt-in not opt-out since there's relatively little you can do under those semantics (and LLVM already does most of those simplifications anyways). To do anything that requires commutativity (and associativity in many cases), you require it not being a floating point number. You can only implement things like Herbie for improving floating point precision, if it doesn't have to keep floating point exact. So you really want something that's a "real number" and not a "floating point number" by default: no CAS has ever made a floating point number.

Now saying that means it's a Real is a little bit odd, I agree. Maybe making it a SymbolicReal <: Real could be a bit more concrete and then it can get more precise rules? But I think you'd want to add many rules onto Real and Number, and explicitly make AbstractFloat opt-out, because by default SymbolicReal, SymbolicComplex, etc. will all want these rules and really AbstractFloat is the odd one out.

Maybe the solution would be to move auto simplification and reordering to some function (cannonicalize()?) but not applying it automatically at creation of <: Sym{Number} expression?

That is impossible because it's part of the data structure itself, i.e. there is no Term(+,[x,y]) but an Add(Set([x,y])) which is an unordered collection and allows a lot of decreases in computational complexity. But that's all done by types, so if the dispatch doesn't do +(x,y) = Add(Set([x,y])) on the types of x and y, it will make a Term. This is how the constructor-based dispatches can be disabled on certain types, though certain things will need to make sure it knows the alternative canonical form.

[shashi]

It strikes me as un-Julian design to have that rule defined on Sym{Real} (or Sym{Number}), and then have a specialization that turns it off for AbstractFloat.

I think that is not very un-Julian, except in this case it would be annoying to end up with a 0 for Real but not for AbstractFloat...

I think there's still a question of what good does using the type Number (e.g. in Sym{Number}) do? Is it to allow tracing through some function defined like f(x::Number) = ...?

That's the idea. In that way, the type parameter has 2 roles -- one as R, the other as the supposed Julia type (arguably they should be separate.). The tracing works by first wrapping a symbolic expression using the Symbolics.wrap function, which results in a wrapper that is a subtype of Real. But we could make wrap(Symbolic{Symbolics.R}) return that same wrapper type as well. Arguably that would be more correct.

I'm open to doing this, but I doubt I will be able to get to it soon myself. So anyone who wants to actually work on such a library of domains, I would appreciate them to open a PR. Maybe we may need some changes to how the wrapper interface works (it's here https://github.com/JuliaSymbolics/Symbolics.jl/blob/master/src/wrapper-types.jl)

[goretkin]

I think that is not very un-Julian, except in this case it would be annoying to end up with a 0 for Real but not for AbstractFloat...

I would be interested in examples where removing a method definition produces the wrong answer instead of causing an error. Removing a more specific definition when a more general definition exists will result in the general definition being used, and I expect that general definition to be correct, or possibly error, but not for it to produce an incorrect answer.