Implement rand()

[zsunberg]

Fix #109

[codecov]

Codecov Report

Merging #112 (482ab5f) into master (60bb050) will increase coverage by 0.21%.
The diff coverage is 94.73%.

@@            Coverage Diff             @@
##           master     #112      +/-   ##
==========================================
+ Coverage   85.80%   86.01%   +0.21%     
==========================================
  Files          30       31       +1     
  Lines        2479     2496      +17     
==========================================
+ Hits         2127     2147      +20     
+ Misses        352      349       -3     

Impacted Files	Coverage Δ
src/applications/random.jl	94.11% <94.11%> (ø)
src/domains/boundingbox.jl	100.00% <100.00%> (ø)
src/domains/cube.jl	96.50% <100.00%> (ø)
src/generic/mapped.jl	83.11% <0.00%> (+5.19%)	⬆️

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[zsunberg]

It took me like 3 hours to understand the type system, but I think this whole implementation is going to be < 50 lines when done. 😂

[daanhb]

Yes, that sounds very Julian: the ratio of time it takes to write code versus length of code can be pretty high (especially with lots of types flying around working together in undocumented ways as is the case here...sorry about that).

On the one hand, one expects there to be an implementation of rand for some or all of the "basic" types. And then it should all "just work" together. For inspiration on that second part you can look at how point_in_domain works. For example, product domains are dealt with using this line:

point_in_domain(d::CompositeDomain) = toexternalpoint(d, map(point_in_domain, components(d)))

The map command will apply point_in_domain to all components of the composite domain. (It is called factors only for product domains, components is slightly more general). The "toexternalpoint" machinery will translate the internal representation of the element to the external one (of type T). That could mean that a tuple of a Float64 and an SVector{2,Float64} is concatenated into an SVector{3}, in case of a product domain of an interval with a disk, say. That is something map would not do on its own. It works like this:

julia> using DomainSets

julia> using DomainSets: ×

julia> d = (1..2.0) × UnitDisk()
(1.0..2.0) × UnitDisk()

julia> eltype(d)
SVector{3, Float64} (alias for StaticArrays.SArray{Tuple{3}, Float64, 1, 3})

julia> x = point_in_domain(d)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 1.5
 0.0
 0.0

julia> DomainSets.tointernalpoint(d, x)
(1.5, [0.0, 0.0])

julia> typeof(ans)
Tuple{Float64, StaticArrays.SVector{2, Float64}}

julia> map(DomainSets.point_in_domain, components(d))
(1.5, [0.0, 0.0])

julia> DomainSets.toexternalpoint(d, ans)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 1.5
 0.0
 0.0

The internal representation works will with map, the external representation doesn't but makes more sense for the user.

Probably rand only makes sense for product domains. It would be difficult to define properly for unions, intersections and differences of domains (the other composite domains), because there is no guarantee about whether or not the composing domains have any overlap - the distribution could be skewed because of that. Same for mapped domains: some kind of random point could be defined but, unless the map is affine, what is its distribution?

[daanhb]

I think here you'd want to put a toexternalpoint(d, outcome_of_map) around it, where d is the product domain

[zsunberg]

@daanhb I implemented the easy ones. I think this is a good enough minimal implementation to be merged, but I would be willing to take a stab at it if you think others should be included.

[daanhb]

@daanhb I implemented the easy ones. I think this is a good enough minimal implementation to be merged, but I would be willing to take a stab at it if you think others should be included.

Just the easy cases is fine! If it is isn't easy, then perhaps it should not be in this package after all.

[daanhb]

boundingbox(d) allocates memory and is not very optimized, perhaps it is better to define box=boundingbox(b) once before the while loop?

[zsunberg]

boundingbox(d) allocates memory and is not very optimized

allocates memory even with StaticArrays? that is unfortunate.

[daanhb]

This looks good. Could you try to add my suggestions in the issues you opened and see if that fixes the broken tests?

[zsunberg]

Ok, everything is fixed - I think it's ready to go!

[daanhb]

I just noticed that

julia> using DomainSets

julia> rand(ChebyshevInterval())
ERROR: ArgumentError: Sampler for this object is not defined

which also causes

julia> rand(UnitBall())
ERROR: ArgumentError: Sampler for this object is not defined

The reason is that rand in IntervalSets is defined for ClosedInterval. I think it could be defined for the more general <:TypedEndpointsInterval{:closed,:closed,T}, and then rand would also work for all of our intervals with fixed endpoints.

[zsunberg]

FYI @daanhb I fixed the issue you mentioned above at JuliaMath/IntervalSets.jl#116