An axiomatically-designed constructor

myFun.jl

using ApproxFun
import ApproxFun: Fun, domain, Space, checkpoints, ArraySpace, zerocfsFun, defaultFun, real

function myFun(f, d::Space)
    #TODO: reuse function values?
    T = eltype(domain(d))
    if T <: Complex
        T = T.parameters[1] #get underlying real representation
    end
    r=checkpoints(d)
    f0=f(first(r))

    if !isa(d,ArraySpace) && isa(f0,Array)
        return zerocfsFun(f,ArraySpace(d,size(f0)...))
    end

    tol =T==Any?eps():eps(T)

    fr=map(f,r)
    maxabsfr=norm(fr,Inf)

    n = 16
    while n ≤ 2^20
        cf = defaultFun(f, d, n)
        maxabsc = norm(cf.coefficients,Inf)
        if maxabsc==0 && maxabsfr==0
            return(zeros(d))
        end

        # we allow for transformed coefficients being a different size
        ##TODO: how to do scaling for unnormalized bases like Jacobi?
        if length(cf) > 8 && norm(slice(cf.coefficients,length(cf)>>2:length(cf)),Inf) < length(cf)*maxabsc*tol && sampletest(cf,r,fr)
            return pad(cf,findend(cf.coefficients,maxabsc*tol/10))
        end
        n = n << 1
    end

    warn("Maximum length "*string(n)*" reached")

    Fun(f,d,n)
end

function findend(cfs::Vector,tol::Real)
    N = length(cfs)
    k = N
    while k > 1
        if abs(cfs[k]) > k*log2(k)*tol break end
        k-=1
    end
    k
end

function sampletest(cf::Fun,r::Vector,fr::Vector)
    ret = true
    for k=1:length(r) ret *= norm(cf(r[k])-fr[k],1)<1E-4norm(fr[k],1) end
    ret
end


## My Tests

ω = logspace(0,4,30);
l1 = zeros(30);
l2 = zeros(30);
l3 = zeros(30);
x = Fun()
for j=1:30
    f1 = Fun(x->sinpi(ω[j]*x),Chebyshev());
    f2 = myFun(x->sinpi(ω[j]*x),Chebyshev());
    f3 = sinpi(ω[j]*x);
    l1[j] = length(f1);
    l2[j] = length(f2);
    l3[j] = length(f3);
end
Main.Plots.plot(ω,l1;xscale=:log10,yscale=:log10,label="Current")
Main.Plots.plot!(ω,l2;xscale=:log10,yscale=:log10,label="Proposed")
Main.Plots.plot!(ω,l3;xscale=:log10,yscale=:log10,label="Override")
Main.Plots.xlabel!("\$\\omega\$")
Main.Plots.ylabel!("\$N\$")
Main.Plots.title!("\$f(x) = \\cos(\\omega x)\$")
Main.Plots.savefig("Oscillatory test.pdf")

ɛ = logspace(0,-4,30);
l1 = zeros(30);
l2 = zeros(30);
l3 = zeros(30);
x = Fun()
for j=1:30
    f1 = Fun(x->ɛ[j]/(x^2+ɛ[j]),Chebyshev());
    f2 = myFun(x->ɛ[j]/(x^2+ɛ[j]),Chebyshev());
    f3 = ɛ[j]/(x^2+ɛ[j]);
    l1[j] = length(f1);
    l2[j] = length(f2);
    l3[j] = length(f3);
end
Main.Plots.plot(ɛ,l1;xscale=:log10,yscale=:log10,label="Current")
Main.Plots.plot!(ɛ,l2;xscale=:log10,yscale=:log10,label="Proposed")
Main.Plots.plot!(ɛ,l3;xscale=:log10,yscale=:log10,label="Override")
Main.Plots.xlabel!("\$\\varepsilon\$")
Main.Plots.ylabel!("\$N\$")
Main.Plots.title!("\$f(x) = \\frac{\\varepsilon}{x^2+\\varepsilon}\$")
Main.Plots.savefig("Singularity test.pdf")


n = map(x->round(Int,x),logspace(0,5,30));
l1 = zeros(30);
l2 = zeros(30);
x = Fun()
for j=1:30
    f1 = Fun(x->cos(n[j]*acos(x)),Chebyshev());
    f2 = myFun(x->cos(n[j]*acos(x)),Chebyshev());
    l1[j] = length(f1);
    l2[j] = length(f2);
end
Main.Plots.plot(n,(l1-1)./n;xscale=:log10,yscale=:log10,label="Current")
Main.Plots.plot!(n,(l2-1)./n;xscale=:log10,yscale=:log10,label="Proposed")
Main.Plots.ylims!(1e-1,1e1)
Main.Plots.xlabel!("\$n\$")
Main.Plots.ylabel!("\$(N-1)/n\$")
Main.Plots.title!("\$f(x) = T_n(x)\$")
Main.Plots.savefig("Basis test.pdf")