support domains defined by a characteristic function

[daanhb]

A domain could be defined by a user-supplied function, as in(x, domain) = d.f(x).

[daanhb]

There is a cool unexpected aspect to indicator functions, they allow for a syntax based on generators. Like this:

julia> d = Domain(x > 0 for x in -1..1)
indicator function bounded by: -1..1

The generator gives information about the bounding box (here it is -1..1) and it also gives the characteristic function, which should evaluate to true or false. So the domain above is the part of the interval -1..1 on which x is positive. We don't know explicitly that the result is the (half-open) interval 0..1, at least I wouldn't know how to produce that knowledge automatically. But still, this is a neat new way to define a domain.

It also works with higher-dimensional domains and product domains, or combinations of them:

julia> d = Domain( x+y+z > 0 for (x,y) in UnitDisk(), z in 0..1)
indicator function bounded by: the 2-dimensional closed unit ball x 0..1

julia> ([0.3,0.2], 0.5) ∈ d
true

This domain would be all points of the form ([x,y],z) for which [x,y] lies in the unit disk, z lies in the unit interval, and the sum of all three elements is positive. It is the part of a cylinder that lies above a hyperplane.

The repeated for expressions give rise to nested tuples, which is a bit awkward to process afterwards. It is easier like this:

julia> d = Domain( x+y+z > 0 for (x,y,z) in productdomain(UnitDisk(), 0..1))
indicator function bounded by: the 2-dimensional closed unit ball x 0..1

julia> [0.3,0.2, 0.5] ∈ d
true

All this is on the branch #78.

[dlfivefifty]

That is amazing! Can we use that also for integration? integrate(exp(x) for x in 0..1) ?

[daanhb]

Yes we can and I'm working on it! See DomainIntegrals.jl. Example (on master):

julia> integral(exp(x+y) for (x,y) in (0..1)^2)
2.9524924420120535

I'm trying to teach DomainIntegrals to evaluate integrals on various domains, using maps from intervals and squares, so I'll be working on improving maps and parameterizations.

[daanhb]

Update, because it gets even better: the generator syntax actually supplies information about the bounding box and this is useful. An indicator function that knows its bounding box is created as follows:

julia> I = Domain( x+y <= 1 for (x,y) in (0..1) × (0..1) )

julia> boundingbox(I)
(0..1) × (0..1)

This is the unit simplex in 2D, bounded by [0,1]^2. With that information, DomainIntegrals can evaluate integrals on this domain:

julia> using DomainIntegrals

julia> integral(x->1, Domain( x+y <= 1 for (x,y) in (0..1) × (0..1) ) )
0.5000123377344884

julia> integral(cos(x[1]+x[2]) for x in Domain( x+y <= 1 for (x,y) in (0..1) × (0..1) ))
0.3819126664757159

It works by extending the integrand by zero to the bounding box and to evaluate on the unit box, because adaptive integration methods can deal with the discontinuity. At the cost of performance though, the values above are not highly accurate unless the number of allowed function evaluations is increased.

julia> integral(x->1, UnitDisk())
3.1396125320956303

julia> abs(ans-pi)
0.001980121494162823