From 7d8caf6f4ec9576fa34a49fe334f88832055722c Mon Sep 17 00:00:00 2001 From: Toby Driscoll Date: Tue, 20 Aug 2024 10:57:43 -0400 Subject: [PATCH] update working with different float types --- .gitignore | 1 + Project.toml | 7 +-- src/RationalFunctionApproximation.jl | 2 +- src/aaa.jl | 30 +++++++------ src/approximate.jl | 4 +- src/methods.jl | 15 ++++++- test/interval.jl | 65 ++++++++++++++++------------ test/runtests.jl | 2 +- 8 files changed, 76 insertions(+), 50 deletions(-) diff --git a/.gitignore b/.gitignore index 8c111c9..1fd4145 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ /docs/Manifest.toml /docs/build/ .DS_Store +.vscode diff --git a/Project.toml b/Project.toml index 96c7063..980c935 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.1.4" [deps] ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463" ComplexValues = "41a84b80-6cf2-11e9-379d-9df124847946" +GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PyFormattedStrings = "5f89f4a4-a228-4886-b223-c468a82ed5b9" @@ -28,10 +29,10 @@ julia = "1.9" [extras] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -docs = ["CairoMakie", "ComplexRegions", "DoubleFloats"] -test = ["Test", "LinearAlgebra", "ComplexRegions"] +docs = ["CairoMakie"] +test = ["Test", "LinearAlgebra", "ComplexRegions", "DoubleFloats"] diff --git a/src/RationalFunctionApproximation.jl b/src/RationalFunctionApproximation.jl index d0f119b..336b14f 100644 --- a/src/RationalFunctionApproximation.jl +++ b/src/RationalFunctionApproximation.jl @@ -1,6 +1,6 @@ module RationalFunctionApproximation -using LinearAlgebra, ComplexRegions +using LinearAlgebra, GenericLinearAlgebra, ComplexRegions using PyFormattedStrings export Barycentric, nodes, weights, degree, rewind, diff --git a/src/aaa.jl b/src/aaa.jl index d7ea50c..0e33c5f 100644 --- a/src/aaa.jl +++ b/src/aaa.jl @@ -53,16 +53,17 @@ julia> r(1im * π / 2) See also [`approximate`](@ref) for approximating a function on a curve or region. """ function aaa(z::AbstractVector{<:Number}, y::AbstractVector{<:Number}; - max_degree = 150, float_type = Float64, tol = 1000*eps(float_type), - lookahead = 10, stats = false + max_degree = 150, + float_type = promote_type(typeof(float(1)), eltype(z), eltype(y)), + tol = 1000*eps(float_type), + lookahead = 10, + stats = false ) - @assert float_type <: AbstractFloat - T = float_type fmax = norm(y, Inf) # for scaling m = length(z) iteration = NamedTuple[] - err = T[] + err = float_type[] besterr, bestidx, best = Inf, NaN, nothing # Allocate space for Cauchy matrix, Loewner matrix, and residual @@ -143,17 +144,20 @@ end function aaa( f::Function; - max_degree=150, float_type=Float64, tol=1000*eps(float_type), - refinement=3, lookahead=10, stats=false + max_degree=150, + float_type = promote_type(typeof(float(1)), real_type(f(11//23))), + tol=1000*eps(float_type), + refinement=3, + lookahead=10, + stats=false ) - @assert float_type <: AbstractFloat - T = float_type - CT = Complex{T} + + CT = Complex{float_type} # arrays for tracking convergence progress - err, nbad = T[], Int[] - nodes, vals, pol, weights = Vector{T}[], Vector{CT}[], Vector{CT}[], Vector{CT}[] + err, nbad = float_type[], Int[] + nodes, vals, pol, weights = Vector{float_type}[], Vector{CT}[], Vector{CT}[], Vector{CT}[] - S = [-one(T), one(T)] # initial nodes + S = [-one(float_type), one(float_type)] # initial nodes fS = f.(S) besterr, bestm = Inf, NaN while true # main loop diff --git a/src/approximate.jl b/src/approximate.jl index 812488b..9ea9ce0 100644 --- a/src/approximate.jl +++ b/src/approximate.jl @@ -60,9 +60,10 @@ end approximate(f::Function, d::ComplexCurve; kw...) = approximate(f, Path(d); kw...) approximate(f::Function, d::ComplexClosedCurve; kw...) = approximate(f, ClosedPath(d); kw...) +# all methods end up here function approximate(f::Function, d::ComplexPath; max_degree = 150, - float_type = promote_type(typeof(float(1)), real_type(d)), + float_type = promote_type(real_type(d), typeof(float(1))), tol = 1000*eps(float_type), isbad = z->dist(z, d) < tol, refinement = 3, @@ -70,7 +71,6 @@ function approximate(f::Function, d::ComplexPath; stats = false ) - @assert float_type <: AbstractFloat err, nbad = float_type[], Int[] iteration = NamedTuple[] diff --git a/src/methods.jl b/src/methods.jl index 3211bc6..f719172 100644 --- a/src/methods.jl +++ b/src/methods.jl @@ -70,10 +70,21 @@ roots(f::Approximation) = roots(f.fun) function roots(r::Barycentric) wf = r.w_times_f m = length(wf) - B = diagm( [0; ones(m)] ) + ty = eltype(nodes(r)) + B = diagm( [ty(0); ones(m)] ) # Thanks to Daan Huybrechs: E = [0 transpose(wf); ones(m) diagm(nodes(r))] - return filter(isfinite, eigvals(E, B)) + if ty in (Float64,) + return filter(isfinite, eigvals(E, B)) + else + # super kludgy; thanks to Daan Huybrechs + μ = 11 - 17im # shift since 0 may well be a root + E -= μ*B + EB = inv(E)*B + pp = eigvals(EB) + large_enough = abs.(pp) .> 1e-10 + return μ .+ 1 ./ pp[large_enough] + end end """ diff --git a/test/interval.jl b/test/interval.jl index 6f6d434..5536c69 100644 --- a/test/interval.jl +++ b/test/interval.jl @@ -1,46 +1,58 @@ -pts = 10 .^ range(-15,0,500); pts = [-reverse(pts); 0; pts] -approx(f; kw...) = approximate(f, unit_interval; kw...) -@testset "Basic functions" begin - f = x -> abs(x + 0.5 + 0.01im); @test pass(f, approx(f), pts, atol=5e-13) - f = x -> sin(1/(1.05-x)); @test pass(f, approx(f), pts, atol=2e-13) - f = x -> exp(-1/(x^2)); @test pass(f, approx(f), pts, rtol=4e-13) - f = x -> exp(-100x^2); @test pass(f, approx(f), pts, rtol=2e-13) - f = x -> exp(-10/(1.2-x)); @test pass(f, approx(f), pts, rtol=1e-12) - f = x -> 1/(1+exp(100*(x+.5))); @test pass(f, approx(f), pts, atol=2e-13) - f = x -> sin(100x) * exp(-10x^2); @test pass(f, approx(f), pts, atol=1e-11) - f = x -> abs(x); @test pass(f, approx(f), pts, atol=1e-8) - f = x -> abs(x - 0.95); @test pass(f, approx(f), pts, atol=1e-6) +@testset "Basic functions in $T" for T in (Float64, Double64) + pts = T(10) .^ range(T(-15), T(0), 500); pts = [-reverse(pts); 0; pts] + tol = 2000*eps(T) + approx(f; kw...) = approximate(f, Segment{T}(-1, 1); kw...) + f = x -> abs(x + 1//2 + 1im//100); @test pass(f, approx(f), pts; rtol=tol) + f = x -> sin(1 / (21//20 - x)); @test pass(f, approx(f), pts; rtol=tol) + f = x -> exp(-1 / x^2); @test pass(f, approx(f), pts; rtol=tol) + f = x -> exp(-100x^2); @test pass(f, approx(f), pts; rtol=tol) + f = x -> exp(-10 / (6//5 - x)); @test pass(f, approx(f), pts; rtol=tol) + f = x -> 1/(1 + exp(100*(x + 1//2))); @test pass(f, approx(f), pts; rtol=tol) + f = x -> sin(100x) * exp(-10x^2); @test pass(f, approx(f), pts; rtol=tol) + f = x -> tanh(100*(x - 1//5)); @test pass(f, approx(f), pts, rtol=tol) + f = x -> tanh(100x); @test pass(f, approx(f), pts, rtol=tol) + f = x -> exp(x); @test pass(f, approx(f), pts, rtol=tol) + f = x -> cis(x); @test pass(f, approx(f), pts, rtol=tol) end -@testset "Low-accuracy" begin +@testset "Low accuracy" begin + pts = 10.0 .^ range(-15, 0, 500); pts = [-reverse(pts); 0; pts] + approx(f; kw...) = approximate(f, unit_interval; kw...) f = x -> exp(3x); - @test !pass(f, approximate(f, unit_interval, tol=1e-4), pts, atol=1e-8) + r = approximate(f, unit_interval, tol=1e-5) + @test !pass(f, r, pts, atol=1e-7) + @test pass(f, r, pts, atol=5e-5) + f = x -> abs(x); @test pass(f, approx(f), pts, atol=1e-8) + f = x -> abs(x - 0.95); @test pass(f, approx(f), pts, atol=1e-6) end -@testset "Poles, zeros, residues" begin +@testset "Poles, zeros, residues" for T in (Float64,) + setprecision(80) + approx(f; kw...) = approximate(f, Segment{T}(-1, 1); kw...) f = z -> (z+1) * (z+2) / ((z+3) * (z+4)) r = approx(f) pol = poles(r) zer = roots(r) - @test isapprox(sum(pol+zer), -10, atol=1e-12) + @test isapprox(sum(pol+zer), -10, rtol=5000*eps(T)) f = x -> 2/(3 + x) + 5/(x - 2im); r = approx(f) - @test isapprox(prod(values(residues(r))), 10, atol=1e-8) + @test isapprox(prod(values(residues(r))), 10, rtol=sqrt(eps(T))) f = x -> sinpi(10x); r = approx(f); - @test isapprox(sort(abs.(roots(r)))[19], 0.9, atol=1e-12) + zer = roots(r) + miss = minimum(abs, zer .- 9//10) + @test miss < 1e3*eps(T) f = z -> (z - (3 + 3im))/(z + 2); r = approx(f) pol,zer = poles(r), roots(r) - @test isapprox(pol[1]*zer[1], -6-6im, atol=1e-12) + @test isapprox(pol[1]*zer[1], -6-6im, rtol=5000*eps(T)) end -@testset "Vertical scaling" begin - f = x -> 1e100*sin(x); @test pass(f, approx(f), pts, rtol=2e-13) - f = x -> 1e-100*cos(x); @test pass(f, approx(f), pts, rtol=2e-13) - # f = x -> 1e100*sin(x); @test pass(f, approx(f, =5, lawson=20), pts, rtol=1e-6) - # f = x -> 1e-100*cos(x); @test pass(f, approx(f, =5, lawson=20), pts, rtol=1e-6) +@testset "Vertical scaling in $T" for T in (Float64, Double64) + pts = T(10) .^ range(T(-15), T(0), 500); pts = [-reverse(pts); 0; pts] + f = x -> T(10)^50*sin(x); @test pass(f, approx(f), pts, rtol=2000*eps(T)) + f = x -> T(10)^(-50)*cos(x); @test pass(f, approx(f), pts, rtol=2000*eps(T)) end # @testset "Lawson" begin @@ -72,10 +84,6 @@ end f = x -> 1/(1+1im*x); @test pass(f, approx(f, max_degree=3), pts, atol=2e-13) f = x -> 1/(3+x+x^2); @test pass(f, approx(f, max_degree=2), pts, atol=2e-13) f = x -> 1/(1.01+x^3); @test pass(f, approx(f, max_degree=3), pts, atol=2e-13) - f = x -> tanh(100x); @test pass(f, approx(f), pts, atol=2e-13) - f = x -> tanh(100*(x-.2)); @test pass(f, approx(f), pts, atol=2e-13) - f = x -> exp(x); @test pass(f, approx(f, tol=1e-13), pts, atol=2e-13) - f = x -> cis(x); @test pass(f, approx(f), pts, atol=2e-13) f = x -> cis(x); r = approx(f); @test sum(@. abs(r(pts))-1)/length(pts) < 2e-13 f = x -> exp(exp(x))/(x - 0.2im); r = approx(f); @@ -83,6 +91,7 @@ end end @testset "Other intervals" begin + pts = 10 .^ range(-15, 0, 500); pts = [-reverse(pts); 0; pts] zz(a,b) = (pts .+ 1)*(b-a)/2 .+ a for (a,b) in ((-2,3), (0,1), (-0.01,0) , (-1e4, 2e6)) for f in ( x -> abs(x + 0.5 + 0.01im), x -> sin(1/(1.05im-x)), x -> exp(-10/(1.2*(b+1)-x)) ) diff --git a/test/runtests.jl b/test/runtests.jl index 7209ee6..060cb3a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -using RationalFunctionApproximation, Test, LinearAlgebra, ComplexRegions +using RationalFunctionApproximation, Test, LinearAlgebra, ComplexRegions, DoubleFloats # ii = 1im * range(-300,400,500) # zz = cispi.(2 * (0:500) / 500)