Add arithmetic & related operations when unambiguous

[aplavin]

I'd like to propose adding basic operations on intervals, such as multiplication by another value (the specific list of operations to be defined).
Multiplication is supported by different Julia objects, not just numbers. For example, arrays – and this is often convenient:

2 * [1,2,3]

using Unitful
[1,2,3]u"m"

I don't see any discussions of such a proposal before. What do you think?

[hyrodium]

I think multiplication such as (1..2) * 5 should throw an error because such an operation is not a common notation in mathematics. This error is consistent with *(::Set, ::Number).

julia> Set([1,2,3])*3
ERROR: MethodError: no method matching *(::Set{Int64}, ::Int64)

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:587
  *(::Missing, ::Number)
   @ Base missing.jl:123
  *(::BigFloat, ::Union{Int16, Int32, Int64, Int8})
   @ Base mpfr.jl:447
  ...

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

The reason 2 * [1,2,3] does not throw an error is because it is just a scalar multiplication on a vector space.

[hyrodium]

I don't see any discussions of such a proposal before.

The discussions about broadcasting are related: #32, #55

[daanhb]

Just pointing this is implemented in DomainSets using broadcast, and it is frequently useful, e.g. for affine maps:

julia> using IntervalSets, DomainSets

julia> 2 .* (1..2) .+ 5
7.0 .. 9.0

The more general notation cos.( 0.0..2.0) also makes sense, but one cannot generically deduce what interval this corresponds to for any function. It results in a mapped domain in DomainSets.

[dlfivefifty]

There was a lot of discussion when the package was created. The problem is one probably wants correctly rounded a la IntervalArithmetic.jl so this isn't easy to add, and makes the package very heavy.

[aplavin]

@hyrodium those linked discussions also include operations with a less clear and more ambiguous meaning. For example, what (1..2) ./ (-1..1) should be – is it -Inf..Inf or what?

I propose to focus only on operations where there are (hopefully!) no doubts what they mean. Like, what else could (x::Real) * int::Interval mean other than Interval((x .* endpoints(int))...)? Or (1..2)u"m" === (1u"m")..(2u"m").

Actually, different kinds of unit conversions is my main motivation here. Would be nice to have stuff like int_km = int_m / 1000, int_deg = rad2deg(int_rad) (with just reals, no unit types), and int = int_m * u"m" (with Unitful).

[aplavin]

@daanhb interesting, thanks for pointing it out! Not great that it involves heavy piracy though...
Also, multiplication doesn't really work with Unitful it seems.

[daanhb]

@daanhb interesting, thanks for pointing it out! Not great that it involves heavy piracy though... Also, multiplication doesn't really work with Unitful it seems.

Good catch re Unitful, it seems DomainSets makes incorrect assumptions about Number types when promoting types. This will be good to fix.

I agree on the piracy comment, and there are many more instances of this. I'd really like to clean this up, but just to clarify: can we make a distinction between changing behaviour of an interval and extending behaviour? The former is a definite no-go and I hope there are no instances of that, but making something work which without loading DomainSets would result in an error seems different.

In this case, DomainSets defines broadcast for Domain. If the definition of Domain moves to DomainSetsCore, so that Intervals.jl no longer "owns" it, would that settle it?

[daanhb]

Meanwhile DomainSets is updated to fix the promotions going on:

julia> using DomainSets, Unitful

julia> (1..3) .* 3u"m"
3 m .. 9 m

julia> (2u"m"..3u"m") .* 4.0
8.0 m .. 12.0 m

julia> (2u"m"..3u"m") .* 4
8 m .. 12 m

There is no special treatment to avoid the rounding issues that were discussed in the linked issues. Somehow integer endpoints got promoted to floats when applying the map, and that is also fixed (an interval with integer endpoints times an integer remains an interval with integer endpoints).

It seems that using just * instead of .* also works, which is probably undesired. It may have been like that for some time. I'll make a separate issue about that.

[aplavin]

can we make a distinction between changing behaviour of an interval and extending behaviour? The former is a definite no-go and I hope there are no instances of that, but making something work which without loading DomainSets would result in an error seems different.

I totally agree with this approach, but mostly for enduser-facing packages or when pirating some rarely used function. I didn't know of heavy piracy of common basic function done by DomainSets before. Tbh, it makes me very cautious to use that package as a dependency anywhere...

In this case, DomainSets defines broadcast for Domain. If the definition of Domain moves to DomainSetsCore, so that Intervals.jl no longer "owns" it, would that settle it?

And then DomainSetsCore would define all these broadcasts? What's the benefit then, seems like the majority of DomainSets code should go there to support broadcasts.

This piracy can easily cause concrete problems, whenever there is another piece of code doing similar things. For example, I have a few definitions like Base.:*(i::Interval, x) = @modify(e -> e * x, endpoints(i)[∗]) that I copy-paste when needed to interactive notebooks. Then, (1..2) * 3 works... unless DomainSets is also loaded by some package! In that case, there's a method ambiguity error.

[aplavin]

Meanwhile DomainSets is updated to fix the promotions going on

Tried DomainSets#master, and (1..2)*1u"m" works indeed. But (1..2)u"m" still doesn't, even though it's also just a multiplication operation.

I looked at my notebooks where I define some interval operations for convenient interactive usage, and all of them are generally related to unit conversion:

Base.:*(i::Interval, x) = @modify(e -> e * x, endpoints(i)[∗])
Base.:*(x, i::Interval) = @modify(e -> e * x, endpoints(i)[∗])
Base.:/(i::Interval, x) = @modify(e -> e / x, endpoints(i)[∗])
Unitful.uconvert(u, i::Interval) = @modify(e -> uconvert(u, e), endpoints(i)[∗])
Unitful.ustrip(u::Unitful.Units, i::Interval) = @modify(e -> ustrip(u, e), endpoints(i)[∗])

@modify is from Accessors.jl, used here for simplicity – of course it's trivial to implement these functions without @modify.
I believe the meaning of these operations is straightforward and unambiguous. Is that not the case?
Any concerns regarding adding them here?

[daanhb]

I totally agree with this approach, but mostly for enduser-facing packages or when pirating some rarely used function. I didn't know of heavy piracy of common basic function done by DomainSets before. Tbh, it makes me very cautious to use that package as a dependency anywhere...

Well, let's work it out then, that's the plan anyway. To me the underlying issue is that, for historical and at the time valid reasons, IntervalSets defines the Domain type which is fundamental to the design of DomainSets. It is simply not possible to provide any functionality for domains outside of IntervalSets without affecting intervals, i.e. "piracy".

In this case, DomainSets defines broadcast for Domain. If the definition of Domain moves to DomainSetsCore, so that Intervals.jl no longer "owns" it, would that settle it?

And then DomainSetsCore would define all these broadcasts? What's the benefit then, seems like the majority of DomainSets code should go there to support broadcasts.

Indeed the solution can not be to have IntervalSets exposed to more of DomainSets. But I think there is no need for that. DomainSetsCore is semantic: it defines the notion of sets that are possibly continuous. It does not define or support anything else.

I see no problem with IntervalSets defining * or .* or any other operation using standard syntax and standard function names (i.e. defined in Base or elsewhere), as long as any operation on a continuous set results in something that behaves like a continuous set and which agrees with the mathematical semantics. Again, this is always the intention in IntervalSets anyway. It is not a problem because semantic correctness almost guaranteest interoperability and, if not, DomainSets can still adapt accordingly: having IntervalSets as a dependency its behaviour can be carefully preserved (as is the case for unions of disjoint intervals).

The benefit of DomainSetsCore is to resolve our strange dependency situation, and to allow other types to inherit from Domain without having to depend on IntervalSets. It's a minimal change now, but in time I hope that it will result in improved and clean interoperability in a much broader ecosystem.

This piracy can easily cause concrete problems, whenever there is another piece of code doing similar things. For example, I have a few definitions like Base.:*(i::Interval, x) = @modify(e -> e * x, endpoints(i)[∗]) that I copy-paste when needed to interactive notebooks. Then, (1..2) * 3 works... unless DomainSets is also loaded by some package! In that case, there's a method ambiguity error.

Hmm, if the rule is that package B cannot extend package A because package C (or a user) might also extend package A in a similar way, doesn't that exclude any extension ever?

Anyway, the ambiguity concern is genuine. To me, by inheriting from Domain a type opts into behaviour that is generically defined for domains, even elsewhere. That is the main advantage of having a common abstract supertype. If that's better known, then these errors will be rare (and always fixable).

[daanhb]

Tried DomainSets#master, and (1..2)*1u"m" works indeed. But (1..2)u"m" still doesn't, even though it's also just a multiplication operation.

Thanks for chasing this, there is still a promotion issue. Multiplication by numbers looks convenient in usage, but I now agree with @hyrodium above and think I would rather deprecate this behaviour (in favour of .*) than fixing it.

[aplavin]

DomainSetsCore is semantic: it defines the notion of sets that are possibly continuous. It does not define or support anything else.
The benefit of DomainSetsCore is to resolve our strange dependency situation, and to allow other types to inherit from Domain without having to depend on IntervalSets.

This can be useful (or not) in general, but doesn't seem related to this specific situation.
Subtyping Domain is great, as is defining functions/methods for those subtypes. But relying on type piracy for fundamental functionality is something to be avoided, especially in packages like IntervalSets that often end up deeply in dependency trees.

Hmm, if the rule is that package B cannot extend package A because package C (or a user) might also extend package A in a similar way, doesn't that exclude any extension ever?

This only applies to "extensions" that are based on piracy. And yes, it's a nice and useful rule, to avoid piracy aside from some corner cases.
Piracy here would bring, and already brings with DomainSets, changes in behavior that depend on some completely unrelated package (that happens to depend on DomainSets) being loaded. This is both extremely confusing, and can easily lead to breakage.
Suppose PkgA depends on PkgB that depends on DomainSets. PkgA uses some operation on intervals, like 2 * (1..2) – it works, everyone is happy. Then, PkgB just happens to drop the DomainSets dependency: it was only used internally, and doesn't affect their actual API. Suddenly, PkgA stops working, with errors that are hard to understand and find the fix for.

For something like IntervalSets, it doesn't seem useful or feasible to rely or encourage type piracy in these ways.

[daanhb]

For something like IntervalSets, it doesn't seem useful or feasible to rely or encourage type piracy in these ways.

To be fair, I don't think that "type piracy" is the right way to describe the situation. Since IntervalSets defines the abstract supertype Domain, it is impossible to provide functionality for domains without affecting intervals. It has been that way for years, it was a voluntary move on both ends, and it actually works quite well.

A closer analogy to our situation would be a package implementing a concrete Vector type, with an abstract supertype AbstractArray. Where does one define matrix operations, of which some also apply to vectors? This kind of problem doesn't usually arise, because people tend to implement the most generic case in the base packages first, and specialise later. From that perspective, our dependency tree is upside down.

Of course nobody argues for having DomainSets as a dependency of IntervalSets. The simple alternative is for IntervalSets to recognize that operations on continuous sets which are not specific to intervals are best implemented elsewhere. This arguably applies to the current issue of arithmetic. It more clearly applies to cartesian products such as 0..1 × 0..1, since IntervalSets has no notion of continuous sets in 2D. Yet the principle is the same.

An argument like "hard to understand and find the fix for" seems more easily managed by adding a suggestion to DomainSets somewhere in the documentation. There is certainly room for making a smaller version of DomainSets, if that helps. I'll happily add anyone here to the repo of DomainSets, by the way.

[aplavin]

Since IntervalSets defines the abstract supertype Domain, it is impossible to provide functionality for domains without affecting intervals.

Sure, and that is totally fine. As long as this functionality does not involve pirating Base methods on objects defined in IntervalSets. Otherwise, benefits and risks should be really carefully weighted.
I didn't know that DomainSets does lots of such piracy before, and it significantly changes my perception of that package.

Often, this piracy doesn't seem reasonably motivated at all. For example, I noticed it pirates LinearAlgebra.cross(interval, interval) to mean their cartesian product! Like, why do that, are there any reasonable cases when one uses cross(a, b) generically, and it makes sense for both vectors and intervals? I honestly doubt.
Also, LinearAlgebra.cross for cartesian product cannot work with Base types like arrays and sets – even though they follow the interface.
Would be better to just define and export DomainSets.×! (maybe with an ASCII synonym cartesianproduct)

Arithmetic operations like domain .+ number are somewhat more motivated. Still, I wonder if there are actual usecases for code that uses a .+ b, does not explicitly know about domains, and still makes sense for both numbers and intervals. And if these cases are strong enough to warrant piracy with its common risks.
A straightforward alternative would be to define +ₘ to mean Minkowski sum in general, for both domain +ₘ number and domain +ₘ domain.

[aplavin]

Anyway, that extended discussion is somewhat tangential to this specific issue. No matter what other packages do, I think it still makes sense to define operations that make sense for intervals and are unambiguous here. What are the drawbacks?

As a minimum, such operations include multiplication and division by a scalar.

[daanhb]

Thanks for the discussion here. Feel free to file issues for things you don't agree with.

Sure, and that is totally fine. As long as this functionality does not involve pirating Base methods on objects defined in IntervalSets. Otherwise, benefits and risks should be really carefully weighted. I didn't know that DomainSets does lots of such piracy before, and it significantly changes my perception of that package.

It does not really seem fine if you still think of it as piracy, even with Domain moving outside of IntervalSets. I understand the separate concern about using methods from Base. Returning to the issue at hand, one can wonder whether multiplication of intervals and numbers could be used generically. Thinking of a domain as a continuous set, which is mathematically well defined and differs from AbstractSet and AbstractArray mostly in that it can't be iterated over algorithmically, I do see a good case for broadcast notation to mean an elementwise operation resulting in the lazy definition of a new set.

Conceptually, iteration does kind of work with generator syntax:

julia> using DomainIntegrals

julia> integral(2x for x in 1..2)
3.0

You expressed interest in converting units. Alternatives to multiplication might be, with T = typeof(1u"m"): (ab)using generatorsT(x for x in 1..2), list comprehension T[1..2], or less intrusively map(T, 1..2) or T.(1..2).