Initial implementation for the inverse of a function

Author: daanhb

Project.toml

@@ -1,7 +1,7 @@
 name = "BasisFunctions"
 uuid = "4343a256-5453-507d-8aad-01a9d7189916"
 authors = ["Daan Huybrechs <daan.huybrechs@kuleuven.be>"]
-version = "0.7"
+version = "0.7.1"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"

src/bases/dictionary/basisfunction.jl

@@ -89,6 +89,8 @@ innerproduct(φ::AbstractBasisFunction, ψ::AbstractBasisFunction, measure; opti
 norm(φ::AbstractBasisFunction; options...) = norm(φ, measure(φ); options...)
 norm(φ::AbstractBasisFunction, p; options...) = dict_norm(dictionary(φ), index(φ), p; options...)
 
+moment(φ::AbstractBasisFunction; options...) = dict_moment(dictionary(φ), index(φ); options...)
+
 # The inner product of a basis function with another function: this is an analysis integral
 # We introduce a separate function name for this for easier dispatch.
 innerproduct(φ::AbstractBasisFunction, g, measure; options...) =

src/bases/dictionary/expansions.jl

@@ -35,8 +35,11 @@ const Expansion2d{S <: Number,T,D,C} = Expansion{SVector{2,S},T,D,C}
 const Expansion3d{S <: Number,T,D,C} = Expansion{SVector{3,S},T,D,C}
 const Expansion4d{S <: Number,T,D,C} = Expansion{SVector{4,S},T,D,C}
 
-@forward Expansion.dictionary domaintype, prectype, dimension,
-                                size, Span, support, interpolation_grid, numtype,
+prectype(F::Type{<:Expansion{S,T}}) where {S,T} = prectype(S,T)
+numtype(F::Type{<:Expansion{S,T}}) where {S,T} = numtype(S,T)
+
+@forward Expansion.dictionary domaintype, dimension,
+                                size, Span, support, interpolation_grid,
                                 ncomponents,
                                 measure
 
@@ -127,7 +130,14 @@ adjoint(f::Expansion) = differentiate(f)
 
 ∫(f::Expansion) = antidifferentiate(f)
 
-roots(f::Expansion) = expansion_roots(dictionary(f), coefficients(f))
+function roots(f::Expansion; all = true)
+    r = expansion_roots(dictionary(f), coefficients(f))
+    if all
+        r
+    else
+        restrict_to_domain(r, support(f))
+    end
+end
 
 show(io::IO, mime::MIME"text/plain", fun::Expansion) = composite_show(io, mime, fun)
 Display.displaystencil(fun::Expansion) = ["Expansion(", dictionary(fun), ", ", coefficients(fun), ")"]

src/computations/approximation.jl

@@ -132,3 +132,38 @@ end
 
 default_leastsquares(Φ::Dictionary, gb::GridBasis, T=operatoreltype(Φ,gb); options...) =
     QR_solver(ArrayOperator(leastsquares_matrix(Φ, grid(gb), T), Φ, gb))
+
+
+########################
+# Inverse of a function
+########################
+
+"""
+    inversefunction(F[, domain; n])
+
+Compute an approximation to the inverse of `F`, using an appropriate basis with `n` degrees of freedom.
+"""
+inversefunction(F::Expansion; args...) = inversefunction(F, support(F); args...)
+
+function inversefunction(F, domain::AbstractInterval; n = 2*length(F))
+    Fd = differentiate(F)
+    if length(roots(Fd; all=false)) > 0
+        throw(ArgumentError("Function is not invertible on its domain."))
+    end
+    
+    # compute the value t such that F(t)=z.
+    function inversevalue(z)
+        r = roots(F - z; all=false)
+        if length(r) != 1
+            throw(ArgumentError("No unique inverse value."))
+        end
+        r[1]
+    end
+    
+    a = leftendpoint(domain)
+    b = rightendpoint(domain)
+    A = F(a)
+    B = F(b)
+    G = ChebyshevT(n) → (A..B)
+    Expansion(G, approximation(G) * inversevalue)
+end

src/util/common.jl

@@ -132,3 +132,14 @@ export Domain1d, Domain2d, Domain3d, Domain4d
 iscompatible(map1::AbstractMap, map2::AbstractMap) = map1==map2
 iscompatible(map1::AffineMap, map2::AffineMap) = (map1.A ≈ map2.A) && (map1.b ≈ map2.b)
 iscompatible(domain1::Domain, domain2::Domain) = domain1 == domain2
+
+"Filter the given points by leaving out all points that are not an element of the given domain."
+function restrict_to_domain(x::AbstractVector, d; tol = DomainSets.domain_tolerance(d))
+    checkdomain(d)
+    filter(t->approx_in(t, d, tol), x)
+end
+
+function restrict_to_domain(x::AbstractVector{Complex{T}}, d::AbstractInterval{T}; tol = DomainSets.domain_tolerance(d)) where {T}
+    x1 = real(filter(t->abs(imag(t))<tol, x))
+    restrict_to_domain(x1, d; tol)
+end