Violation of the (informal) `Number` interface

[baggepinnen]

I have been hitting my head at #330 for quite some time now, and I just realized that it's due to the type TaylorN claiming to be a <: Number, but deviates from the interface for getindex for numbers:

julia> x, = set_variables("x")
1-element Vector{TaylorN{Float64}}:
  1.0 x + 𝒪(‖x‖²¹)

julia> x
 1.0 x + 𝒪(‖x‖²¹)

julia> x isa Number
true

julia> x === x[]
true

julia> x === x[1]
false

For any scalar x::Number, the operations x[] and x[1] are expected to be identical, but for TaylorN they aren't.

[PerezHz]

Discussion related to #224 might be relevant here.

[baggepinnen]

Yeah, that issue did bring up this particular issue, but has been closed without the actual issue being resolved :/

From #224

just let me understand what is the benefit that such a change could bring.

The benefit of respecting the Number interface is correctness, that you can pass TaylorN into user and library defined functions, otherwise you cannot rely on the result being correct

[lbenet]

I agree with your point as well as that in #224. Let me first say that the fact that we have that TaylorN{T} is a Number is for certain convenience... no mathematical implications, just in the same sense that Real means floating point, and not $\mathbb{R}$.

Now, in your example, x[] == x shows that it behaves as a "scalar" quantity, just as y = 1.0; y == y[] does. Once agreeing with this point, it is so comfortable to refer to the first-order coefficient of a polynomial as pol[1], or to the zeroth as pol[0], or the fifth as pol[5], that we simply do that. That's all we do, with no implications beyond that.

[baggepinnen]

I understand that it is convenient, but the implication is that you cannot pass TaylorN into a library function and rely on it doing the correct thing, since the library function might index a scalar like scalar[1] and that should be okay, but it has unintended consequences for TaylorN. It may sound silly that a library function should index a scalar like that, but it's in fact rather common and the reason scalars can be indexed at all is that it simplifies writing generic code that operates on both scalars and arrays.

I was facing similar problems in MonteCarloMeasurements.jl where indexing was used to represent a sample of a quantity that behaved like a number. After having encountered silently incorrect results a number of times I concluded that this was really dangerous and had to abandon the indexing convenience. My documentation still has this remainder to this day:

[lbenet]

I truly see your point, and indeed convenience is the property challenged here. Yet, you can argue just the same but in the opposite direction, about a program which uses scalar[] or scalar[1] and on the hope that it works always, specially when it is standard/cleaner to refer to a scalar quantity as scalar. (Incidentally, does this affect somehow performance? My guess is no, but I'm wondering about type stability.)

After saying this, let me argue that the fact that TaylorN is essentially a vector of HomogenoeousPolynomials suggests such "convenient" referencing... polynomials are isomorphic to vectors by using the first component indexing.

I guess this discussion boils down to "a matter of taste" (and tastes are not discussed 😄): If I would like that a function works with vectors and scalars, I would exploit (multiple) dispatch, in the sense of writing two or more functions that differ in some types; your choice is to have a single function that works universally.

[baggepinnen]

The problem is that my taste does not matter if I call a function inside a library I have not written. To make use of TaylorN, I now have to inspect the entire call stack to verify that problematic scalar indexing isn't used anywhere, because if it is, I will get the wrong result when passing TaylorN in there.

The problem extends to iteration as well:

julia> t = TaylorN(1)
 1.0 x₁ + 𝒪(‖x‖⁷)

julia> for i in t
       println(i)
       end
 0.0
 1.0 x₁
 0.0
 0.0
 0.0
 0.0
 0.0

julia> for i in 1
       println(i)
       end
1

your choice is to have a single function that works universally.

This is not only my choice, there is an infamous PR that tried to change this behavior so that scalars were not iterable and indexable, but so much of Base julia relied on this behavior that it was abandoned.

• RFC: eliminate getindex(::Number) JuliaLang/julia#19700

It has been observed many times that this functionality is used all over the place

Ref: https://discourse.julialang.org/t/numbers-as-single-element-collections/1621/5?u=baggepinnen

and any time a TaylorN is passed into such code, the user gets the wrong result without any error

[lbenet]

From JuliaLang/julia#19700, my take-away is to avoid indexing scalars. We do index though, with the meaning it is described, and this is explained in the documentation.

Again, this is a matter of taste, in this case, our tastes are distinct. Sorry for the inconveniences.