A common "thick numbers" package?

[timholy]

CC @giordano, @devmotion (and others?)

I've not been Watching this package but it looks like interesting and important changes are afoot, and I wonder if it's an opportunity for some shared development among packages like Measurements and possibly even dual-number packages like ForwardDiff (and probably more, feel free to CC others). All of these packages share a somewhat-awkward desire to do math on objects that act like numbers but also have some "thickness" to them. (In ForwardDiff the thickness is infinitesimal, so maybe it's not quite as relevant to this discussion.) All of these packages have been tempted to inherit from Real, but many have recognized that this has some issues. I wonder if it's time to create a core package---let me throw out ThickNumbers.jl as a possible name---that defines abstract types and exports common operators, traits, fallbacks, errors, etc., but leaves all concrete implementations to downstream packages.

To make this a bit concrete, here's an example from the development (master) branch of IntervalArithmetic:

julia> using IntervalArithmetic, IntervalArithmetic.Symbols

julia> 1..3 == 1..3
ERROR: ArgumentError: == is purposely not supported, use ≛ instead
Stacktrace:
 [1] ==(::Interval{Float64}, ::Interval{Float64})
   @ IntervalArithmetic ~/.julia/dev/IntervalArithmetic/src/intervals/interval_operations/boolean.jl:31
 [2] top-level scope
   @ REPL[2]:1

julia> 1..3 ≛ 1..3
true

julia> 1..3 ≛ 1..3.1
false

Note the use of a new operator ≛. My understanding is that IntervalArithmetic is trying to adhere more strictly to the fundamental principle that f(X, Y) (here f is ==) includes all possible values of f(x, y) for any x ∈ X, y ∈ Y; since almost any independent choice of x, y are not equal, testing interval equality is nearly useless. One could make a case for Interval{Bool}(false, true), but having the operation be undefined is at least safe, and gets you to the same result (an error) that putting a == b inside an if statement would generate even if it did return Interval{Bool}(false, true) (because if requires a Bool). This is just one example among many.

I tend to think of Measurement as kind of like an interval, but less strict and less affected by the dependency problem. I actually have a practical use case where I'd like to be able to support users choosing whether they want to perform a computation using IntervalArithmetic (guaranteed but probably overly-conservative bounds) or Measurements (approximate but probably tighter bounds) and hence it would be interesting to consider making the two packages more compatible.

While I haven't actually written a line of code, I wonder if a shared "abstract package" might standardize some of the difficult decisions we have in common, or at least provide a social gathering point for discussing such issues. Julia's own base/numbers.jl and base/operators.jl might serve as inspiration for what could go in there, as well as noticing and sharing relevant patterns that are currently in the packages that might eventually use ThickNumbers.jl.

[OlivierHnt]

Thanks for the opening this issue!
It would be nice indeed to have a common interface. But I am not sure on the design either; to do this we would need to write down what are the expected features of "thick numbers" and hopefully see if we can reach a consensus.

To add context: interval arithmetic has stringent specifications which play poorly with generic code intended for real numbers (in the mathematical sense). To illustrate, let me copy-paste @Kolaru's recent example

f(x, y) =  (x == y) ? 1 : exp(x - y)

which fails for the reasons mentioned already in @timholy's comment.

So to prevent silent failures, we have to fight quite a bit against Julia's flexibility (that's the reason for the recent changes on master) and this is why having Interval <: Real is not great. On master, we still have Interval <: Real but we removed all comparisons ==, <, etc. and conversion/promotion with Number; in PR #572 we are also discussing if intersect, union, etc. should follow the Base definition or just error to disambiguates with their interval arithmetic meaning.

It may also be worth mentioning here that recently we opened #584 as a way to recover interoperability with the ecosystem by exploiting the decoration system. This came up when thinking about conversion, see recently #580. One scenario could be to define

struct Interval{T} <: Real
    interval   :: BareInterval{T}
    decoration :: Decoration 
end

struct BareInterval{T}
    lo :: T
    hi :: T
end

Then, Interval can be made more permissive since a decoration could flag "unsafe operations" (converting a Number to an Interval for instance). On the other hand, BareInterval will be very conservative.

[Kolaru]

To compare with other similar packages having the same problem (Unitful.jl, Measurements.jl and ForwardDiff.jl for what I know), I think it may be useful to explicit what we want or not from the Real supertype. That way we may get a clearer view of what we have in common.

What we want:

• Arithmetic.
• Mathematical functions.
• Interoperability with the ecosystem ^^'.

What we do not want:

• Automatic promotion to or from float or other Real types. Correctness can simply not be guaranteed in general.
• Number as 1 element 0-size arrays (broadcastable, valid for things like iteration, and set operation like intersect). We don't want that because it is ambiguous since intervals are sets by their own rights, which conflict with Base definition of e.g. intersect for a Real. Being broadcastable alone would be fine, maybe.

What we can kind of live with:

• Ordering (< and such) and comparison (==). We can in principle use ternary logic for that, that would be guaranteed, but it is still slightly ambiguous (especially for ==).
• Mixing different Real types. Float64 * Interval has unclear interpretation, because we don't know if the Float64 is an exact trustable number. As @OlivierHnt mentioned we could tackle that by extending the concept of DecoratedInterval, and tag computation that may have lost validity that way.

I am pretty sure that I am missing some elements, but that's what come to mind right now.

Does it leaves enough in common with other packages ? I am not sure.

[timholy]

Update (sorry for the delay, I've been immersed in other things): I've realized that Measurements.jl doesn't implement the arithmetic rules I want:

julia> using Measurements

julia> (0 ± 1) * (0 ± 1)
0.0 ± 0.0

Instead I want something consistent with the algebra of random variables (what I think of as "soft intervals"). AFAICT there isn't really a package doing what I want, so I may start one (locally for now), and as I develop it I'll look for common interests. I may have a more informed view on this in a couple of weeks.

[giordano]

CC @baggepinnen for MonteCarloMeasurements.jl .

I don't think it's particularly easy to define a common generic interface. For example I feel quite strongly about making Measurements.Measurement <: AbstractFloat, which I don't think would make sense for IntervalArithmetic.

@timholy Measurements.jl implements linear error propagation, which is a simple and commonly used method in Physics, but it relies on the first derivative of the function being the "dominant" one, hypothesis which can locally fail for simple functions like sin/cos (think of their stationary points). One could extend the propagation to second-order derivatives to work around this limitation, another approach is Monte Carlo sampling, which should be equivalent to linear error propagation when its condition is satisfied (albeit much slower), and provide more "meaningful" results otherwise.

[timholy]

There's an initial implementation at https://github.com/timholy/ThickNumbers.jl. Durng this preliminary phase, I'll host it in my personal GitHub account, but if we decide to adopt it I'll transfer it to JuliaMath.

It's preliminary in the sense that I've not actually tried it in practical usage, and I don't plan to register it until I've created or rebased at least two packages on it, so we know if the generalizations "work" in real-world practice.

However, it's got a lot going for it already:

• extensive documentation (still waiting for it to be deployed at the time of writing, so I haven't checked availability)
• default implementations for a subset of IEEE Std 1788-2015 (more presumably to come)
• a subdir package, ThickNumbersInterfaceTests, which can be used to determine whether a concrete subtype satisfies the interface.
• support for ForwardDiff across all ThickNumber subtypes, even though they are <: Number rather than <: Real.

Note that it mirrors the "new" direction for IntervalArithmetic, but with a couple changes in terms of unicode operators. I decided to aim for "dot" consistency, riffing off of ⪽. So most of the operators I chose have a "dot" in them. I'm happy to discuss and modify these choices.

[timholy]

There were some initial issues with doc depoloyment, but those have been fixed and additional improvements made.

[timholy]

Toy package serving as a proof-of-principle: https://github.com/HolyLab/GaussianRandomVariables.jl

[lbenet]

CC To @AnderGray, who has experience and may add something to the discussion.

[AnderGray]

I think a package implementing a base uncertain number type, with common interface, would be very useful. I would have used it when making ProbabilityBoundsAnalysis.jl and MomentArithmetic.jl.

Note: if you're looking for an easy-ish way to reduce puffiness in interval arithmetic, centred from interval arithmetic is still rigorous but uses the first derivative (we can use ForwardDiff.jl) to reduce puffiness. I think when you use higher orders, it generalises to Taylor Models.

Unfortunately gaussian random variables aren't closed under arithmetic operations (except linear transformations), so the distribution's shape will change for anything non-linear. Essentially why first-order error propagation is only accurate for linear-enough functions.

You can however perform computations with moments, or at least bound them in non-linear transformations. So unless you additionally track the shapes, the transformed moments become intervals. Afterwards however, probabilities can be bounded from the moment, using e.g. the Chebyshev inequality.

If useful, we did publish open access article on the topic, and a sort of repository/package. Inside the paper you can find some tables containing the arithmetic rules we worked out (similar to the algebra of random variables)

Reducing puffiness is all about correlations:

julia> using MomentArithmetic, IntervalArithmetic

julia> a = Moments(0, 1, -1 .. 1)  # mean 0, var 1, range [-1, 1]
julia> b = Moments(0, 1, -1 .. 1)

julia> a - b   # Default is unknown correlation
moment: 	  ~ ( mean  = 0, var = [0.0,4.0] , range = [-2.0,2.0] )

julia> subIndep(a, b)  # Independence gives precise variance
moment: 	  ~ ( mean  = 0, var = 2 , range = [-2.0,2.0] )

At the end of the paper there is some discussion on using covariances to reduce puffiness in arithmetic like in the above.

[timholy]

Thanks, this looks useful. In case you're not aware, I tried building MomentArithmetic on 1.10 and got this:

julia> using MomentArithmetic
Precompiling MomentArithmetic
        Info Given MomentArithmetic was explicitly requested, output will be shown live
WARNING: could not import ProbabilityBoundsAnalysis.plot into MomentArithmetic
WARNING: Method definition (::Type{MomentArithmetic.Moments})() in module MomentArithmetic at /home/tim/.julia/packages/MomentArithmetic/x2Mva/src/Moments.jl:43 overwritten at /home/tim/.julia/packages/MomentArithmetic/x2Mva/src/Moments.jl:52.
ERROR: LoadError: Method overwriting is not permitted during Module precompile.
Stacktrace:
 [1] top-level scope
   @ ~/.julia/packages/MomentArithmetic/x2Mva/src/Moments.jl:52
...

1.10 is pickier about turning things that may have been warnings on previous releases into errors, because it was found that in some cases unfixed warnings were the source of serious, very difficult-to-track-down bugs. (In other cases they are fairly harmless, and this may be such a case.)

[timholy]

I also like this line:

There are no infinitely massive body weights (despite recent trends in western dietary health).

😆