Add support for domains with vector elements

[daanhb]

This pull request is a work in progress. It started by addressing #53. It also (partially) addresses #44 and #43.

1) We add support for balls, spheres and simplices with Vector elements, e.g.:

julia> d = VectorUnitBall(10)
the 10-dimensional closed unit ball

julia> eltype(d)
Array{Float64,1}

julia> rand(10) ∈ d
false

Similar to intervals, balls and simplices can now be closed or open.

A product domain has a Vector element type by invoking the constructor with a vector of domains (rather than a list of arguments or a tuple of domains):

julia> ProductDomain([0..rand() for i in 1:10])
0.0..0.24528845369688734 x 0.0..0.6088588002463498 x 0.0..0.554007387962554 x 0.0..0.8423193868801728 x 0.0..0.5512596358694362 x 0.0..0.08741952576511558 x 0.0..0.024360588101769665 x 0.0..0.46943098768146463 x 0.0..0.6181070274644074 x 0.0..0.2949387321672272

julia> eltype(ans)
Array{Float64,1}

2) A ProductDomain has two constructors: the unqualified ProductDomain attempts to simplify the element type to an SVector as before, but ProductDomain{T} simply makes a product domain with the specified element type T. The translation between nested tuples and flat vectors is now done explicitly in productdomain.jl. As a result, the machinery of spaces and embeddings could be removed from this package, addressing #43 and simplifying the code. Rather than storing an "internal and external" element type, the ProductDomain now keeps the dimensions of the domains:

julia> p = cylinder()
the 2-dimensional closed unit ball x 0.0..1.0 (Unit)
julia> typeof(p)
ProductDomain{SArray{Tuple{3},Float64,1,3},(2, 1),Tuple{DomainSets.FixedUnitBall{SArray{Tuple{2},Float64,1,2},:closed},UnitInterval{Float64}}}

The type parameter DIM is (2,1) in this example, which is enough information to do the conversion from a vector [1,2,3] to ([1,2], 3) and back efficiently and generically.

3) Addressing #44, the code is simplified somewhat by considering ProductDomain, WrappedDomain, DifferenceDomain, UnionDomain and IntersectionDomain as subtypes of a LazyDomain. A LazyDomain implements its in method in three steps:

• preprocess the x (this can be about applying a map, or converting between [1,2,3] and ([1,2],3),
• distributing x over its member domains, invoking their in function
• combining the output of in of the member domains to a final result. For example, a UnionDomain will pass x unaltered to all its member domains, and do a mapreduce on the outcomes with |. A ProductDomain will convert x from flat to nested format, distribute the elements of x separately to its member domains, and combines the outcomes with &. This does not go as far as #44 proposes, but it gets us closer.

4) Currently, the in implementation of the generic Domain attempts to do some conversions as before (possibly promoting both x and the domain), but it will also automatically do a limited number of conversions. It allows any AbstractVector for a VectorDomain (which is a Domain{Vector{T}}). It allows any AbstractVector for a EuclideanDomain (which is a Domain{SVector{N,T}}) by converting to an SVector. It also allows an NTuple{N,T} for a EuclideanDomain{N,T}. It returns false if the types do not match, but currently prints a warning if that happens.

[daanhb]

The errors in Julia 1.0 and 1.1 are related to mapreduce. A ProductDomain previously had some hardcoded implementations of in for a small number of dimensions and a reduce(&, map(...)) fallback for higher dimensions that would allocate memory. I found that mapreduce works allocation-free and I could remove the special cases, but apparently that is only possible in more recent Julia.

[daanhb]

We're getting there, I'm writing some more tests. Some preliminary conclusions:

1. I've rewritten product domains some more, more closely following the suggestion of #44. The notion of a VcatDomain nicely captures the "flattening magic" of the previous ProductDomain constructor. There are now several product domains, depending on whether the points of the member domains should be flattened to a static vector, collected in a Vector, or collected in a tuple. The first two promote their member domains to a common numeric type, the latter one doesn't. The ProductDomain constructor still exists and decides on which concrete domain to return based on the types of the arguments, but VcatDomain, VectorProductDomain and TupleProductDomain can also be invoked directly. The tuple product domain is the most general one of the three:

   
   julia> p = ProductDomain(1:10, [true,false], 0..1)
   DomainSets.WrappedDomain{Int64,UnitRange{Int64}}(1:10) x DomainSets.WrappedDomain{Bool,Array{Bool,1}}(Bool[1, 0]) x 0..1

julia> eltype(p) Tuple{Int64,Bool,Int64}

julia> (4, true, 0.5) ∈ p true

julia> p = ProductDomain([0..i for i in 1:5]) 0..1 x 0..2 x 0..3 x 0..4 x 0..5

julia> typeof(p) VectorProductDomain{Int64,Interval{:closed,:closed,Int64}}

julia> p = ProductDomain(UnitDisk(), UnitDisk()) the 2-dimensional closed unit ball x the 2-dimensional closed unit ball

julia> eltype(p) SArray{Tuple{4},Float64,1,4}

2. Identifying tuples with an `SVector` was not a good idea because, well, `(0,1) != SVector(0,1)`. I included it at first but then removed it. On the other hand, identifying abstract vectors with other vectors is a good idea, because e.g. `1:2 == SVector(1,2)`.

3. Compared to #52, the `in` function in this pull request behaves like `in` already did. The notion of a domain remains that `x` is in `d` if there is a `y` in `d` such that `x==y`, and this very often requires promotion of both `x` and `d`. Two differences are that fewer conversion errors are thrown (owing to more careful treatment of vectors and the way they promote) and in a few cases a warning is displayed before `false` is returned to flag a likely error (e.g. if the lenghts of vectors don't match).

[daanhb]

I think this pull request is ready. Can you check @dlfivefifty for unexpected compatibility problems?

I've tried to minimize breaking changes (exports etcetera). Support for vectors is much improved, and in particular the issue leading to #58 should work. One conflict with ApproxFun is probably that the boundary is redefined (the boundary of an open interval now contains the endpoints even if they are not included in the domain, see #56). This can be reverted for the time being if that's a problem.

[dlfivefifty]

This is creating several strange test failures and I can't figure out why....

[daanhb]

I ran the ApproxFun test suite a while back and concluded at the time that at least the definition of the boundary made a difference, but I did not pursue it. In contrast I've had to make only very minor changes in my own packages because functionally nothing changed much, it was related to syntax change of spheres. Alright, I'll look into it, but it may not be right away.

[dlfivefifty]

Can you fix the conflicts and I'll try again?

[daanhb]

I've looked at it but I don't think there is anything I can merge. The added test works on the "vector" branch:

julia> using StaticArrays, DomainSets

julia> struct Basis3Vector <: StaticVector{3,Float64} end

julia> Base.getindex(::Basis3Vector, k::Int) = k == 1 ? 1.0 : 0.0

julia> Basis3Vector() in UnitSphere()
true

The conversion in circle.jl that was extended no longer exists. To avoid hardcoding static vectors everywhere EuclideanHyperSphere{N,T} is now just an alias for a HyperSphere with element type SVector{N,T}.

How did you arrive at a StaticVector? Perhaps another issue was at play here: sometimes an SVector was converted into a mutable MArray, because of the way we use conversion to AbstractVector{T} to update an element type. I've chased this issue both in StaticArrays (see https://github.com/JuliaArrays/StaticArrays.jl/issues/746) and in Julia itself, but the fix is pending. On the vector branch I've temporarily overridden this behaviour. Could that be how you ended up with a StaticVector rather than an SVector?

[dlfivefifty]

I ended up with a static vector because I made a special type to represent a point in spherical coordinates:

https://github.com/JuliaApproximation/SphericalHarmonics.jl/blob/f296d8b976c3ea2e33f8cc069584cc9b4c5b5367/src/SphericalHarmonics.jl#L24

[daanhb]

That's good to know of, I did not check for that use case.