Merge 4753d82 into 9b31b53

Project.toml

@@ -17,6 +17,7 @@ GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GridArrays = "356ae075-8bc9-5ef9-a342-46c9ba26438c"
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
 PGFPlotsX = "8314cec4-20b6-5062-9cdb-752b83310925"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

src/BasisFunctions.jl

@@ -23,6 +23,7 @@ using Base: IteratorSize
 using DSP: conv
 using IterativeSolvers: lsqr, lsmr
 using SpecialFunctions: gamma
+using MacroTools: @forward
 import Calculus: derivative
 
 

src/bases/dictionary/dictionary.jl

@@ -35,8 +35,6 @@ of the features for a description of that interface and the syntax.
 abstract type Dictionary{S,T}
 end
 
-FunDict = Dictionary
-
 # Useful abstraction for special cases
 const Dictionary1d{S <: Number,T} = Dictionary{S,T}
 # Warning: not all 2d function sets have NTuple{2,S} type, they could have (S1,S2) type

src/bases/dictionary/expansions.jl

@@ -3,6 +3,7 @@
 compatible_coefficients(dict::Dictionary, coefficients) =
     length(dict) == length(coefficients)
 
+export Expansion
 """
 An `Expansion` describes a function using its expansion coefficients in a certain
 dictionary.
@@ -14,58 +15,56 @@ Parameters:
 - D is the dictionary.
 - C is the type of the expansion coefficients
 """
-struct Expansion{D,C}
+struct Expansion{S,T,D<:Dictionary{S,T},C}
     dictionary      ::  D
     coefficients    ::  C
 
-    function Expansion{D,C}(dict::D, coefficients::C) where {D,C}
+    function Expansion{S,T,D,C}(dict::D, coefficients::C) where {S,T,D<:Dictionary{S,T},C}
         @assert compatible_coefficients(dict, coefficients)
         new(dict, coefficients)
     end
 end
 
 Expansion(dict::Dictionary) = Expansion(dict, zeros(dict))
-Expansion(dict::Dictionary, coef) = Expansion{typeof(dict),typeof(coef)}(dict, coef)
+Expansion(dict::Dictionary, coef) =
+    Expansion{domaintype(dict),codomaintype(dict),typeof(dict),typeof(coef)}(dict, coef)
 
+export Expansion1d, Expansion2d, Expansion3d, Expansion4d
+# Warning: not all 2d function sets have SVector{2,T} type, they could have (S,T) type
+const Expansion1d{S <: Number,T,D,C} = Expansion{S,T,D,C}
+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,
+                                numelements
+
+@forward Expansion.coefficients length, eachindex, getindex, setindex!
+codomaintype(e::Expansion) = promote_type(codomaintype(dictionary(e)), eltype(coefficients(e)))
+isreal(e::Expansion) = isreal(e.dictionary) && isreal(eltype(coefficients(e)))
+
+export expansion
+"""
+    expansion(dict::Dictionary, coefficients)
+
+An expansion of a dictionary with given coefficients.
+"""
 expansion(dict::Dictionary, coefficients) =
     Expansion(dict, native_coefficients(dict, coefficients))
 
 similar(e::Expansion, coefficients) = Expansion(dictionary(e), coefficients)
 
-eltype(::Type{Expansion{D,C}}) where {D,C} = eltype(C)
-
 dictionary(e::Expansion) = e.dictionary
 coefficients(e::Expansion) = e.coefficients
 
-Span(e::Expansion) = Span(dictionary(e))
-
+export random_expansion
 random_expansion(dict::Dictionary) = Expansion(dict, rand(dict))
 
 # For expansions of composite types, return a Expansion of a subdict
 element(e::Expansion, i) = Expansion(element(e.dictionary, i), element(e.coefficients, i))
-
 elements(e::Expansion) = map(Expansion, elements(e.dictionary), elements(e.coefficients))
 
-# Delegation of methods
-for op in (:length, :size, :support, :interpolation_grid)
-    @eval $op(e::Expansion) = $op(dictionary(e))
-end
-
-# Delegation of property methods
-for op in (:numtype, :dimension, :numelements)
-    @eval $op(s::Expansion) = $op(dictionary(s))
-end
-
-hasbasis(e::Expansion) = isbasis(dictionary(e))
-hasframe(e::Expansion) = isframe(dictionary(e))
-
-eachindex(e::Expansion) = eachindex(coefficients(e))
-
-getindex(e::Expansion, i...) = e.coefficients[i...]
-
-setindex!(e::Expansion, v, i...) = setindex!(e.coefficients, v, i...)
-
-
 # This indirect call enables dispatch on the type of the dict of the expansion
 (e::Expansion)(x) = call_expansion(e, dictionary(e), coefficients(e), x)
 (e::Expansion)(x, y) = call_expansion(e, dictionary(e), coefficients(e), SVector(x, y))
@@ -77,12 +76,12 @@ call_expansion(e::Expansion, dict::Dictionary, coefficients, x) =
     eval_expansion(dict, coefficients, x)
 
 function differentiate(e::Expansion, order = difforder(dictionary(e)); options...)
-    op = differentiation(eltype(e), dictionary(e), order; options...)
+    op = differentiation(codomaintype(e), dictionary(e), order; options...)
     Expansion(dest(op), apply(op,e.coefficients))
 end
 
 function antidifferentiate(e::Expansion, order = 1; options...)
-    op = antidifferentiation(eltype(e), dictionary(e); options...)
+    op = antidifferentiation(codomaintype(e), dictionary(e); options...)
     Expansion(dest(op), apply(op,e.coefficients))
 end
 
@@ -114,19 +113,18 @@ roots(f::Expansion) = roots(dictionary(f), coefficients(f))
 
 # Delegate generic operators
 for op in (:extension, :restriction, :evaluation, :approximation, :transform)
-    @eval $op(src::Expansion, dest::Expansion) = $op(promote_type(eltype(src),eltype(dest)), dictionary(src), dictionary(dest))
+    @eval $op(src::Expansion, dest::Expansion) = $op(promote_type(codomaintype(src),codomaintype(dest)), dictionary(src), dictionary(dest))
 end
 
 for op in (:interpolation, :approximation)
     @eval $op(s::Expansion) = $op(dictionary(s))
 end
 
-differentiation(e::Expansion; options...) = differentiation(eltype(e), dictionary(e); options...)
-differentiation(e::Expansion, order; options...) = differentiation(eltype(e), dictionary(e), order; options...)
+differentiation(e::Expansion; options...) = differentiation(codomaintype(e), dictionary(e); options...)
+differentiation(e::Expansion, order; options...) = differentiation(codomaintype(e), dictionary(e), order; options...)
 
 
 show(io::IO, fun::Expansion) = show_setexpansion(io, fun, dictionary(fun))
-
 function show_setexpansion(io::IO, fun::Expansion, fs::Dictionary)
     println(io, "A ", dimension(fun), "-dimensional Expansion with ", length(coefficients(fun)), " degrees of freedom.")
     println(io, "Basis: ", name(fs))
@@ -177,8 +175,10 @@ end
 
 (*)(op::DictionaryOperator, e::Expansion) = apply(op, e)
 
-(*)(a::Number, e::Expansion) = Expansion(dictionary(e), a*coefficients(e))
-(*)(e::Expansion, a::Number) = a*e
+for op in (:+, :-, :*)
+    @eval Base.$op(a::Number, e::Expansion) = Expansion(dictionary(e), $op(a, coefficients(e)))
+    @eval Base.$op(e::Expansion, a::Number) = Expansion(dictionary(e), $op(coefficients(e), a))
+end
 
 apply(op::DictionaryOperator, e::Expansion) = Expansion(dest(op), op * coefficients(e))
 

src/bases/dictionary/gram.jl

@@ -106,7 +106,7 @@ default_gram(Φ::Dictionary, args...; options...) =
 	default_gram(operatoreltype(Φ), Φ, args...; options...)
 
 function default_gram(::Type{T}, Φ::Dictionary, μ::Measure = measure(Φ); options...) where {T}
-    @debug "Slow computation of Gram matrix entrywise."
+    @debug "Slow computation of Gram matrix entrywise of $Φ with measure $μ) ."
     A = grammatrix(Φ, μ, T; options...)
     R = ArrayOperator(A, Φ, Φ)
 end

src/bases/dictionary/gridbasis.jl

@@ -25,9 +25,7 @@ GridBasis(d::Dictionary, g::AbstractGrid = interpolation_grid(d)) = GridBasis{co
 
 grid(b::GridBasis) = b.grid
 
-for op in (:size, :eachindex)
-    @eval $op(b::GridBasis) = $op(grid(b))
-end
+@forward GridBasis.grid size, eachindex
 
 dimension(b::GridBasis) = GridArrays.dimension(grid(b))
 

src/bases/modifiers/tensorproduct_dict.jl

@@ -271,12 +271,9 @@ end
 stencil_parentheses(dict::TensorProductDict) = true
 object_parentheses(dict::TensorProductDict) = true
 
-
 evaluation(::Type{T}, dict::TensorProductDict, gb::GridBasis, grid::ProductGrid; options...) where {T} =
     tensorproduct(map( (d,g) -> evaluation(T, d, g; options...), elements(dict), elements(grid))...)
 
-
-
 dual(dict::TensorProductDict, measure::Union{ProductMeasure,DiscreteProductMeasure}=measure(dict); options...) =
     TensorProductDict([dual(dicti, measurei; options...) for (dicti, measurei) in zip(elements(dict),elements(measure))]...)
 

src/computations/conversion.jl

@@ -28,6 +28,11 @@ conversion
 @deprecate restriction_operator restriction
 @deprecate differentiation_operator differentiation
 @deprecate antidifferentiation_operator antidifferentiation
+@deprecate DictFun Expansion
+@deprecate DictFun1d Expansion1d
+@deprecate DictFun2d Expansion2d
+@deprecate DictFun3d Expansion3d
+@deprecate DictFun4d Expansion4d
 
 operatoreltype(Φ::Dictionary...) = promote_type(map(coefficienttype, Φ)...)
 

src/computations/evaluation.jl

@@ -37,7 +37,7 @@ evaluation(::Type{T}, Φ::Dictionary, grid::AbstractSubGrid; options...) where {
      restriction(T, supergrid(grid), grid; options...) * evaluation(T, Φ, supergrid(grid); options...)
 
 function evaluation(::Type{T}, Φ::Dictionary, gb::GridBasis, grid; options...) where {T}
-    @debug "No fast evaluation available in $grid, using dense evaluation matrix instead."
+    @debug "No fast evaluation available in $grid for dictionary $Φ, using dense evaluation matrix instead."
     dense_evaluation(T, Φ, gb; options...)
 end
 

src/operator/operator.jl

@@ -263,7 +263,12 @@ end#matrix related features
 
 function ≈(op1::DictionaryOperator,op2::DictionaryOperator)
     r = rand(src(op1))
-    op1*r≈op2*r
+    if op1*r≈op2*r
+		return true
+	else
+		@debug "Approx gives difference of $(norm(op1*r-op2*r))"
+		return false
+	end
 end
 
 

src/operator/solvers.jl

@@ -142,30 +142,92 @@ SVD_solver(op::DictionaryOperator; options...) =
 regularized_SVD_solver(op::DictionaryOperator; options...) =
     GenericSolverOperator(op, regularized_svd_factorization(op; options...))
 
+for (damp, LS) in ((:λ, :LSMR), (:damp, :LSQR))
+    OPTIONS = Symbol(LS,"Options")
+    @eval begin
+        export $OPTIONS
+        mutable struct $OPTIONS{T}
+            atol::T
+            btol::T
+            conlim::T
+            $damp::T
+            maxiter::Int
+            verbose::Bool
+
+            function $OPTIONS{T}(;atol=T(1e-6), btol=T(1e-6), conlim=T(Inf), λ=T(0), maxiter=-1, verbose=false) where T
+                new{T}(atol, btol, conlim, λ, maxiter, verbose)
+            end
+            function $OPTIONS(;atol=1e-6, btol=1e-6, conlim=Inf, λ=0, maxiter=-1, verbose=false)
+                T = promote_type(map(typeof, (atol, btol, conlim, λ))...)
+                new{T}(T(atol), T(btol), T(conlim), T(λ), maxiter, verbose)
+            end
+            $OPTIONS{T}(itoptions::$OPTIONS{T}) where T =
+                itoptions
+            $OPTIONS{T}(itoptions::$OPTIONS) where T =
+                $OPTIONS{T}(;convert(NamedTuple, itoptions)...)
+        end
+
+        convert(::Type{<:NamedTuple}, opts::$OPTIONS) =
+            (atol=real(opts.atol), btol=real(opts.btol), conlim=real(opts.conlim), $damp=real(opts.$damp), verbose=real(opts.verbose), maxiter=real(opts.maxiter))
+    end
+end
+
+export ItOptions
+const ItOptions = LSMROptions
+LSQROptions{T}(options::LSMROptions) where T =
+    LSQROptions{T}(;convert(NamedTuple,options)...)
+
 struct LSQR_solver{T} <: VectorizingSolverOperator{T}
     op      ::  DictionaryOperator{T}
-    options ::  NamedTuple
+    options ::  LSQROptions{T}
 
     src_linear  ::  Vector{T}
     dest_linear ::  Vector{T}
-    LSQR_solver(op::DictionaryOperator{T}; atol=1e-6, btol=1e-6, conlim=Inf, damp = 0, maxiter=100*max(size(op)...),verbose=false, options...) where T =
-        new{T}(op, (damp=damp, atol=atol, btol=btol, conlim=conlim, verbose=verbose),
-            Vector{T}(undef, length(dest(op))), Vector{T}(undef, length(src(op))))
+    LSQR_solver(op::DictionaryOperator{T}; itoptions=LSQROptions{T}(), options...) where T =
+        LSQR_solver(op, LSQROptions{T}(itoptions); options...)
+
+
+    function LSQR_solver(op::DictionaryOperator{T}, opts::LSQROptions{T}; verbose=false, options...) where T
+        if opts.maxiter == -1
+             opts.maxiter=100*max(size(op)...)
+        end
+        opts.verbose = verbose
+        new{T}(op, opts,
+            Vector{T}(undef, length(dest(op))), Vector{T}(undef, length(src(op)))
+        )
+    end
+end
+
+function linearized_apply!(op::LSQR_solver, coef_dest::Vector, coef_src::Vector)
+    keys = convert(NamedTuple, op.options)
+    op.options.verbose && @info "LSQR with $keys"
+    copy!(coef_dest, lsqr(op.op, coef_src; keys...))
 end
 
-linearized_apply!(op::LSQR_solver, coef_dest::Vector, coef_src::Vector) =
-    copy!(coef_dest, lsqr(op.op, coef_src; op.options...))
 
 struct LSMR_solver{T} <: VectorizingSolverOperator{T}
     op      ::  DictionaryOperator{T}
-    options ::  NamedTuple
+    options ::  LSMROptions{T}
 
     src_linear  ::  Vector{T}
     dest_linear ::  Vector{T}
-    LSMR_solver(op::DictionaryOperator{T}; atol=1e-6, btol=1e-6, conlim=Inf, damp = 0, maxiter=100*max(size(op)...),verbose=false, options...) where T =
-        new{T}(op, (λ=damp, atol=atol, btol=btol, conlim=conlim, verbose=verbose),
+
+    LSMR_solver(op::DictionaryOperator{T}; itoptions=LSMROptions{T}(), options...) where T =
+        LSMR_solver(op, LSMROptions{T}(itoptions); options...)
+
+    function LSMR_solver(op::DictionaryOperator{T}, opts::LSMROptions{T}; verbose=false, options...) where T
+        if opts.maxiter == -1
+             opts.maxiter=100*max(size(op)...)
+        end
+        opts.verbose = verbose
+        new{T}(op, opts,
             Vector{T}(undef, length(dest(op))), Vector{T}(undef, length(src(op)))
         )
+    end
 end
 
-linearized_apply!(op::LSMR_solver, coef_dest::Vector, coef_src::Vector) = copy!(coef_dest, lsmr(op.op, coef_src; op.options...))
+function linearized_apply!(op::LSMR_solver, coef_dest::Vector, coef_src::Vector)
+    keys = convert(NamedTuple, op.options)
+    op.options.verbose && @info "LSMR with $keys"
+    copy!(coef_dest, lsmr(op.op, coef_src; keys...))
+end

src/util/display/recipes.jl

@@ -4,7 +4,10 @@
 # - evaluate the expansion in the gridpoints (fast if possible)
 # - postprocess the data
 
-@recipe function f(S::Expansion; plot_complex = false, n=200)
+@recipe function f(S::Expansion; plot_complex = false, n=200, plot_ext=false)
+    if plot_ext
+        @warn "Deprecation warning: `plot_ext=true` not longer supported"
+    end
     legend --> false
     # title --> "Expansion"
     grid = plotgrid(dictionary(S), n)
@@ -17,7 +20,7 @@ end
 @recipe function f(S::Expansion, target::Function; n=200)
     grid = plotgrid(dictionary(S), n)
     # Determine the return type so we know where to sample
-    origvals = sample(grid, target, eltype(S))
+    origvals = sample(grid, target, codomaintype(S))
     vals = abs.(origvals - S(grid))
     dictionary(S), grid, postprocess(dictionary(S), grid, vals)
 end
@@ -118,3 +121,20 @@ end
     end
     nothing
 end
+
+
+@recipe function f(dom::Domain2d; n=300, distance=false, xlim=[-1,1], ylim=[-1,1])
+    seriescolor --> :tempo
+     seriestype --> :heatmap
+    aspect_ratio --> 1
+    cbar --> false
+    # xrange = linspace(xlim[1],xlim[2],n)
+    xrange = EquispacedGrid(n,xlim[1],xlim[2])
+    # yrange = linspace(ylim[1],ylim[2],n)
+    yrange = EquispacedGrid(n,ylim[1],ylim[2])
+    # grid = [SVector(i,j) for i in xrange , j in yrange]
+    grid = xrange×yrange
+
+    plotdata = distance ? dist.(grid, Ref(dom)) : in.(grid, Ref(dom))
+    collect(xrange),collect(yrange),plotdata
+end