Use DomainSets.jl?

[cscherrer]

Our domains.jl is currently kind of hacked together. When I wondered about replacing this with ManifoldsBase, @gdalle reminded me about DomainSets.jl:

Just popping by to mention https://github.com/JuliaApproximation/DomainSets.jl as part of the discussion of AbstractDomains.

Originally posted by @gdalle in #7 (comment)

Funnily enough, we had some brief discussion about this three years ago here.

There are a few places domains come up, though we don't currently have a uniform interface for managing them

1. For determining a bijection to a real vector space
2. To anticipate error conditions in logdensity and replace these cases with -Inf
3. With ⊆, for ≪ computations (nested supports are necessary for absolute continuity)
4. With ∩ for supports of pointwise product (⊙)
5. With ∪ for superposition (+)

Some potential concerns:

1. Can this play nice with ManifoldMeasures? (cc @sethaxen)
2. DomainSets currently depends on StaticArrays, which is a heavy dependency for a core library like this. Would it be possible to refactor DomainSets somewhat to use Static.jl instead, possibly with an optional StaticArrays dependency? Or maybe there's a better approach? (This part should probably be a new issue in DomainsSets, cc @daanhb @dlfivefifty)

[cscherrer]

☝️ 😕 Well, that didn't work

[gdalle]

AFAIC, the main benefit of using domains for me was being able to integrate a measure over a set. This can be useful to normalize it, or to simulate Poisson processes for instance.

[daanhb]

Just pointing out DomainIntegrals.jl here, a package for (numerically) evaluating integrals over domains, possibly with respect to a measure. The included measures are heavily biased towards approximation theory though, e.g. weight functions for orthogonal polynomials, and not very well developed. DomainSets itself also has some of that bias, to be honest.

[cscherrer]

Thanks @daanhb , DomainIntegrals looks really nice!

Off hand, challenges to address are

1. Keep MeasureBase dependencies relatively small
2. But inclusive enough to give a consistent feel to the ecosystem
3. Play nicely with packages like Manifolds.jl

We'll also need to get to things like measures on function spaces, so other packages from JuliaApproximation could well come into play.

From what I'm seeing so far, it seems like it could make sense to move toward having DomainSets as a dependency for MeasureBase. The fact that it's easily extensible should give us lots of flexibility. The main interest is not so much for its included measures (though that could be interesting as well) but as a way to describe supports of measures we build.

We could consider adding DomainIntegrals as a dependency for MeasureTheory (not MeasureBase), though we should first look at advantages of this over just keeping it separate. I mean, if we start with DomainSets, it might be easy enough for DomainIntegrals to "just work" without the need to have it as a dependency. Then packages like @gdalle is building could still take advantage of it.

[daanhb]

That sounds good. Small comment: measures are only defined in DomainIntegrals.jl, not in DomainSets. At some point they'll probably move to a package that is yet-to-be-named, or even better, be replaced with an existing package.

[cscherrer]

Simple proof of concept here:
#19

There's really not much to it, more interesting will be to see how this change would affect MeasureTheory.jl

[cscherrer]

This is really helpful (thanks @daanhb ):
JuliaApproximation/DomainSets.jl#91 (comment)

So one possibility would be for a Measure to be parameterized by a Domain. It would be really nice to not have to worry about domains, and to have a clean connection to the JuliaApproximation tools. My biggest concern at this point is with how this could affect @sethaxen's work on
https://github.com/JuliaManifolds/ManifoldMeasures.jl

Daan or Seth (or others), what do you think?

[sethaxen]

So one possibility would be for a Measure to be parameterized by a Domain. It would be really nice to not have to worry about domains, and to have a clean connection to the JuliaApproximation tools. My biggest concern at this point is with how this could affect @sethaxen's work on
https://github.com/JuliaManifolds/ManifoldMeasures.jl

I haven't really been following this discussion, but this is similar to the solution I have been using for ManifoldMeasures.jl. There, each measure is associated with a manifold. Because some measures can be written generically for multiple manifolds (e.g. VonMisesFisher), many of the measures actually store a manifold object, but this is not currently part of the API; instead, we overload the Manifolds.base_manifold method from Manifolds to retrieve the manifold.

If one wanted to do the same thing with MeasureTheory but using Domains, this would be fine unless core functionality of MeasureTheory unnecessarily relied on a domain being defined as as a subtype of some type in DomainSets. An option is to now have a type constraint, which could allow swapping a Manifold as a "domain".

[cscherrer]

If one wanted to do the same thing with MeasureTheory but using Domains, this would be fine unless core functionality of MeasureTheory unnecessarily relied on a domain being defined as as a subtype of some type in DomainSets. An option is to now have a type constraint, which could allow swapping a Manifold as a "domain".

If I'm understanding correctly, this should be ok as long as every Manifold has in andeltype methods. Here's from the bottom of the DomainSets README:

The domain interface

A domain is any type that implements the functions eltype and in. If d is an instance of a type that implements the domain interface, then the domain consists of all x that is an eltype(d) such that x in d returns true.

Domains often represent continuous mathematical domains, for example, a domain d representing the interval [0,1] would have eltype(d) == Int but still have 0.2 in d return true.

The Domain type

DomainSets.jl contains an abstract type Domain{T}. All subtypes of Domain{T} must implement the domain interface, and in addition support convert(Domain{T}, d).

DomainSets seems to have very nice handling of unions and products. So I think we can have some simple laws like

the domain of a ProductMeasure is the ProductDomain of its components

and

the domain of a SuperpositionMeasure is the union of its components.

Other interesting things:

• DomainIntegrals includes some measures, so hopefully we can interface cleanly with this
• Lots of the JuliaApproximation packages deal with infinite arrays and function spaces. I'm not sure yet about how these connect to DomainSets, but we'll definitely want measures on spaces like these

[daanhb]

Just to confirm, we like to avoid requiring domains be a subtype of Domain. This is not always convenient or possible, but it does lead to flexibility where it works. Even a vector is a domain, because it defines membership. I'd support a design in which the Domain type is not explicitly assumed, only something that behaves like a domain. I did not yet look at the manifolds being mentioned here.

[sethaxen]

If I'm understanding correctly, this should be ok as long as every Manifold has in andeltype methods. Here's from the bottom of the DomainSets README:

While all Manifolds have an in method, they do not have an eltype method by design, since there are many potential representations for points on a manifold, including different array types one might want to use as well as different coordinate systems. We of course have a documented "default representation," but even that is never type-constrained, it's simply assumed to be the argument unless one adds a custom method for a different representation.

In ManifoldMeasures, for sampling, we have had to make this more explicit by defining a default_point, but that is just to generate a point in rand that is overwritten by rand!. One can pass a different point type to rand! if they like.

[cscherrer]

Probably worth noting again that Distributions.jl has very weird eltype behavior:

julia> using Distributions

julia> eltype(Normal())
Float64

julia> eltype(MvNormal(ones(3)))
Float64

This doesn't compose well, but it's probably much to late to fix. I had thought this meant we should avoid eltype (this led us to sampletype), but maybe we should have it after all? It seems easy enough to have it be abstract, for example

julia> eltype(AbstractArray{AbstractVector{Real}})
AbstractVector{Real} (alias for AbstractArray{Real, 1})

Also, there's a default method:

julia> abstract type T end

julia> eltype(T)
Any

So in the DomainSets requirement,

the domain consists of all x that is an eltype(d) such that x in d returns true.

I guess this just becomes a universal quantifier?