Rethink `Domain2d`

[dlfivefifty]

Heterogenous types (I.e., Domain{Tuple{Int,Float64}}. are also 2-dimensional domains.

[daanhb]

Currently, Domain2D{T} is short for EuclideanDomain{2,T}. It is what I use most and it simply felt silly to make a longer shorthand such as EuclideanDomain2{T}, even though that would be a more correct name.

In general, it is not so easy to define the dimension of a domain. That is because we expect the dimension of a domain with SVector{2,Float64} elements to have dimension 2. On the other hand, the convention is that product domains have a tuple type, like your example Tuple{Int,Float64}. But that doesn't mean that Domain{Tuple{A,B}} always has dimension two. If we accept that a single type T can have some dimension higher than 1, like it does in the SVector{2,Float64} case, things become fuzzy: what is the dimension of Domain{Tuple{SVector{2,Int},SVector{2,Float64}}}? Is it two or four?

Both answers lead to problems: the notion of dimension itself is simply not sufficiently flexible to allow for the genericity we have. I'd avoid using it, except in cases where we clearly know what it means.

There are several consistent definitions of a dimension:

1. The static one. Domain{T} has dimension 1, unless T is a tuple type. In that case, its dimension is the length of the tuple. This is the easiest to get right generically, but it leads to the weird conclusion that the dimension of any homogeneous domain is 1.
2. The recursive one. Domain{T} has some dimension associated with T (as returned by dimension_type(T)). We acknowledge that there exist "container types" that have a dimension higher than 1. The dimension when T is a tuple type is the sum of the dimensions of the elements of the tuple. This is difficult to implement: the seemingly simple operation of iterating over the dimensions can be really hard. (This is solved by writing a proper isomorphism to a simpler space).
3. The current one. The dimension of a domain is only defined for some types where we clearly know what the user is expecting, and it fails with an error in other cases. The use of "dimension" is considered harmful in generic code. I guess the case Domain{Tuple{S,T}} where S and T are both numbers could be added. The user can get the answer he wants by defining dimension_type for his preferred type.

[dlfivefifty]

The only usage for dimension I can think of is number of function arguments, so Tuple{SVector{2,Int}, SVector{2,Float64}} would be 2-dimensional, since it would work for functions of the form f(SVector(1,2), SVector(2,3))

[daanhb]

That would lead to option 1 above, and possibly the alias Domain2D{S,T} = Domain{Tuple{S,T}} or something of the sort. It is indeed consistent. I wonder whether it is very useful, i.e., do we even need to bother defining this alias. It still seems confusing to me to define dimension this way. On the other hand, it still seems very useful to me to define the dimension of a Euclidean domain (perhaps with a different name, I don't particularly mind).

[MikaelSlevinsky]

If dimension is determined by the evaluation signature, which usually calls in, then the 2-sphere's dimension is 3. But what if a function is evaluated on the sphere with two arguments? Should they be interpreted as parametric arguments or the third ghost argument of 0? or throw an error?

Should there be a distinction between ambient dimension and parametric dimension?

[daanhb]

As things stand, a domain has a unique element type T and that is the type that in expects. So a sphere would have to choose between ambient space and parametric space. Of course any domain is free to add more implementations of in for convenience. It is T that determines dimension.

Conversion between spaces can be handled by defining maps between the spaces (by defining an embedding from the sphere to the ambient space and a restriction in the other direction). That makes it possible to define a 2-sphere with a 2d parameterization, and express that it is embedded in 3D. That would actually be a nice test case for this mechanism.

[daanhb]

The only usage for dimension I can think of is number of function arguments, so Tuple{SVector{2,Int}, SVector{2,Float64}} would be 2-dimensional, since it would work for functions of the form f(SVector(1,2), SVector(2,3))

An argument against that is that you can make a function that accepts a single argument of type Tuple{SVector{2,Int},SVector{2,Float64}}. It only becomes a function of two variables if you splat the tuple. But you can also splat a vector. Perhaps the dimension should be the number of elements that you get when something of type T is splatted?

[MikaelSlevinsky]

It is T that determines dimension.

So in this framework, the 2-sphere could have dimension 3 if the input to evaluation contains three terms (T == Tuple{3,S} where S) or dimension 2 if input to evaluation is a parametric representation with (as we know, the minimal requirement of) two terms, and there is conversion between the different representations of the 2-sphere.

But the 2-sphere has no thickness! How can a 3-dimensional 2-sphere be legitimate?

[MikaelSlevinsky]

I'm sorry, I don't mean to be critical. It's just that as some of my repositories are downstream, I'll de facto be a user of this repository, so I may as well ask my questions now, even if they're silly.

[MikaelSlevinsky]

With the minimal parameterization of the 2-sphere, in is practically free. This is not the case for the Cartesian parameterization, which is more expensive.

[dlfivefifty]

A domain implements In and the only definition that makes any sense for a sphere is norm(x) = 1.

You can make your own domain TopologicalSphere that is a rectangle with the topology of a sphere

[MikaelSlevinsky]

Ok, but what does the 2 stand for in 2-sphere?

Cheers,

Mikael

On Oct 19, 2018, at 11:21 AM, Sheehan Olver notifications@github.com<mailto:notifications@github.com> wrote:

A domain implements In and the only definition that makes any sense for a sphere is norm(x) = 1.

You can make your own domain TopologicalSphere that is a rectangle with the topology of a sphere

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/JuliaApproximation/Domains.jl/issues/31#issuecomment-431418874, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AHzBpZm9oz0rqdhLndNxEIIsarUa6KlSks5umfwIgaJpZM4XjxXk.

[dlfivefifty]

https://en.wikipedia.org/wiki/Sphere

Here dimension is used for the dimension of the ambient space, not the dimension of the surface.

Maybe dimension(::Domain) should be removed in favour of dimension(eltype(d)) to avoid confusion.

[dlfivefifty]

Or renamed nargs(d)

[daanhb]

Now is a better time to ask questions than later :-)

The current interpretation is that a Domain{T} is a collection of elements of type T. That is pretty much it. That leads to norm(x) == 1 as @dlfivefifty indicates: it is the only sensible implementation of in. It is also, I'm afraid, practically useless due to rounding errors (that is why we also have approx_in).

The circle is a simpler example. If you think of it as the map of [0,1] using a cos+sin parameterisation then you can have it your way too. This actually works with the latest master:

julia> using Domains, StaticArrays

julia> param_domain = 0.0..1.0
0.0..1.0

julia> circle = Domains.BidirectionalMappedDomain(param_domain, Domains.UnitCircleMap{Float64,SVector{2,Float64}}(), Domains.AngleMap{SVector{2,Float64},Float64}())
A mapped domain based on 0.0..1.0

julia> dimension(param_domain)
1

julia> dimension(circle)
2

julia> in(SVector(cos(pi/3), sin(pi/3)), circle)
true

The syntax isn't exactly pretty and I'm open to suggestions on how to improve things! But now the in for the circle is based on mapping the point in the plane back to [0,1] and checking for membership there. It actually allows for some rounding errors in this case (it would be nearly impossible not to).

[MikaelSlevinsky]

@dlfivefifty , I'm reading the Wikipedia page, and perhaps you can be more precise?

For example, in the third paragraph it says "...a sphere, which is a two-dimensional closed surface, embedded in a three-dimensional Euclidean space..."

[MikaelSlevinsky]

@daanhb , that example is sensible except when it comes to dimension(circle) == 2. I understand the reasoning, based on the parameter T, and it's certainly easier from a software point of view. But it's just counterintuitive to me.

Why can't a circle always be a parameterized domain, and if in(x, d) receives an SVector{2,T} as first argument, it can dispatch to checking norm(x) == 1 or a reasonable approximation?

[dlfivefifty]

Well just do

` struct RadialCoordinate{T} <: AbstractStaticVector{2,T}
   r
   theta
end

norm(z::RadialCoordinate) = z.r

[daanhb]

@MikaelSlevinsky The goal of the game is consistency: generic code only works if there is a consistent semantic meaning of all the types involved. If you think of the sphere as part of 3D space, then use T = SVector{3,Float64}. If you think of it as being 2D, then you implement the parameterisation and choose T=SVector{2,Float64}. If you want to map between the two, then you implement the maps (like UnitCircleMap and AngleMap in the example above). They are not elegant, but they simply reflect the complexity of the geometric situation. You can interpret a sphere in a million different ways. The best we can do is provide software such that you can implement your particular point of view (and make the least number of assumptions about it, and if we do make assumptions, be very clear about them). Also, we can provide methods to convert between different software representations of what is mathematically the same thing.

Here is another problem. The domains package allows you to make unions of domains. What if you want to represent the union of a sphere and a point? You can only do that in 3D. Hence, even the dimension of the point would be 3 (because it has type T=SVector{3,Float64}). What is the dimension of a point? It depends on T. Perhaps @dlfivefifty is right and we shouldn't be calling it the dimension of the domain. What we really mean is the dimension of T.

[daanhb]

Well just do ` struct RadialCoordinate{T} <: AbstractStaticVector{2,T} r theta end

norm(z::RadialCoordinate) = z.r

Cool example! Yes, that suddenly makes the implementation of in for spheres more realistic. The way you see a domain is really encoded into the type T.

[MikaelSlevinsky]

Yes, I agree: using the term "dimension" to refer to the length of T, which is unambiguous, is perhaps not a good choice of name. nargs (or perhaps some other alternative, but not dimension) is better, IMO.

[daanhb]

Note for the record that I'm still struggling too with the current implementation of a sphere, even as I try to defend it. The only interface for a Domain we ask is an implementation of in, and even for the important special case of a sphere that implementation isn't particularly helpful. I guess it really comes down to a practical choice. Having a simple consistent interpretation helps. Having a single type parameter T also helps in writing type-stable code.

Trying to shoehorn all interpretations of a sphere into one Julia type seems infeasible. Two domains are equal if they describe the same set of points -- period. That is simple. They do not necessarily have the same type, we'll have to live with that. It means that you can't rely on dispatch to catch all possible spheres (even an interval is a 1d sphere). You could have an issphere method, or you could write a program that only ever uses one representation and dispatch on that.

[dlfivefifty]

There’s the same issue in 1D: it’s much better to represent points in [-1,1] as

struct CosPoint{T} <: Real
theta
end
convert(::Type{T} , x::CosPoint) = cos(x.theta)

This makes Clenshaw more accurate

[daanhb]

Did you use that somewhere yet? That is hilarious :-)

[MikaelSlevinsky]

Based on https://en.wikipedia.org/wiki/Dimension, the dimension of an object "is the minimum number of coordinates needed to specify any point within it."

So, @daanhb, if I take a stab at your union domain problem of a point and a sphere, I'd say the dimension is 2, because the point is 0-dimensional (in any number of ambient dimensions).

By the sounds of this conversation, I think we've moved on from claiming to support the notion of dimension in the code, but I could be wrong.

[daanhb]

@MikaelSlevinsky That makes sense, thanks. I think we've agreed that all references to "dimension" in the code actually refer to the dimension of T, i.e., of the ambient space, not to that of the domain. I think it is difficult to generically discover the dimension of a domain when given only the implementation of in, especially after a few manipulations. It seems the definition Domain{T} naturally leads to a setting with an ambient space, which we arrived at because it is practical.

It remains to be decided what to do with the Domain2D alias which started this discussion.

There could be a way for the user (or the implementer of a type) to provide more information such as dimension. Now that the Julia infrastructure is there, I'd like to look into the possibility of having some methods return true, false or Missing if we don't know the answer for sure. That would allow us to return results when we're sure, while never returning incorrect answers.

I also think we should implement a Klein bottle, algebraic varieties, spacetime, and other funny topological domains with strange parameter types, just to make sure our implementation is sufficiently generic to cover everything under the sun.

[dlfivefifty]

I think those make more sense in a Topologies.jl setting: A Möbius strip is probably best represented as a rectangle with a special topology, since an embedding is artificial.

Spheres are different since they are “natural” embedded in R^3, the debate is really about representing points in R^3 on the sphere, not about the sphere (and in ) representation.

I guess a topology would be a able to generate an open set around any point in the domain.

[MikaelSlevinsky]

Did you use that somewhere yet? That is hilarious :-)

P.S. it's not only that Clenshaw gets more accurate, it's that you can come up with different evaluation algorithms for angles near the endpoints that also work best with angular input. (Section 4.2, https://arxiv.org/abs/1602.02618)

[MikaelSlevinsky]

@dlfivefifty , I feel like there should be a "CoordinateSystems.jl" package in Julia. I guess this one has the name but has been untouched for a while https://github.com/drewkett/CoordinateSystems.jl (and is not very fleshed out).

Then, we could represent the 11 (resp. 13) known coordinate systems where the Helmholtz (Laplace) equation is separable (http://mathworld.wolfram.com/LaplacesEquation.html).

Then, in for domains could dispatch to all sorts of fun coordinates!

[MikaelSlevinsky]

Looking at the code and returning to the first issue at hand, I'd suggest to remove exports of controversial aliases. As presented, someone else may have another use for Domainxd in their code and it seems unreasonable to me to increase/impose on the namespace.

This shouldn't be a big issue because they're not that frequently used internally either. Also, I don't think it's type piracy to declare type aliases in a different repository (I could be wrong), so even if the aliases figured prominently in dependent repositories, the aliases could be transferred.

[MikaelSlevinsky]

Re: ambient dimension, the 2-sphere may be embedded in Rⁿ for any n>2, and arguably the type information available in T is insufficient for a generic implementation of in. In R⁴, for example, is the 2-sphere embedded in the first three of four coordinates? or more generically in some three-dimensional subspace unaligned with any of the canonical coordinate axes?

I think this is an argument for keeping track of canonical parameterizations and topologies instead of the ambience.

[daanhb]

With the release behind us, let's just remove Domain2D and friends?

[dlfivefifty]

Now that the name of the package is settled it's a lot easier to tag new releases, we can always have a few 0.0.x releases with big changes. I think the goals for 0.1 should be:

• [ ] Remove Domain2d and friends since its really a property of T.
• [ ] Remove unnecessary restrictions to ::Domain: one should be able to do UnionDomain([1,2,3], 4..5, 6) since each of the arguments supports in
• [ ] Come up with an alternative to convert(Domain{T}, d) that supports [1,2,3], 6, etc. Something like promote_eltype(T, d) (unfortunately Base.promote_eltype exists but does something different so that's not a good name.)

For now I need to get ApproxFun (depending on DomainSets.jl), BlockBandedMatrices, etc. updated and tagged so I can fix some outstanding issues, so I won't be able to work on this for a while.

[daanhb]

Meanwhile, Domain2d and friends have been removed in master with commit c044484. I did not add any alternative to compute the dimension of a domain: this is indeed a property of T and any user can define the dimension of T however they like.

[dlfivefifty]

Sounds good. We'll have to put an upper bound on ApproxFun in METADATA but that's fine.

There's two options for the next version:

1. 0.0.2, since both Semvar and "Julia Semvar" say that 0.0.x is any changes go
2. 0.1.0, which is standard Semvar for any changes go, but in "Julia Semvar" means versions with breaking changes should be 0.2.0. So one should only tag 0.1.0 when things are somewhat settled.

[daanhb]

Hmm, good point, I did not consider how breaking this was. Our code elsewhere relies on Domain2d being absent (we redefined it in another package), but we can change that. This seems like too small a change to warrant a new release. Revert for the time being?

[dlfivefifty]

Just leave it and use the tagged version in other packages

[daanhb]

Ok, I'll leave it in.

About the alternative to convert(Domain{T}, d): as you mention we can't use promote_eltype because it is already defined. The name does indicate what we would want it to be: something that allows to make sure that different objects that implement in use a common element type. What about promote_space?

I can have a go, we can change the name later.

[daanhb]

I went with convert_domain and promote_domain, see #33 for a start.

The goal is not to promote the types of the domains themselves, but the element type: it turns out this is easy to implement in terms of the existing promotion system.

We now have:

julia> DomainSets.promote_domain([1,2,3], 6)
([1, 2, 3], 6)

julia> DomainSets.promote_domain([1,2,3], 6.0)
([1.0, 2.0, 3.0], 6.0)

julia> DomainSets.promote_domain([1,2+im,3], Point(6.0))
(Complex{Float64}[1.0+0.0im, 2.0+1.0im, 3.0+0.0im], Point{Complex{Float64}}(6.0 + 0.0im))

[daanhb]

Note for completeness that the spaces are for different things: they allow to easily express that you want to identify certain types even if Julia doesn't. For example, identifying complex numbers with the plane. We probably do not want this by default, and it is cleaner here to just rely on the existing promotion system. Hence, I did not pursue using promote_space.

[dlfivefifty]

Be careful here since its convention to treat 1+0im and 1+0.0im as "real", but no errors are allowed:

julia> (1+0im) ∈ 0:5
true

julia> (1+0.0im) ∈ 0:5
true

julia> (1+nextfloat(0.0)im) ∈ 0:5
false

[daanhb]

I don't think this is an issue currently: Julia itself also promotes Complex{Int} to Complex{Float64} if you mix it with floats. (Or am I missing your point?)

Most importantly, any promotion that takes place here does not change the set:

julia> a, b = DomainSets.promote_domain([1,2+im,3], Point(6.0))
(Complex{Float64}[1.0+0.0im, 2.0+1.0im, 3.0+0.0im], Point{Complex{Float64}}(6.0 + 0.0im))

julia> 1 ∈ a
true

julia> (1+nextfloat(0.0)im) ∈ a
false

[dlfivefifty]

Most importantly, any promotion that takes place here does not change the set:

Great! that's my point I think.

[daanhb]

We've removed the DomainXD aliases, so I'll close this issue. The code does have a function dimension now which is non-recursive and corresponds to the length of a static vector or the length of a tuple, mostly for use in product domains. If there are any issues with that definition, we can discuss it separately.