Core

Project.toml

@@ -1,10 +1,11 @@
 name = "DomainIntegrals"
 uuid = "cc6bae93-f070-4015-88fd-838f9505a86c"
 authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
-version = "0.4.4"
+version = "0.4.5"
 
 [deps]
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"
+DomainSetsCore = "b5e7cfa8-5ebe-46e7-951c-e7d99cb94c6d"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 GaussQuadrature = "d54b0c1a-921d-58e0-8e36-89d8069c0969"
@@ -15,15 +16,15 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+julia = "1.6"
 CompositeTypes = "0.1.3"
-DomainSets = "0.6.2"
+DomainSets = "0.6.7,0.7"
 FastGaussQuadrature = "0.5"
 GaussQuadrature = "0.5.7"
 HCubature = "1.4.0"
 IntervalSets = "0.7"
 QuadGK = "2.3.1"
 StaticArrays = "1.0"
-julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/DomainIntegrals.jl

@@ -21,6 +21,7 @@ import DomainSets:
 
 using DomainSets: ×,
     EmptySpace, FullSpace,
+    domain_prectype, domain_numtype,
     euclideandimension
 
 

src/integral.jl

@@ -8,7 +8,7 @@ The suggested quadrature strategy based on the `integral` arguments.
 
 By default, adaptive quadrature is chosen.
 """
-suggestedstrategy(domain::Domain, args...) = QuadAdaptive{prectype(domain)}()
+suggestedstrategy(domain, args...) = QuadAdaptive{prectype(AsDomain(domain))}()
 
 returntype(integrand, ::Type{S}) where {S} = Base.Core.Compiler.return_type(integrand, (S,))
 
@@ -46,15 +46,7 @@ associated_domain(μ::Measure) = support(μ)
 associated_domain(μ::DiscreteWeight{T}) where {T} = DummyDiscreteDomain{T}()
 
 # associate a measure with a domain
-associated_measure(domain::Domain{T}) where {T} = Lebesgue{T}()
-
-# Process the arguments until there is a domain, a measure and zero or more property objects.
-process_arguments(measure::Measure, properties::Property...) =
-    process_arguments(associated_domain(measure), measure, properties...)
-process_arguments(domain::Domain, properties::Property...) =
-    process_arguments(domain, associated_measure(domain), properties...)
-process_arguments(domain::Domain, measure::LebesgueDomain{T}, properties::Property...) where {T} =
-    process_arguments(domain ∩ support(measure), Lebesgue{T}(), properties...)
+associated_measure(domain) = Lebesgue{domaineltype(domain)}()
 
 promote_domain_and_measure(domain::Domain{T}, measure::Measure{T}) where {T} =
     domain, measure
@@ -63,14 +55,22 @@ function promote_domain_and_measure(domain::Domain{S}, measure::Measure{T}) wher
     convert(Domain{U}, domain), convert(Measure{U}, measure)
 end
 
+
+# Process the arguments until there is a domain (not necessarily of Domain type),
+# a measure and zero or more property objects.
+process_arguments(measure::Measure, properties::Property...) =
+    process_arguments(associated_domain(measure), measure, properties...)
+process_arguments(domain, properties::Property...) =
+    process_arguments(checkdomain(domain), associated_measure(domain), properties...)
+
 function process_arguments(domain::Domain{S}, measure::Measure{T}, properties::Property...) where {S,T}
     d, m = promote_domain_and_measure(domain, measure)
     process_arguments(d, m, properties...)
 end
-
-# all good now
 process_arguments(domain::Domain{T}, measure::Measure{T}, properties::Property...) where {T} =
     (domain, measure, properties...)
+process_arguments(domain, measure::Measure{T}, properties::Property...) where {T} =
+    (domain, measure, properties...)
 
 process_arguments(args...) = error("Arguments to integral or integrate functions not understood.")
 
@@ -79,31 +79,31 @@ integrate(qs::ChebyshevIntervalRule{T}, integrand) where {T} =
 integrate(qs::HalfLineRule{T}, integrand) where {T} =
     integrate(qs, integrand, HalfLine{T}())
 integrate(qs::RealLineRule{T}, integrand) where {T} =
-    integrate(qs, integrand, FullSpace{T}())
+    integrate(qs, integrand, RealLine{T}())
 
 ## The use of generators
 integrate(gen::Base.Generator, args...) =
     integrate(process_generator(gen)..., args...)
 # - catch something like f(x) for x in domain
-process_generator(gen::Base.Generator{<:Domain}) = (gen.f, gen.iter)
+process_generator(gen::Base.Generator{<:AnyDomain}) = (gen.f, gen.iter)
 # - catch something like f(x,y) for x in domain1, y in domain2
 process_generator(gen::Base.Generator{<:Base.Iterators.ProductIterator}) =
     process_generator(gen, gen.iter.iterators)
-function process_generator(gen, iterators::Tuple{Vararg{Domain}})
-    domain = productdomain(iterators...)
+function process_generator(gen, iterators::Tuple{Vararg{AnyDomain}})
+    domain = productdomain(map(domain, iterators)...)
     dims = map(dimension, iterators)
     (gen.f, domain)
 end
 
 integrate(integrand, args...) =
-    integrate(integrand, process_arguments(args...)...)
+    integrate1(integrand, process_arguments(args...)...)
+
+integrate1(integrand, domain, measure, properties...) =
+    integrate(suggestedstrategy(domain, measure, properties...), integrand, domain, measure, properties...)
 
 integrate(qs::QuadratureStrategy, integrand, args...) =
     integrate_start(qs, integrand, process_arguments(args...)...)
 
-integrate(integrand, domain::Domain, measure::Measure, properties::Property...) =
-    integrate(suggestedstrategy(domain, measure, properties...), integrand, domain, measure, properties...)
-
 
 # The process is as follows:
 # - The argument list is completed (missing domain or measure are added)
@@ -118,7 +118,7 @@ integrate(integrand, domain::Domain, measure::Measure, properties::Property...)
 
 
 integrate_start(qs, integrand, domain, measure, properties...) =
-    integrate_property(qs, FunIntegrand{promote_type(numtype(domain),codomaintype(measure))}(integrand), domain, measure, properties...)
+    integrate_property(qs, FunIntegrand{promote_type(numtype(AsDomain(domain)),codomaintype(measure))}(integrand), domain, measure, properties...)
 
 integrate_start(qs, integrand::Integrand, domain, measure, properties...) =
     integrate_property(qs, integrand, domain, measure, properties...)
@@ -161,7 +161,11 @@ function integrate_domain(qs, integrand, domain, measure, properties...)
         sum_integrals(qs, integrand, domains, measure, properties...)
     else
         transformed_arguments = process_domain(qs, integrand, domain, measure, properties...)
-        integrate_done(qs, transformed_arguments...)
+        if domain_splits(transformed_arguments[2])
+            integrate_domain(qs, transformed_arguments...)
+        else
+            integrate_done(qs, transformed_arguments...)
+        end
     end
 end
 
@@ -191,7 +195,7 @@ fallback_integrate(qs, integrand, domain, measure, properties...) =
 # For quadgk, we only know how to compute intervals
 function apply_quad(qs::Q_quadgk, integrand, domain::AbstractInterval, measure::LebesgueMeasure, properties...)
     if isempty(domain)
-        zero_result(integrand, prectype(domain))
+        zero_result(integrand, prectype(AsDomain(domain)))
     else
         quadgk(integrand, extrema(domain)...; atol = qs.atol, rtol = qs.rtol, maxevals = qs.maxevals)
     end
@@ -261,7 +265,7 @@ function apply_quad(qs::HalfLineRule, integrand, domain::HalfLine, measure::Lebe
 end
 
 # Make sure that a real line rule is applied to a real line integral
-function apply_quad(qs::RealLineRule, integrand, domain::FullSpace, measure::LebesgueMeasure, properties...)
+function apply_quad(qs::RealLineRule, integrand, domain::RealLine, measure::LebesgueMeasure, properties...)
     x = points(qs)
     w = weights(qs)
     z = sum(w[i]*integrand(x[i]) for i in 1:length(x))

src/processing/domain.jl

@@ -1,6 +1,10 @@
 
 ## Process domains
 
+using DomainSets:
+	simplifies,
+	simplify
+
 "A canonical domain for the purposes of evaluating integrals."
 struct CanonicalInt <: DomainSets.CanonicalType
 end
@@ -108,14 +112,20 @@ function process_domain(qs, integrand, domain::DomainSets.TupleProductDomain, me
     error("Don't know how to integrate on a tuple product domain.")
 end
 
-function process_domain(qs, integrand, domain::Domain, measure::Lebesgue, properties...)
-    if has_integration_domain(domain)
-        paramdomain = integration_domain(domain)
-        fmap = mapfrom_integration_domain(domain)
-        process_domain(qs, diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{eltype(paramdomain)}(), properties...)
-    else
-        (integrand, domain, measure, properties...)
-    end
+function process_domain(qs, integrand, domain, measure::Lebesgue, properties...)
+	# first simplify the domain
+	if simplifies(domain)
+		process_domain(qs, integrand, simplify(domain), measure, properties...)
+	else
+		# next find a parametric description of the domain
+	    if has_integration_domain(domain)
+	        paramdomain = integration_domain(domain)
+	        fmap = mapfrom_integration_domain(domain)
+	        process_domain(qs, diffvolume(fmap) * (integrand ∘ fmap), paramdomain, Lebesgue{domaineltype(paramdomain)}(), properties...)
+	    else
+	        (integrand, domain, measure, properties...)
+	    end
+	end
 end
 
 

src/processing/measure.jl

@@ -1,6 +1,9 @@
 
 ## Process measures
 
+integrate_measure(qs, integrand, domain, μ::LebesgueDomain{T}, properties...) where T =
+    integrate_measure(qs, integrand, domain ∩ support(μ), Lebesgue{T}(), properties...)
+
 function integrate_measure(qs, integrand, domain, δ::DiracWeight, properties...)
     x = point(δ)
     if x ∈ domain
@@ -72,7 +75,7 @@ end
 
 
 # Truncate an infinite domain to a finite one for numerical evaluation of Hermite integrals
-function process_measure(qs::AdaptiveStrategy, integrand, domain::FullSpace{T}, measure::HermiteWeight{T}, properties...) where {T}
+function process_measure(qs::AdaptiveStrategy, integrand, domain::Union{RealLine{T},FullSpace{T}}, measure::HermiteWeight{T}, properties...) where {T}
     U = sqrt(-log(eps(T)))
     (hermite_weightfun * integrand, -U..U, Lebesgue{T}(), properties...)
 end
@@ -120,7 +123,7 @@ function process_measure(qs::AdaptiveStrategy, integrand, domain::AbstractInterv
 end
 
 function measure_fetch_transformations(qs, domain, weight, properties...)
-    T = promote_type(numtype(domain),numtype(weight))
+    T = promote_type(numtype(AsDomain(domain)),numtype(weight))
     I = TransformationLogIntegrand{T}()
     I_trans, domain_trans, weight_trans = process_measure(qs, I, domain, weight)
     I_trans.prefactor, I_trans.fun, domain_trans, weight_trans
@@ -161,5 +164,5 @@ integrate_measure(qs::BestRule, integrand, domain::ChebyshevInterval, μ::Jacobi
 integrate_measure(qs::BestRule, integrand, domain::HalfLine, μ::LaguerreWeight{T}, properties...) where {T} =
     integrate_measure(Q_GaussLaguerre(qs.n, μ.α), integrand, domain, Lebesgue{T}(), properties...)
 
-integrate_measure(qs::BestRule, integrand, domain::FullSpace{T}, μ::HermiteWeight{T}, properties...) where {T} =
+integrate_measure(qs::BestRule, integrand, domain::Union{FullSpace{T},RealLine{T}}, μ::HermiteWeight{T}, properties...) where {T} =
     integrate_measure(Q_GaussHermite(qs.n), integrand, domain, Lebesgue{T}(), properties...)

src/processing/properties.jl

@@ -81,7 +81,7 @@ end
 
 "Split the domain according to a singularity at the point `x`."
 function splitdomain_point(x, domain::AbstractInterval)
-    T = promote_type(typeof(x),prectype(domain))
+    T = promote_type(typeof(x),prectype(AsDomain(domain)))
     tol = tolerance(T)
     a, b = extrema(domain)
     if (x > a+tol) && (x < b-tol)
@@ -114,7 +114,7 @@ splitdomain_sing(sing::SingularDiagonal, domain) =
 splitdomain_diagonal(domain::ProductDomain) = splitdomain_diagonal(domain, components(domain)...)
 
 function splitdomain_diagonal(domain::ProductDomain, domain1::AbstractInterval, domain2::AbstractInterval)
-    T = prectype(domain)
+    T = prectype(AsDomain(domain))
     diff1 = domain1 \ domain2
     diff2 = domain2 \ domain1
     overlap = domain1 ∩ domain2

src/weights.jl

@@ -253,6 +253,7 @@ HermiteWeight() = HermiteWeight{Float64}()
 hermite_weightfun(x) = exp(-x^2)
 
 similar(μ::HermiteWeight, ::Type{T}) where {T <: Real} = HermiteWeight{T}()
+support(μ::HermiteWeight{T}) where {T} = RealLine{T}()
 unsafe_weightfun(μ::HermiteWeight, x) = hermite_weightfun(x)
 
 weightfunction(μ::HermiteWeight) = hermite_weightfun
@@ -284,6 +285,7 @@ Display.object_parentheses(μ::GaussianWeight) = true
 lebesguemeasure(domain::UnitInterval{T}) where {T} = LebesgueUnit{T}()
 lebesguemeasure(domain::ChebyshevInterval{T}) where {T} = LegendreWeight{T}()
 lebesguemeasure(domain::FullSpace{T}) where {T} = Lebesgue{T}()
+lebesguemeasure(domain::RealLine{T}) where {T} = Lebesgue{T}()
 lebesguemeasure(domain::Domain{T}) where {T} = LebesgueDomain{T}(domain)
 # avoid propagating Int for intervals with integer eltype
 lebesguemeasure(domain::AbstractInterval{T}) where {T<:Integer} =

test/test_gauss.jl

@@ -30,7 +30,7 @@ function test_gausshermite()
     n = 5
     qs = Q_GaussHermite(n)
     I = integral(x->cos(x)*exp(-x^2), -5..5)
-    @test abs(integral(qs, cos, DomainSets.FullSpace{Float64}()) - I) < 1e-5
+    @test abs(integral(qs, cos, RealLine{Float64}()) - I) < 1e-5
     @test abs(integral(BestRule(n), cos, HermiteWeight()) - I) < 1e-5
 end