Initial commit.

.gitignore

@@ -0,0 +1,4 @@
+*Manifest.toml
+
+*data*
+*build*

Project.toml

@@ -0,0 +1,6 @@
+[deps]
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/Project.toml

@@ -0,0 +1,9 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/make.jl

@@ -0,0 +1,11 @@
+push!(LOAD_PATH,"../src/")
+using Pkg
+Pkg.activate("./")
+using Documenter, Literate, STrAW
+
+# # Generating some documentation files with `Literate.jl`
+# rm("./src/example.md")
+# Literate.markdown("../examples/example.jl", "./src/"; flavor=Literate.DocumenterFlavor())
+
+# Generating documentation with `Documenter.jl`
+makedocs(sitename="STrAW Documentation.")

docs/src/index.md

@@ -0,0 +1,63 @@
+# Features implemented
+
+## Circular waves
+
+```@docs
+b̃p
+βp
+bp
+solution_surrogate
+```
+
+## Herglotz densities
+
+```@docs
+ãp
+wz
+αp
+ap
+```
+
+## Generalized plane waves
+
+```@docs
+gpw
+approximation_set
+```
+
+## Regularized SVD approximation
+
+```@docs
+RegularizedSVDPseudoInverse
+condition_number
+solve_via_regularizedSVD
+```
+
+## Parametric sampling
+```@docs
+TruncKernel
+kernel_diagonal
+probabilitydensityfunction
+cumulativedensityfunction
+Sampler
+inversion
+weight
+random_sampling
+uniform_sampling
+sobol_sampling
+```
+
+## Dirichlet sampling
+```@docs
+samples_from_nodes
+boundary_sampling_nodes
+number_of_boundary_sampling_nodes
+Dirichlet_sampling
+```
+
+## Convenience functions
+
+```@docs
+ϵsupport
+adaptive_G_quad
+```

src/STrAW.jl

@@ -0,0 +1,34 @@
+module STrAW
+
+using LinearAlgebra, SparseArrays, SpecialFunctions
+using FastGaussQuadrature, Roots
+using QuasiMonteCarlo
+
+include("adaptive-quad-rule.jl")
+export ϵsupport, adaptive_G_quad
+
+include("circular-waves.jl")
+export b̃p, βp, bp, solution_surrogate
+
+include("herglotz-densities.jl")
+export ãp, αp, ap, wz, logwz
+
+include("plane-waves.jl")
+export gpw, approximation_set
+
+include("parametric_sampling.jl")
+export TruncKernel, kernel_diagonal
+export probabilitydensityfunction, probabilitydensityfunction_weighted
+export cumulativedensityfunction
+export Sampler, weight, inversion
+export random_sampling, uniform_sampling, sobol_sampling
+
+include("regularizedSVD.jl")
+export RegularizedSVDPseudoInverse, condition_number, solve_via_regularizedSVD
+
+include("dirichlet_sampling.jl")
+export samples_from_nodes, boundary_sampling_nodes
+export number_of_boundary_sampling_nodes
+export Dirichlet_sampling
+
+end

src/adaptive-quad-rule.jl

@@ -0,0 +1,50 @@
+"""
+    function ϵsupport(f, ϵ; δ=1e-3)
+
+Computes the ϵ-support of a function.
+
+This is useful to compute quadratures of compactly ϵ-supported functions on
+unbounded domains.
+"""
+function ϵsupport(f, ϵ; δ=1e-3)
+    a = 0; fa = Inf; while fa > ϵ fa = abs(f(a)); a -= δ end
+    b = 0; fb = Inf; while fb > ϵ fb = abs(f(b)); b += δ end
+    return (a, b)
+end
+
+
+"""
+    adaptive_G_quad(f; quadrule=gausslegendre, a=-1, b=1)
+    
+Computes the integral of ``f`` on ``[a,b]`` using adaptive quadrature rule.
+
+Here adaptivity is only understood with respect to the number of nodes.
+There is no-subdivision of the interval.
+
+!!! note "Note!"
+    The default implementation relies on the
+    [FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)
+    package through the `gausslegendre` quadrature rule function.
+"""
+function adaptive_G_quad(f; quadrule=gausslegendre, a=-1, b=1,
+                         tol=1e-12, Qmin=10, Qstep=1, Qmax=10^3)
+    val = 0
+    err = 1
+    q = Qmin
+    # Estimate
+    x, w = quadrule(q)
+    val = sum((b-a)/2 .* w .* f.((b-a)/2 .* x .+ (a+b)/2))
+    while err > tol && q < Qmax
+        q += Qstep
+        # Reference
+        x, w = quadrule(q)
+        ref = sum((b-a)/2 .* w .* f.((b-a)/2 .* x .+ (a+b)/2))
+        # Error
+        err = abs(val - ref) / ref
+        # Update
+        val = ref
+    end
+    if err > tol @warn "Tolerance not reached: error = $(err)." end
+    if q == Qmax @warn "Maximum number $(Qmax) of nodes reached." end
+    return val, err, q
+end

src/circular-waves.jl

@@ -0,0 +1,89 @@
+"""
+    b̃p(p; k=1)
+    
+Closure for the circular wave function in polar coordinates.
+
+The closure captures the values of the mode number ``p`` as well as the
+wavenumber ``k`` and does not include any normalization.
+The function can then be evaluated at any ``(r,θ)`` point.
+
+The circular waves in polar coordinates are defined by
+```math
+b̃_p := (r,θ) ↦ J_{p}(kr) e^{ı p θ}, \\qquad ∀ p ∈ \\mathbb{Z}.
+```
+"""
+b̃p(p; k=1) = (r,θ) -> besselj(p, k*r) * exp(im * p * θ)
+
+"""
+    βp(p; k=1)
+    
+Computation of the normalization constant of the circular waves.
+
+The circular wave is defined by [`b̃p`](@ref) and the normalization is done
+using the ``k``-weighted ``H^1`` norm.
+
+By definition
+```math
+    \\| b̃_p \\|_{H^{1}}^2 = \\| b̃_p \\|_{L^{2}}^2 + k^{-2} \\| \\nabla b̃_p \\|_{L^{2}}^2.
+```
+Using integration by parts,
+```math
+    \\| \\nabla b̃_p \\|_{L^{2}}^2 = k^2 \\| b̃_p \\|_{L^{2}}^2 + (\\partial_n b̃_p , b̃_p).
+```
+On the one hand
+```math
+    \\| b̃_p \\|_{L^{2}}^2 = \\pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k)).
+```
+On the other hand
+```math
+    (\\partial_n b̃_p , b̃_p) = \\pi k (J_{p-1}(k) - J_{p+1}(k))J_p(k).
+```
+Hence
+```math
+    \\| b̃_p \\|_{H^{1}}^2 = 2\\pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k))
+    + \\pi k^{-1} (J_{p-1}(k) - J_{p+1}(k))J_p(k).
+```
+"""
+function βp(p; k=1)
+    L2nrm2 = π * (besselj(p,k)^2 - besselj(p-1,k) * besselj(p+1,k))
+    H1nrm2 = 2 * L2nrm2 + (π/k) * (besselj(p-1,k) - besselj(p+1,k)) * besselj(p,k)
+    return 1 / sqrt(H1nrm2)
+end
+
+"""
+    bp(p; k=1)
+    
+Closure for the normalized circular waves.
+
+The closure captures the values of the mode number ``p`` and the wavenumber
+``k``.
+The normalization is the default normalization defined in the function
+[`βp`](@ref) and is precomputed once for all.
+The function can then be evaluated at any ``(r,θ)`` point.
+
+The normalized circular waves are defined by
+```math
+b_p := β_p b̃_p, \\qquad ∀ p ∈ \\mathbb{Z}.
+```
+"""
+function bp(p; k=1)
+    β = βp(p; k=k)
+    return (r,θ) -> β * b̃p(p; k=k)(r,θ)
+end
+
+"""
+    solution_surrogate(U; k=1)
+    
+Computes a solution surrogate from a set of coefficients `U`.
+
+The parameter `U` is a vector of coefficients ``u_p`` for ``p`` ranging from
+``-P`` to ``P`` (hence of size ``2P+1``).
+The solution surrogate is then
+```math
+\\mathbf{x} ↦ \\sum_{|p| \\leq P} u_p b_{p}(\\mathbf{x}).
+```
+"""
+function solution_surrogate(U; k=1)
+    P = (length(U) - 1) // 2
+    return (r, θ) -> sum([u * bp(p; k=k)(r, θ) for (p, u) in zip(-P:P, U)])
+end

src/dirichlet_sampling.jl

@@ -0,0 +1,84 @@
+"""
+    samples_from_nodes(f, X)
+    samples_from_nodes(f::Vector, X)
+    
+Evaluate `f` at the sampling nodes `X`.
+
+If `f` is a vector, a matrix is computed.
+"""
+samples_from_nodes(f, X) = [f(x...) for x in X]
+function samples_from_nodes(F::Vector{T}, X) where T
+    return [f(x...) for x in X, f in F]
+end
+
+"""
+    boundary_sampling_nodes(S)
+    
+The sampling nodes on the boundary of the unit disk.
+
+The samplings nodes are then ``(1, θ_s)`` with
+```math
+θ_s := 2π s / S, \\qquad s = 0,…,S-1.
+```
+"""
+boundary_sampling_nodes(S) = [(1, 2π * s / S) for s in 0:S-1]
+
+"""
+    number_of_boundary_sampling_nodes(M; η=2, P=0)
+    
+Number of sampling nodes necessary for an approximation set of dimension ``M``.
+    
+The number of sampling points on the boundary of the unit disk is given by
+``S := \\max(ηM, 2P+1)``.
+The (overdetermined) linear system that will be solved is then of size ``S`` by
+``M``.
+
+The oversampling parameter ``η`` defaults to ``2`` but can be provided by the
+user as an optional parameter `η`.
+The mode number ``P`` defaults to ``0`` but can be provided by the user as an
+optional parameter `P`.
+"""
+number_of_boundary_sampling_nodes(M; η=2, P=0) = max(η * M, 2P+1)
+
+"""
+    Dirichlet_sampling(k, N, U; smpl_type=nothing, η=2, ϵ=1e-8, Q=(length(U)-1)//2)
+
+Example of reconstruction of a solution surrogate via Dirichlet sampling.
+
+The solution surrogate is defined by its vector of coefficients `U`.
+The approximation set of dimension `N` can be composed of propagative plane
+waves (default) or evanescent plane waves (depending on the value of
+`smpl_type` and `Q` the maximum mode number assumed to be in the target).
+The number of sampling points on the boundary can be controlled via the
+oversampling ration `η`.
+The amount of regularization in the SVD can be controlled by the regularization
+parameter `ϵ`.
+"""
+function Dirichlet_sampling(k, U, N; smpl_type=nothing, η=2, ϵ=1e-8,
+                            Q=Integer((length(U)-1)//2))
+    # Maximum mode number in approximation target
+    P = Integer((length(U) - 1) // 2)
+    # Approximation set
+    if smpl_type == nothing # Propagative plane waves
+        Φ = approximation_set(N; k=k)
+    else
+        Φ = approximation_set(N, Q, smpl_type; k=k)
+    end
+    # Sampling nodes
+    S = number_of_boundary_sampling_nodes(N; η=η, P=P)
+    X = boundary_sampling_nodes(S)
+    # Matrix and its factorization
+    A = samples_from_nodes(Φ, X)
+    iA = RegularizedSVDPseudoInverse(A; ϵ=ϵ)
+    # Right-hand-side
+    u = solution_surrogate(U; k=k)
+    b = samples_from_nodes(u, X)
+    # Solution
+    ξ = solve_via_regularizedSVD(iA, b)
+    ũ = (r,θ) -> sum([ξi * ϕi(r,θ) for (ξi, ϕi) in zip(ξ, Φ)])
+    # Errors and stability estimate
+    res = norm(A * ξ - b) / norm(b) # Relative residual
+    nrm = norm(ξ) / norm(U)         # Stability measure
+    e = (r,θ) -> ũ(r,θ) - u(r,θ)    # Absolute error function
+    return u, ξ, ũ, res, nrm, e
+end