juliaapproximation / domainsets.jl Goto Github PK

The way a domain is constructed from a set is currently inconsistent:

julia> s = Set([ [2,3], [4,5]])
Set(Array{Int64,1}[[2, 3], [4, 5]])

julia> d = Domain(s)
a union of 2 domains:
	1.	: DomainSets.WrappedDomain{Array{Int64,1},Int64}([2, 3])
	2.	: DomainSets.WrappedDomain{Array{Int64,1},Int64}([4, 5])


julia> 2 ∈ d
true

julia> 2 ∈ s
false

julia> [2,3] ∈ s
true

julia> [2,3] ∈ d
ERROR: MethodError: Cannot `convert` an object of type Array{Int64,1} to an object of type Int64

It seems safer to wrap a set into a wrapped domain as a whole, rather than to make a union of wrapped domains out of it. An alternative is to wrap each element of a set into a Point domain.

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

Since × means cross which means the cross product, I'm not sure it's good style to overload to also mean cartesianproduct. In particular, in my case I'm using AbstractVectors as domains and wanted to form the Cartesian product, and was surprised by the error.

eltype should be reserved for the iterator interface: otherwise you get the following inconsistency:

julia> Set(1)
Set([1])

julia> Set(interval(0,1))
ERROR: MethodError: no method matching start(::Domains.Interval{:closed,:closed,Float64})
Closest candidates are:
  start(::SimpleVector) at essentials.jl:258
  start(::Base.MethodList) at reflection.jl:560
  start(::ExponentialBackOff) at error.jl:107
  ...
Stacktrace:
 [1] union!(::Set{Float64}, ::Domains.Interval{:closed,:closed,Float64}) at ./set.jl:63
 [2] Set(::Domains.Interval{:closed,:closed,Float64}) at ./set.jl:10

julia> 5 in ℝ
false

Also, its weird that ℝ is a type, not an instance

I would expect these both to be SVector{2,Int}:

julia> DomainSets.ProductDomain(0..1, 2..3) |> eltype
StaticArrays.SArray{Tuple{2},Int64,1,2}

julia> DomainSets.ProductDomain(0..1, [2,3]) |> eltype # Why Tuple?
Tuple{Int64,Int64}

People (other than me) don't like non-ASCII, so let's avoid the complaints and make ∂ shorthand for boundary, not the other way around.

A Domain has an eltype, but it is not iterable. This leads to some problems when using functionality that is defined with sets in mind.

Here is a simple problem:

julia> using DomainSets

julia> Set(0..1)
ERROR: MethodError: Cannot `convert` an object of type Interval{:closed,:closed,Int64} to an object of type Int64
Closest candidates are:
  convert(::Type{T<:Number}, ::T<:Number) where T<:Number at number.jl:6
  convert(::Type{T<:Number}, ::Number) where T<:Number at number.jl:7
  convert(::Type{T<:Integer}, ::Ptr) where T<:Integer at pointer.jl:23

You can not create a Set that contains a Domain this way. The way to do it is:

julia> Set([0..1])
Set(Interval{:closed,:closed,Int64}[0..1])

That is fair enough, but it leads to other issues. For example, we can not compute the intersection of three domains:

julia> using DomainSets

julia> I = 0..1
0..1

julia> intersect(I, I)
0..1

julia> intersect(I, I, I)
ERROR: MethodError: Cannot `convert` an object of type Interval{:closed,:closed,Int64} to an object of type Int64
Closest candidates are:
  convert(::Type{T<:Number}, ::T<:Number) where T<:Number at number.jl:6
  convert(::Type{T<:Number}, ::Number) where T<:Number at number.jl:7
  convert(::Type{T<:Integer}, ::Ptr) where T<:Integer at pointer.jl:23
  ...
Stacktrace:
 [1] setindex!(::Dict{Int64,Nothing}, ::Nothing, ::Interval{:closed,:closed,Int64}) at ./dict.jl:373
 [2] push!(::Set{Int64}, ::Interval{:closed,:closed,Int64}) at ./set.jl:48
 [3] union!(::Set{Int64}, ::Interval{:closed,:closed,Int64}) at ./abstractset.jl:81
 [4] Set{Int64}(::Interval{:closed,:closed,Int64}) at ./set.jl:10
 [5] _Set(::Interval{:closed,:closed,Int64}, ::Base.HasEltype) at ./set.jl:23
 [6] Set(::Interval{:closed,:closed,Int64}) at ./set.jl:21
 [7] _shrink(::Function, ::Interval{:closed,:closed,Int64}, ::Tuple{Interval{:closed,:closed,Int64},Interval{:closed,:closed,Int64}}) at ./array.jl:2396
 [8] intersect(::Interval{:closed,:closed,Int64}, ::Interval{:closed,:closed,Int64}, ::Interval{:closed,:closed,Int64}) at ./array.jl:2400
 [9] top-level scope at none:0

The two-argument function succeeds because we defined it ourselves. The three-argument version fails because the fallback in Julia tries to make a Set out of the domains.

Also, we currently have this inconsistent behaviour of union:

julia> I = 0..1
0..1

julia> union(I)
a union of 1 domains:
        1.      : 0..1


julia> union(I,I)
0..1

julia> union(I,I,I)
a union of 3 domains:
        1.      : 0..1
        2.      : 0..1
        3.      : 0..1

Of course we could fix intersect and union ourselves, but perhaps there is a better change that would make Domain integrate better with the existing Julia functions.

I'm thinking instead of the chebyshev interval [-1,1], maybe -1 .. 1 (Chebyshev)

@JuliaRegistrator register()

• Get ApproxFun to compile using Domains.jl
• Move any "general purpose" domains over (Arc, UnionDomain)
• Delete redundant domains in ApproxFun
• #1 Support complex domains
• #2
• #3 Support derivatives of maps
• #4 Publish to METADATA
• #6

I'm considering moving ApproxFun over now: the only scalable way to do JuliaApproximation/ApproxFun.jl#531 I can see is by adding support for Chebyhev(ChebyshevInterval()), which returns points ChebyshevGrid(n), while Chebyshev(Interval(a,b)) has points MappedGrid(ChebyshevGrid(n), AffineMap(...)). This change to ChebyshevInterval would also simplify a lot of the code.

The notion of a product of domains is really about concatenation of the inputs. Renaming it to VcatDomain makes this clear, and opens up the possibility of HcatDomain and HvcatDomain, so that:

[1; 2] in VcatDomain(1..2, 2..3)
[1 2] in HcatDomain(1..2, 2..3)
[1 2; 3 4] in HvcatDomain([1..2 1...2; 3..4 3..4])

Hi, this package seems like it could be really useful for Bayesian modeling, specifically the kind of transformations Stan implements. It would be really helpful to have some simple examples as a jumping-off point. I'd like to try to implement things like logistic and log transforms, but I need to get my head around the code first. Anything you can point me to?

I've done this before so I can take care of it

It's better to use Set since this avoids redundant domains, and also makes == easy to implement.

The concept of isomorphisms is useful, but is orthogonal to the concept of domains. The way the package interweaves these two ideas makes the code hard to modify. I'd propose removing all references to isomorphisms from DomainSets.jl and moving them to another package, say ArrayIsomorphisms.jl. (Isomorphisms.jl seems to broad.)

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

There are too many convenience functions, which goes against a good practice code design principle that there should only be one way to write anything (to make find-and-replace easier, etc.). I think all the lowercase convenience functions can be removed: UnitDisk is just as good as disk, etc.

Defining const Point{N,T} = SVector{N,T} prevents the useful definition

struct Point{N,T} <: Domain{T}
  x::T
end

:open, :closed should apply to any interval, so I don't see why this wouldn't be part of AbstractInterval.

I am completely new to Julia. I have started it because I have heard some pretty good things about the package ApproxFun. I am eager to try but unfortunately so far I cannot make it run on my laptop.

I am using Julia 1.0.2
My computer has 8,00 RAM
My operating system is Windows 7

I have no problem adding Approxfun:

input:
import Pkg;
Pkg.add(“ApproxFun”)
Pkg.status()

output:
Resolving package versions…
Updating C:\Users\charlierch.julia\environments\v1.0\Project.toml
[28f2ccd6] + ApproxFun v0.11.1
Updating C:\Users\charlierch.julia\environments\v1.0\Manifest.toml
[621f4979] + AbstractFFTs v0.4.1
[28f2ccd6] + ApproxFun v0.11.1
[fbd15aa5] + ApproxFunBase v0.1.3
[59844689] + ApproxFunFourier v0.1.2
[b70543e2] + ApproxFunOrthogonalPolynomials v0.2.0
[f8fcb915] + ApproxFunSingularities v0.1.0
[aae01518] + BandedMatrices v0.9.3
[8e7c35d0] + BlockArrays v0.9.1
[ffab5731] + BlockBandedMatrices v0.4.4
[00ebfdb7] + CSTParser v0.6.1
[49dc2e85] + Calculus v0.5.0
[717857b8] + DSP v0.5.2
[5b8099bc] + DomainSets v0.1.0
[fa6b7ba4] + DualNumbers v0.6.2
[7a1cc6ca] + FFTW v0.2.4
[442a2c76] + FastGaussQuadrature v0.3.3
[057dd010] + FastTransforms v0.5.0
[1a297f60] + FillArrays v0.6.3
[7c893195] + HierarchicalMatrices v0.1.4
[4858937d] + InfiniteArrays v0.1.1
[8197267c] + IntervalSets v0.3.1
[5078a376] + LazyArrays v0.9.0
[898213cb] + LowRankApprox v0.2.3
[1914dd2f] + MacroTools v0.5.1
[a3b82374] + MatrixFactorizations v0.0.4
[4d1e1d77] + Nullables v0.0.8
[f27b6e38] + Polynomials v0.5.2
[92933f4c] + ProgressMeter v1.0.0
[c751599d] + ToeplitzMatrices v0.5.0
[0796e94c] + Tokenize v0.5.5
[9abbd945] + Profile
Status C:\Users\charlierch.julia\environments\v1.0\Project.toml
[28f2ccd6] ApproxFun v0.11.1
[7073ff75] IJulia v1.14.1
[91a5bcdd] Plots v0.26.0
[276daf66] SpecialFunctions v0.7.2
[37e2e46d] LinearAlgebra

However I have a problem with “using ApproxFun”

input:
using ApproxFun

output:
┌ Info: Precompiling ApproxFun [28f2ccd6-bb30-5033-b560-165f7b14dc2f]
└ @ Base loading.jl:1192
ERROR: LoadError: LoadError: StackOverflowError:
in expression starting at C:\Users\charlierch.julia\packages\DomainSets\D8lcW\src\spaces\productspace.jl:107
in expression starting at C:\Users\charlierch.julia\packages\DomainSets\D8lcW\src\DomainSets.jl:152
ERROR: LoadError: Failed to precompile DomainSets [5b8099bc-c8ec-5219-889f-1d9e522a28bf] to C:\Users\charlierch.julia\compiled\v1.0\DomainSets\ksPnB.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:311 [inlined]
[3] compilecache(::Base.PkgId, ::String) at .\loading.jl:1190
[4] macro expansion at .\logging.jl:309 [inlined]
[5] _require(::Base.PkgId) at .\loading.jl:947
[6] require(::Base.PkgId) at .\loading.jl:858
[7] macro expansion at .\logging.jl:309 [inlined]
[8] require(::Module, ::Symbol) at .\loading.jl:840
[9] include at .\boot.jl:317 [inlined]
[10] include_relative(::Module, ::String) at .\loading.jl:1044
[11] include(::Module, ::String) at .\sysimg.jl:29
[12] top-level scope at none:2
[13] eval at .\boot.jl:319 [inlined]
[14] eval(::Expr) at .\client.jl:393
[15] top-level scope at .\none:3
in expression starting at C:\Users\charlierch.julia\packages\ApproxFun\f2ab0\src\ApproxFun.jl:4

Failed to precompile ApproxFun [28f2ccd6-bb30-5033-b560-165f7b14dc2f] to C:\Users\charlierch.julia\compiled\v1.0\ApproxFun\jGqLz.ji.

Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:311 [inlined]
[3] compilecache(::Base.PkgId, ::String) at .\loading.jl:1190
[4] macro expansion at .\logging.jl:309 [inlined]
[5] _require(::Base.PkgId) at .\loading.jl:947
[6] require(::Base.PkgId) at .\loading.jl:858
[7] macro expansion at .\logging.jl:309 [inlined]
[8] require(::Module, ::Symbol) at .\loading.jl:840
[9] top-level scope at In[37]:1

It seems the StackOverflowError is related to DomainSets. I have no problem with the installation of this package separately:

input:
Pkg.add(“DomainSets”)
Pkg.status()

output:
Resolving package versions…
Updating C:\Users\charlierch.julia\environments\v1.0\Project.toml
[5b8099bc] + DomainSets v0.1.0
Updating C:\Users\charlierch.julia\environments\v1.0\Manifest.toml
[no changes]
Status C:\Users\charlierch.julia\environments\v1.0\Project.toml
[28f2ccd6] ApproxFun v0.11.1
[5b8099bc] DomainSets v0.1.0
[7073ff75] IJulia v1.14.1
[91a5bcdd] Plots v0.26.0
[276daf66] SpecialFunctions v0.7.2
[37e2e46d] LinearAlgebra

But again I face the problem with “using DomainSets”

input:
using DomainSets

output:
┌ Info: Precompiling DomainSets [5b8099bc-c8ec-5219-889f-1d9e522a28bf]
└ @ Base loading.jl:1192
ERROR: LoadError: LoadError: StackOverflowError:
in expression starting at C:\Users\charlierch.julia\packages\DomainSets\D8lcW\src\spaces\productspace.jl:107
in expression starting at C:\Users\charlierch.julia\packages\DomainSets\D8lcW\src\DomainSets.jl:152

Failed to precompile DomainSets [5b8099bc-c8ec-5219-889f-1d9e522a28bf] to C:\Users\charlierch.julia\compiled\v1.0\DomainSets\ksPnB.ji.

Stacktrace:
[1] error(::String) at .\error.jl:33
[2] macro expansion at .\logging.jl:311 [inlined]
[3] compilecache(::Base.PkgId, ::String) at .\loading.jl:1190
[4] macro expansion at .\logging.jl:309 [inlined]
[5] _require(::Base.PkgId) at .\loading.jl:947
[6] require(::Base.PkgId) at .\loading.jl:858
[7] macro expansion at .\logging.jl:309 [inlined]
[8] require(::Module, ::Symbol) at .\loading.jl:840
[9] top-level scope at In[39]:1

In the error messages I have written above, there is
ERROR: LoadError: Failed to precompile DomainSets [5b8099bc-c8ec-5219-889f-1d9e522a28bf] to C:\Users\charlierch.julia\compiled\v1.0\DomainSets\ksPnB.ji.

But when I go in the folder Users\charlierch.julia\compiled\v1.0\DomainSets, there is no file at all. I guess I should have ksPnB.ji but there is nothing.

I have tried several times Pkg.remove("ApproxFun"), Pkg.add("DomainSets"), and then using DomainSets, but the same problem persists.

Many thanks in advance for your help !

The problem is the use of broadcast in setdiff for a UnionDomain and a Domain:

julia> UnionDomain(Interval(0.0, 0.5), Interval(1.0,1.5)) \ Interval(0.0, 1.5)
ERROR: MethodError: no method matching iterate(::Interval{:closed,:closed,Float64})
Closest candidates are:
  iterate(::Core.SimpleVector) at essentials.jl:589
  iterate(::Core.SimpleVector, ::Any) at essentials.jl:589
  iterate(::ExponentialBackOff) at error.jl:171
  ...

The offending line is this one: broadcast is not appropriate here.

I would like to use the new naming conventions that are in master.

This is the evolution of the ideas that have been floating around about domains being an interface, not a type.

I think to support domains from other packages and Base (like IntervalSets, or Set) we should re-structure this package a bit to use traits. That is, we'd have the following:

abstract type DomainStyle end
struct IntervalStyle <: DomainStyle end
struct PointStyle <: DomainStyle end
struct SetStyle{T} <: DomainStyle end

DomainStyle(::Type{<:IntervalSets.Interval}) = IntervalStyle()
DomainStyle(::Type{<:Number}) = PointStyle()
DomainStyle(::Type{<:AbstractArray{T}}) where T = SetStyle{T}()
DomainStyle(::Type{<:AbstractSet{T}}) where T = SetStyle{T}()

I think the best way to implement this proposal is to do this in a separate branch, side-by-side with ApproxFun's adoption of Domains.jl. The result will be that, for example, we can do f = Fun(exp, [1,2,3]) for a function supported on the points [1,2,3] instead of g = Fun(exp, union(Point(1), Point(2), Point(3)), and have the property that domain(f) === [1,2,3]. And also domain(Fun(exp, 1..2)) === 1..2.

This will ease the "wrapper bloat" of trying to make everything a subtype Domain.

Julia v0.7 has much improved support for broadcasting, for example, one can now more easily capture complex function calls like exp.( log.(x .+ 1) .+ 2)). We can exploit this for domains: if d is a domain then ``exp.( log.(d .+ 1) .+ 2))` should be the mapped domain, and we can capture each function call and turn it into a composition of maps.

Here's a list of packages that overlap with this and might influence the design:

https://github.com/JuliaGeometry/GeometryTypes.jl
https://github.com/zenna/AbstractDomains.jl
https://github.com/tpapp/ContinuousTransformations.jl/blob/master/README.md
https://github.com/JuliaReach/LazySets.jl

I like overriding ctranspose so that

m = LinearMap(5)
m' # returns `ConstantMap` or something

Would it be possible for VectorProductDomain to have a type parameter for the array type, to allow allow for using a CuArray, StructArray, etc.? If so, adding support for Adapt.jl might make sense in the same context.

One would most likely to expect that arithmetic with intervals is interval arithmetic, not union.

Bad:

julia> ChebyshevInterval()
the interval [-1.0, 1.0]

julia> Interval(-1,1.)
the interval [-1.0, 1.0]

For higher-dimensional problems it would be great to have support for non-static dimensionality and argument types like Vector in addition to SVector. E.g. a ProductDomain based on AbstractVector{IntervalSet} should allow for in(x::AbstractVector, d::ProductDomain). N-dimensional hyper-spheres with N not fixed at compile time would also come on very handy.

Is this actually needed in practice? I think overriding broadcast is less confusing as domains don't really form an "algebra".

This seems wrong to me, it should be an SVector{3}:

julia> using DomainSets, StaticArrays

julia> struct Disk <: Domain{SVector{2,Float64}} end

julia> d = (0.0..1) × Disk();

julia> eltype(d)
Tuple{Float64,SArray{Tuple{2},Float64,1,2}}

We are converging to Domain as an interface: it is any object type that supports eltype and in. We make methods with which one can combine domains, where the individual domains don't have to be Domain types. This is clearly very flexible. Perhaps we can do the same for maps and make the framework more generic/applicable?

A map is really just a function: y = f(x). Thus, perhaps any type that supports the calling operator can be seen as a map. Or perhaps any type that supports * can be a map. We could treat a matrix A as a map, the map then being the matrix-vector product y = A*x. Even a number could be a map: 2 represents multiplication by 2. As long as * is associative (which I think Julia already assumes in many cases), the interpretation of * as applying a map is consistent with composition: (A*B)*x = A*(B*x). With this approach, I think we would not even have to implement affine maps. We can just use the existing linear algebra functions and vectors/matrices.

Julia functions have a known output type if you give it a (combination of) input type(s). The current definition of AbstractMap{S,T} really is like a typed function. It corresponds to a function for which you restrict the input type to be S. The output type T is perhaps superfluous. From this point of view, a better name might be TypedMap{S} or TypedFunction{S}. It could have an implementation like:

(*)(f::TypedFunction{S}, x) where {S} = *(f, convert(S, x))

Any concrete subtype would implement multiplication for arguments of type S only.

In order to generically support maps like 2, cos and rand(2,2), it is somewhat restrictive to assume that the input type is known. We can in general assume that, given the map and the input type, the output type is known. (Worst case it is Any.) Thus, we could have TypedMap{S} <: AbstractMap where AbstractMap has no type parameters at all. What the map does, depends on the type of the input and is expressed by implementing (*)(m::MyMap, x).

Ideally, eventually, computing Jacobians etcetera is done generically via automatic differentiation. Composition of maps ideally would use fusion (like broadcast), rather than layered CompositeMap's or CompositeMap's with many elements.

As with domains, as long as you are manipulating variables of type AbstractMap we can use standard function names. If you intend to manipulate any object that implements the map interface, whatever we decide it to be, the intention has to be implicit in the name of the function that is being called (like convert_map and promote_map) or the operation has to be consistent with the generic meaning of the function name (as with *).

I'm sure there'll be some problems I'm not seeing at the moment too.

I don't have an immediate problem with it, but currently we have:

julia> (-1..1) \ (0..1)
-1..0

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

The result of the difference of two closed intervals should probably be a half-open one, i.e., the correct result of [-1,1] \ [0,1] arguably should be:

julia> Interval{:closed,:open,Float64}(-1, 0)
-1.0..0.0 (closed–open)

That would mean changing the setdiff function in interval.jl here.

In the package IntervalSets, isopen and isclosed will soon be renamed to isopenset and isclosedset (see JuliaMath/IntervalSets.jl#59). We will have to follow suit.

One issue I noticed is that the current definition of the boundary of an interval is unexpected. It does not include endpoints if the interval is open:

julia> boundary(ClosedInterval(0.,1.))
a union of 2 domains:
	1.	: Point{Float64}(0.0)
	2.	: Point{Float64}(1.0)


julia> boundary(OpenInterval(0.,1.))
the empty space with eltype Float64

If an interval is seen as a subset of the real line, then its boundary includes both endpoints, regardless of whether the interval is open or not. I.e., I would expect:

boundary(d::AbstractInterval) = Point(leftendpoint(d)) ∪ Point(rightendpoint(d))

Then, one could define the closure of an open interval as the corresponding closed interval, and the interior of the closed interval as the open one - the boundary would be the difference between the two.

Is there a reason for the current implementation? Now seems like a good time to change it.

I think the current implementation of in adds more complexity than it removes:

function in(x, d::Domain)
   T = promote_type(typeof(x), eltype(d))
   T == Any && return false
   in(convert(T, x), convert(Domain{T}, d))
end

I think it would be simpler if each domain just overloads in, and can handle types however it sees fit...

There is no rush, we could wait to see whether the package will ultimately prove to be useful, but it certainly doesn't have to stay here long term.

It's not clear why UnitSphere refers to an N-sphere while Sphere is a 3-sphere. I'd expect a Sphere to also have a radius and centre (down the line). So maybe UnitSphere -> UnitHyperSphere while Sphere -> UnitSphere.

Before this happens, there should be a discussion with https://github.com/zenna/AbstractDomains.jl about unifying, to avoid multiple packages that do the same thing.

I'm not sure what the conclusion was what to call it, but there should be some version of arclength for an interval.

It could be called lebesguemeasure, though that's a bit confusing for intervals embedded in the complex plane, as the measure with respect to dx*dy would be zero....

@JuliaRegistrator register()

“Domains” is a little bit ambiguous as it is, ahem, a domain specific term. I suspect when it comes to tagging a release there’ll be some quibbles over the name.

I like “InfiniteSets.jl” as it accurately reflects that this is a generalisation of Set .

Also, I like the rhyming with InfiniteArrays.jl.

This is currently a bit unwieldy ... :-)

ProductDomain([0..i for i in 1:500])

Should probably use a ... or so to shorten it, like show for long arrays.