Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implementing domains as an interface #141

Merged
merged 13 commits into from
Oct 9, 2023
8 changes: 4 additions & 4 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "DomainSets"
uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
version = "0.6.7"
version = "0.7.0"

[deps]
CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"
Expand All @@ -11,15 +11,15 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[compat]
julia = "1.6"
CompositeTypes = "0.1.2"
IntervalSets = "0.7.4"
StaticArrays = "0.12.2, 1"
StableRNGs = "1"
julia = "1.6"
StaticArrays = "0.12.2, 1"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test", "StableRNGs"]
Empty file added docs/src/interface.md
Empty file.
6 changes: 4 additions & 2 deletions src/DomainSets.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,6 @@ export ..

import CompositeTypes: component, components


################################
## Exhaustive list of exports
################################
Expand Down Expand Up @@ -81,14 +80,16 @@ export AffineMap, Translation, LinearMap,

## Generic domains

@deprecate point_in_domain(d) choice(d)

# from generic/domain.jl
export Domain, EuclideanDomain, VectorDomain,
dimension,
approx_in,
isopenset, isclosedset, iscompact,
boundary, ∂,
interior, closure,
point_in_domain
choice

# from generic/geometry.jl
export boundingbox,
Expand Down Expand Up @@ -192,6 +193,7 @@ include("maps/basic.jl")
include("maps/affine.jl")
include("maps/arithmetics.jl")

include("generic/core.jl")
include("generic/domain.jl")
include("generic/geometry.jl")
include("generic/canonical.jl")
Expand Down
8 changes: 4 additions & 4 deletions src/domains/ball.jl
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ isequaldomain(d1::Ball, d2::Ball) = isclosedset(d1)==isclosedset(d2) &&
radius(d1)==radius(d2) && center(d1)==center(d2)
hash(d::Ball, h::UInt) = hashrec("Ball", isclosedset(d), radius(d), center(d), h)

function issubset(d1::Ball, d2::Ball)
function issubset_domain(d1::Ball, d2::Ball)
if dimension(d1) == dimension(d2)
if center(d1) == center(d2)
if radius(d1) < radius(d2)
Expand Down Expand Up @@ -84,7 +84,7 @@ normal(d::Ball, x) = normal(boundary(d), x)
distance_to(d::Ball, x) = x ∈ d ? zero(prectype(d)) : norm(x-center(d))-radius(d)

# We choose the center of the ball here. Concrete types should implement 'center'
point_in_domain(d::Ball) = center(d)
choice(d::Ball) = center(d)


"The unit ball."
Expand Down Expand Up @@ -119,7 +119,7 @@ approx_indomain(x, d::ClosedUnitBall, tolerance) = norm(x) <= 1+tolerance
isequaldomain(d1::UnitBall, d2::UnitBall) = isclosedset(d1)==isclosedset(d2) &&
dimension(d1) == dimension(d2)

issubset(d1::UnitBall, d2::UnitBall) =
issubset_domain(d1::UnitBall, d2::UnitBall) =
dimension(d1) == dimension(d2) && (isclosedset(d2) || isopenset(d1))

convert(::Type{SublevelSet}, d::UnitBall{T,C}) where {T,C} =
Expand Down Expand Up @@ -304,7 +304,7 @@ normal(d::Sphere, x) = (x-center(d))/norm(x-center(d))

distance_to(d::Sphere, x) = abs(norm(x-center(d))-radius(d))

point_in_domain(d::Sphere) = center(d) + unitvector(d, 1)
choice(d::Sphere) = center(d) + unitvector(d, 1)

isequaldomain(d1::Sphere, d2::Sphere) =
radius(d1)==radius(d2) && center(d1)==center(d2)
Expand Down
34 changes: 18 additions & 16 deletions src/domains/boundingbox.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,52 +3,54 @@ boundingbox(d::Vector{T}) where {T <: Number} = minimum(d)..maximum(d)
boundingbox(d::Set{T}) where {T<:Number} = minimum(d)..maximum(d)

"Return the bounding box of the union of two or more bounding boxes."
unionbox(d::Domain) = d
unionbox(d1::Domain, d2::Domain) = unionbox(promote_domains(d1, d2)...)
unionbox(d1::Domain, d2::Domain, domains...) =
unionbox(d) = d
unionbox(d1, d2) = unionbox1(promote_domains(d1, d2)...)
unionbox(d1, d2, domains...) =
unionbox(unionbox(d1,d2), domains...)

unionbox(d1::Domain{T}, d2::Domain{T}) where {T} = unionbox1(d1, d2)
unionbox1(d1, d2) = unionbox2(d1, d2)
unionbox2(d1, d2) = fullspace(d1)
unionbox1(d1::EmptySpace, d2) = d2
unionbox1(d1::FullSpace, d2) = d1
unionbox2(d1, d2::EmptySpace) = d1
unionbox2(d1, d2::FullSpace) = d2

unionbox(d1::D, d2::D) where {D<:FixedInterval} = d1
unionbox1(d1::D, d2::D) where {D<:FixedInterval} = d1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's unionbox1 and below intersectbox1?

Copy link
Member Author

@daanhb daanhb Sep 5, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

unionbox(d1,d2) is a function of two arguments, which sometimes simplifies for certain d1 (independently of d2 or for some possibilities of d2). Similarly, it may simplify for certain d2. The call chain goes:
unionbox(d1,d2) -> unionbox1(d1,d2) -> unionbox2(d1,d2) -> some_default_algorithm

The idea is that unionbox1 can dispatch safely (without ambiguity) solely on d1, and unionbox2 can dispatch on d2. I've used this pattern in many places and it works rather well. If a combination of d1 and d2 is specific enough, then one can dispatch on it using just the unionbox function.


function unionbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T}
function unionbox1(d1::AbstractInterval, d2::AbstractInterval)
a, b = endpoints(d1)
c, d = endpoints(d2)
A = min(a, c)
B = max(b, d)
isinf(A) && isinf(B) ? fullspace(T) : A..B
isinf(A) && isinf(B) ? fullspace(d1) : A..B
end

unionbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T} =
function unionbox(d1::HyperRectangle, d2::HyperRectangle)
@assert dimension(d1) == dimension(d2)
T = promote_type(eltype(d1),eltype(d2))
Rectangle{T}(map(unionbox, components(d1), components(d2)))
end

"Return the bounding box of the intersection of two or more bounding boxes."
intersectbox(d::Domain) = d
intersectbox(d1::Domain, d2::Domain) = intersectbox(promote_domains(d1, d2)...)
intersectbox(d1::Domain, d2::Domain, domains...) =
intersectbox(d) = d
intersectbox(d1, d2) = intersectbox1(promote_domains(d1, d2)...)
intersectbox(d1, d2, domains...) =
intersectbox(intersectbox(d1,d2), domains...)

intersectbox(d1::Domain{T}, d2::Domain{T}) where {T} = intersectbox1(d1, d2)
intersectbox1(d1, d2) = intersectbox2(d1, d2)
intersectbox2(d1, d2) = fullspace(d1)
intersectbox1(d1::EmptySpace, d2) = d1
intersectbox1(d1::FullSpace, d2) = d2
intersectbox2(d1, d2::EmptySpace) = d2
intersectbox2(d1, d2::FullSpace) = d1

intersectbox(d1::D, d2::D) where {D<:FixedInterval} = d1
intersectbox1(d1::D, d2::D) where {D<:FixedInterval} = d1

intersectbox(d1::AbstractInterval{T}, d2::AbstractInterval{T}) where {T} =
intersectbox1(d1::AbstractInterval, d2::AbstractInterval) =
intersectdomain(d1, d2)

function intersectbox(d1::HyperRectangle{T}, d2::HyperRectangle{T}) where {T}
function intersectbox1(d1::HyperRectangle, d2::HyperRectangle)
@assert dimension(d1) == dimension(d2)
T = eltype(d1)
d = Rectangle{T}(map(intersectbox, components(d1), components(d2)))
isempty(d) ? emptyspace(d) : d
end
Expand Down
4 changes: 2 additions & 2 deletions src/domains/cube.jl
Original file line number Diff line number Diff line change
Expand Up @@ -351,8 +351,8 @@ ProductDomain(domains::NTuple{N,D}) where {N,D <: FixedInterval} =
ProductDomain(domains::SVector{N,<:FixedInterval}) where {N} =
FixedIntervalProduct(domains)
ProductDomain{T}(domains::D...) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
FixedIntervalProduct(convert_eltype.(Ref(S), domains))
ProductDomain{T}(domains::NTuple{N,D}) where {N,S,T<:SVector{N,S},D <: FixedInterval} =
FixedIntervalProduct(convert.(Ref(Domain{S}), domains))
FixedIntervalProduct(convert_eltype.(Ref(S), domains))
ProductDomain{T}(domains::SVector{N,<:FixedInterval}) where {N,T<:SVector{N}} =
FixedIntervalProduct(domains)
25 changes: 13 additions & 12 deletions src/domains/indicator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ implementing `in` directly.
abstract type AbstractIndicatorFunction{T} <: Domain{T} end

"The indicator function of a domain is the function `f(x) = x ∈ D`."
indicatorfunction(d::Domain) = x -> x ∈ d
indicatorfunction(d) = x -> x ∈ d

indomain(x, d::AbstractIndicatorFunction) = _indomain(x, d, indicatorfunction(d))
_indomain(x, d::AbstractIndicatorFunction, f) = f(x)
Expand All @@ -35,8 +35,8 @@ indicatorfunction(d::IndicatorFunction) = d.f
similardomain(d::IndicatorFunction, ::Type{T}) where {T} = IndicatorFunction{T}(d.f)

convert(::Type{IndicatorFunction}, d::AbstractIndicatorFunction) = d
convert(::Type{IndicatorFunction}, d::Domain{T}) where {T} =
IndicatorFunction{T}(indicatorfunction(d))
convert(::Type{IndicatorFunction}, d) =
IndicatorFunction{domaineltype(d)}(indicatorfunction(checkdomain(d)))

isequaldomain(d1::IndicatorFunction, d2::IndicatorFunction) =
indicatorfunction(d1)==indicatorfunction(d2)
Expand All @@ -45,18 +45,19 @@ intersectdomain1(d1::IndicatorFunction, d2) = BoundedIndicatorFunction(d1.f, d2)
intersectdomain2(d1, d2::IndicatorFunction) = BoundedIndicatorFunction(d2.f, d1)

"An indicator function with a known bounding domain."
struct BoundedIndicatorFunction{F,D,T} <: AbstractIndicatorFunction{T}
struct BoundedIndicatorFunction{T,F,D} <: AbstractIndicatorFunction{T}
f :: F
domain :: D
end

BoundedIndicatorFunction(f, domain::Domain{T}) where {T} =
BoundedIndicatorFunction{typeof(f),typeof(domain),T}(f, domain)

function BoundedIndicatorFunction(f, domain)
T = eltype(domain)
BoundedIndicatorFunction{typeof(f),typeof(domain),T}(f, domain)
end
BoundedIndicatorFunction(f, domain) =
BoundedIndicatorFunction{domaineltype(domain)}(f, domain)
BoundedIndicatorFunction{T}(f, domain::Domain{T}) where {T} =
BoundedIndicatorFunction{T,typeof(f),typeof(domain)}(f, domain)
BoundedIndicatorFunction{T}(f, domain) where {T} =
_BoundedIndicatorFunction(T, f, checkdomain(domain))
_BoundedIndicatorFunction(::Type{T}, f, domain) where {T} =
BoundedIndicatorFunction{T,typeof(f),typeof(domain)}(f, domain)

indicatorfunction(d::BoundedIndicatorFunction) = d.f

Expand All @@ -70,7 +71,7 @@ hash(d::BoundedIndicatorFunction, h::UInt) =
hashrec(indicatorfunction(d), boundingdomain(d), h)

similardomain(d::BoundedIndicatorFunction, ::Type{T}) where {T} =
BoundedIndicatorFunction(d.f, convert(Domain{T}, d.domain))
BoundedIndicatorFunction(d.f, convert_eltype(T, d.domain))

boundingbox(d::BoundedIndicatorFunction) = boundingbox(boundingdomain(d))

Expand Down
30 changes: 17 additions & 13 deletions src/domains/interval.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
iscompact(d::TypedEndpointsInterval{:closed,:closed}) = true
iscompact(d::TypedEndpointsInterval) = false

isinterval(d::Domain) = false
isinterval(d) = false
isinterval(d::AbstractInterval) = true

hash(d::AbstractInterval, h::UInt) =
Expand All @@ -23,15 +23,15 @@ approx_indomain(x, d::TypedEndpointsInterval{:open,:open}, tolerance) =
(x > leftendpoint(d)-tolerance) && (x < rightendpoint(d)+tolerance)


function point_in_domain(d::AbstractInterval)
function choice(d::AbstractInterval)
isempty(d) && throw(BoundsError())
mean(d)
end

center(d::AbstractInterval) = mean(d)

# For an interval of integers, try to find an integer point
function point_in_domain(d::AbstractInterval{T}) where {T<:Integer}
function choice(d::AbstractInterval{T}) where {T<:Integer}
isempty(d) && throw(BoundsError())
x = round(T, mean(d))
if x ∈ d
Expand All @@ -54,9 +54,9 @@ end
boundary(d::AbstractInterval) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))
corners(d::AbstractInterval) = [leftendpoint(d), rightendpoint(d)]

normal(d::AbstractInterval, x) = (abs(minimum(d)-x) < abs(maximum(d)-x)) ? -one(eltype(d)) : one(eltype(d))
normal(d::AbstractInterval, x) = (abs(minimum(d)-x) < abs(maximum(d)-x)) ? -one(domaineltype(d)) : one(domaineltype(d))

distance_to(d::AbstractInterval, x) = x ∈ d ? zero(eltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))
distance_to(d::AbstractInterval, x) = x ∈ d ? zero(domaineltype(d)) : min(abs(x-supremum(d)), abs(x-infimum(d)))

boundingbox(d::AbstractInterval) = d

Expand All @@ -79,6 +79,10 @@ and `ChebyshevInterval`.
abstract type FixedInterval{L,R,T} <: TypedEndpointsInterval{L,R,T} end
const ClosedFixedInterval{T} = FixedInterval{:closed,:closed,T}

convert(::Type{AbstractInterval{T}}, d::FixedInterval{L,R,T}) where {L,R,T} = d
convert(::Type{AbstractInterval{T}}, d::FixedInterval{L,R,S}) where {L,R,S,T} =
similardomain(d, T)

closure(d::AbstractInterval) = ClosedInterval(endpoints(d)...)
closure(d::ClosedInterval) = d

Expand Down Expand Up @@ -244,8 +248,8 @@ function similar_interval(d::HalfLine{T,C}, a::S, b::S) where {T,S,C}
HalfLine{promote_type(float(T),S),C}()
end

point_in_domain(d::NonnegativeRealLine) = zero(eltype(d))
point_in_domain(d::PositiveRealLine) = one(eltype(d))
choice(d::NonnegativeRealLine) = zero(domaineltype(d))
choice(d::PositiveRealLine) = one(domaineltype(d))


"The negative halfline `(-∞,0]` or `(-∞,0)`, right-closed or right-open."
Expand Down Expand Up @@ -276,15 +280,15 @@ function similar_interval(d::NegativeHalfLine{T,C}, a::S, b::S) where {S,T,C}
NegativeHalfLine{promote_type(S,float(T)),C}()
end

point_in_domain(d::NegativeRealLine) = -one(eltype(d))
point_in_domain(d::NonpositiveRealLine) = zero(eltype(d))
choice(d::NegativeRealLine) = -one(domaineltype(d))
choice(d::NonpositiveRealLine) = zero(domaineltype(d))


"The real line `(-∞,∞)`."
struct RealLine{T} <: FixedInterval{:open,:open,T} end
RealLine() = RealLine{Float64}()

point_in_domain(d::RealLine) = zero(eltype(d))
choice(d::RealLine) = zero(eltype(d))

similardomain(::RealLine, ::Type{T}) where {T} = RealLine{T}()

Expand Down Expand Up @@ -318,11 +322,11 @@ intersect(d1::FixedInterval, d2::FixedInterval) = intersectdomain(d1, d2)

# Promotion to joint type T
uniondomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) =
uniondomain(promote(d1,d2)...)
uniondomain(promote_domains(d1,d2)...)
intersectdomain(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) =
intersectdomain(promote(d1,d2)...)
intersectdomain(promote_domains(d1,d2)...)
setdiffdomain(d1::AbstractInterval, d2::AbstractInterval) =
setdiffdomain(promote(d1,d2)...)
setdiffdomain(promote_domains(d1,d2)...)



Expand Down
2 changes: 1 addition & 1 deletion src/domains/point.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ supremum(d::Point) = pointval(d)
interior(d::Point) = emptyspace(d)
closure(d::Point) = d

point_in_domain(d::Point) = pointval(d)
choice(d::Point) = pointval(d)

distance_to(d::Point, x) = norm(x-pointval(d))

Expand Down
2 changes: 1 addition & 1 deletion src/domains/simplex.jl
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ closure(d::VectorUnitSimplex{T}) where {T} = VectorUnitSimplex{T,:closed}(dimens

# We pick the center point, because it belongs to the domain regardless of
# whether it is open or closed.
point_in_domain(d::UnitSimplex) = center(d)
choice(d::UnitSimplex) = center(d)

similardomain(d::StaticUnitSimplex{S,C}, ::Type{T}) where {S,T,C} =
StaticUnitSimplex{T,C}()
Expand Down
13 changes: 8 additions & 5 deletions src/domains/trivial.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ EmptySpace(::Type{T}) where {T} = EmptySpace{T}()
similardomain(::EmptySpace, ::Type{T}) where {T} = EmptySpace{T}()

"Return the empty space with the same element type as the given domain."
emptyspace(d) = emptyspace(eltype(d))
emptyspace(d) = emptyspace(domaineltype(d))
emptyspace(::Type{T}) where {T} = EmptySpace{T}()

indomain(x::T, d::EmptySpace{T}) where {T} = false
Expand Down Expand Up @@ -64,14 +64,17 @@ struct FullSpace{T} <: Domain{T} end
const AnyFullSpace = FullSpace{Any}

FullSpace() = FullSpace{Float64}()
FullSpace(d) = FullSpace{eltype(d)}()
FullSpace(d) = FullSpace{domaineltype(d)}()

"Return the full space with the same element type as the given domain."
fullspace(d) = fullspace(eltype(d))
fullspace(d) = fullspace(domaineltype(d))
fullspace(::Type{T}) where {T} = FullSpace{T}()

isfullspace(d::FullSpace) = true
isfullspace(d::Domain) = false
isfullspace(d) = _isfullspace(d, DomainStyle(d))
_isfullspace(d, ::IsDomain) = false
_isfullspace(d, ::NotDomain) = error("isfullspace invoked on non-domain type")

similardomain(::FullSpace, ::Type{T}) where {T} = FullSpace{T}()

Expand All @@ -84,7 +87,7 @@ approx_indomain(x, d::FullSpace, tolerance) = in(x, d)
show(io::IO, d::FullSpace) = print(io, "{x} (full space)")

# We choose the origin as a point in the full space
point_in_domain(d::FullSpace) = zero(eltype(d))
choice(d::FullSpace) = zero(eltype(d))

isempty(::FullSpace) = false

Expand Down Expand Up @@ -125,7 +128,7 @@ type `T`.
struct TypeDomain{T} <: Domain{T} end

"Return the domain for the element type of the given domain."
typedomain(d) = typedomain(eltype(d))
typedomain(d) = typedomain(domaineltype(d))
typedomain(::Type{T}) where {T} = TypeDomain{T}()

iscompatiblepair(x::T, d::TypeDomain{T}) where {T} = true
Expand Down
Loading