Initial commit

Project.toml

@@ -1,13 +1,21 @@
 name = "BaryRational"
 uuid = "91aaffc3-5777-4842-85b7-5d3d5d6a3494"
-authors = ["Don MacMillen"]
+authors = ["Don MacMillen <[email protected]> and contributors"]
 version = "0.1.0"
 
+[deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+
 [extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 julia = "1.5"
 
 [targets]
-test = ["Test"]
+test = ["Random", "SpecialFunctions", "Test"]

README.md

@@ -2,3 +2,131 @@
 
 [![Build Status](https://github.com/macd/BaryRational.jl/workflows/CI/badge.svg)](https://github.com/macd/BaryRational.jl/actions)
 [![Coverage](https://codecov.io/gh/macd/BaryRational.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/macd/BaryRational.jl)
+
+"You want poles with that?"
+
+This small package contains both one dimensional barycentric rational
+approximation, using the AAA algorithm [1], and one dimensional
+barycentric rational interpolation with the Floater-Hormann weights
+[2]. The derivatives are calculated using the formulas of [3].
+
+The AAA approximation algorithm can model the poles of a function, if
+present. The FH interpolation is guaranteed to not contain any poles 
+inside of the interpolation interval.
+
+NOTE: Currently, the derivatives are not working.
+
+## Usage
+
+    julia> using BaryRational
+    julia> x = [-3.0:0.1:3.0;];
+    julia> f = x -> sin(x) + 2exp(x)
+    julia> fh = FHInterp(x, f.(x), order=3, grid=true)
+    julia> fh(1.23)
+    7.78493669233287
+
+Note that the default order is 0. The best choice of the order
+parameter appears to be dependent on the number of points (see Table 2
+of [1]) So for smaller data sets, order=3 or order=4 can be good
+choices. This algorithm is not adaptive so you will have to try and see
+what works best for you
+
+If you know that the x points are on an even grid, use grid=true
+
+For approximation using aaa:
+
+    julia> a = aaa(x, f.(x))
+    julia> a(1.23)
+    7.784947874510929
+
+and finally the exact result
+
+    julia> f(1.23)
+    7.784947874511044
+
+The AAA algorithm is adaptive in the subset of support points that it
+chooses to use.
+
+## Examples
+
+Here is an example of fitting f(x) = abs(x) with both FH and AAA. Note
+that because the first derivative is discontinuous at x = 0, we can
+achieve only linear convergence. (Note that systems like Chebfun and
+ApproxFun engineer around this by breaking up the interval at the
+points of discontinuity.)  While the convergence order is the same for
+both algorithms, we see that the AAA has an error that is about a factor
+of 1.6 smaller than the Floater-Hormann scheme.
+
+    using PyPlot
+    using BaryRational
+    function plt_err_abs_x()
+        pts = [40, 80, 160, 320, 640]
+        fh_err  = Float64[]
+        aaa_err = Float64[]
+        order = 3
+        for p in pts
+            xx = collect(range(-5.0, 5.0, length=2p-1))
+            xi = xx[1:2:end]
+            xt = xx[2:2:end]
+            yy = abs.(xi)
+            fa = aaa(xi, yy)
+            fh = FHInterp(xi, yy, order=order, grid=true)
+            push!(aaa_err, maximum(abs.(fa.(xt) .- abs.(xt))))
+            push!(fh_err, maximum(abs.(fh.(xt) .- abs.(xt))))
+        end
+        plot(log.(pts), log.(fh_err), ".-", label="FH Error")
+        plot(log.(pts), log.(aaa_err), ".-", label="AAA Error")
+        xlabel("Log(Number of points)")
+        ylabel("Log(Error)")
+        legend()
+        axis("equal")
+        title("Error in approximating Abs(x)")
+    end
+    plt_err_abs_x()
+
+![image](images/abs_x_error.png)
+
+Since both of these can approximate / interpolate on regular as well as irregular grid
+points they can be used to create ApproxFun Fun's.  ApproxFun needs to be able to evaluate,
+or have evaluated, a function on the Chebyshev points (1st kind here, 2nd kind for Chebfun),
+mostly if you have function values on a regular grid you are out of luck.  Instead, use the
+AAA approximation algorithm to generate an approximation, use that to generate the values on
+the Chebyshev grid, use ApproxFun.transform to transform the function values to coefficients
+and then construct the Fun.  The following shows how.
+
+    using ApproxFun
+    import BaryRational as br
+
+    # our function
+    f(x) = tanh(4x - 1)
+
+    # a regular grid
+    xx = [-1.0:0.01:1.0;];
+
+    # and evaluated on a regular grid
+    yy = f.(xx);
+
+    # and then approximated with AAA
+    faaa = br.aaa(xx, yy);
+
+    # but ApproxFun needs to be evaluated on the Chebyshev points
+    n = 129
+    pts = chebyshevpoints(n);
+    S = Chebyshev();
+
+    # construct the Fun using the aaa approximation on the Chebyshev points
+    pn = Fun(S, ApproxFun.transform(S, faaa.(pts)));
+
+    # now compare it to the "native" fun
+    x = Fun();
+    fapx = tanh(4x - 1);
+    println(norm(fapx - pn))
+
+which yields an error norm of 2.955569189697878e-14. Pretty nice.
+
+
+[1] [The AAA algorithm for rational approximation](http://arxiv.org/abs/1612.00337)
+
+[2] [Barycentric rational interpolation with no poles and high rates of approximation](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.475.3902&rep=rep1&type=pdf)
+
+[3] [Some New Aspects of Rational Interpolation](https://www.ams.org/journals/mcom/1986-47-175/S0025-5718-1986-0842136-8/S0025-5718-1986-0842136-8.pdf)

src/BaryRational.jl

@@ -1,5 +1,32 @@
 module BaryRational
+export FHInterp, bary, chebpts, chebwts, floater_weights, lagrange_weights, aaa, prz, differ, Dcheb
+using LinearAlgebra
+using Statistics
+using SparseArrays
+
+abstract type BRInterp <: Function end
+
+struct FHInterp{T <: AbstractArray} <: BRInterp
+    x::T
+    f::T
+    w::T
+    order::Int
+end
+
+include("weights.jl")
+include("bary.jl")
+include("aaa.jl")
+include("derivatives.jl")
+
+function FHInterp(x::Vector{T}, f::Vector{T}; order::Int=0, grid=false) where {T}
+    if grid
+        FHInterp(collect(x), f, floater_weights(length(x), order), order)
+    else
+        FHInterp(collect(x), f, floater_weights(x, order), order)
+    end
+end
+
+(y::FHInterp)(z) = bary(z, y.f, y.x, y.w)
 
-# Write your package code here.
 
 end

src/aaa.jl

@@ -0,0 +1,170 @@
+# AAA algorithm from the paper "The AAA Alorithm for Rational Approximation"
+# by Y. Nakatsukasa, O. Sete, and L.N. Trefethen, SIAM Journal on Scientific Computing
+# 2018 
+
+struct AAAapprox{T <: AbstractArray} <: BRInterp
+    x::T
+    f::T
+    w::T
+    errvec::T
+end
+
+(a::AAAapprox)(zz) = reval(zz, a.x, a.f, a.w)
+
+# Handle function inputs as well
+function aaa(Z::AbstractArray{T,1}, F::S;  tol=1e-13, mmax=100, verbose=false, clean=true) where {T, S<:Function}
+    aaa(Z, F.(Z), tol=tol, mmax=mmax, verbose=verbose, clean=clean)
+end
+
+"""aaa  rational approximation of data F on set Z
+        r = aaa(Z, F; tol, mmax, verbose, clean)
+
+ Input: Z = vector of sample points
+        F = vector of data values at the points in Z
+        tol = relative tolerance tol, set to 1e-13 if omitted
+        mmax: max type is (mmax-1, mmax-1), set to 100 if omitted
+        verbose: print info while calculating default = false
+        clean: detect and remove Froissart doublets default = true
+
+ Output: r = an AAA approximant as a callable struct with fields
+         z, f, w = vectors of support pts, function values, weights
+         errvec = vector of errors at each step
+
+ NOTE: changes from matlab version:
+         switched order of Z and F in function signature
+         added verbose and clean boolean flags
+         pol, res, zer = vectors of poles, residues, zeros are now only calculated on demand
+         by calling prz(z::AAAapprox)
+"""
+function aaa(Z::AbstractArray{T,1}, F::AbstractArray{S,1}; tol=1e-13, mmax=100,
+             verbose=false, clean=true) where {S, T}
+    # filter out any NaN's or Inf's in the input
+    keep = isfinite.(F)
+    F = F[keep]
+    Z = Z[keep]
+
+    M = length(Z)                    # number of sample points
+    mmax = min(M, mmax)              # max number of support points
+
+    reltol = tol*norm(F, Inf)
+    SF = spdiagm(M, M, 0 => F)       # left scaling matrix
+
+    F, Z = promote(F, Z)
+    P = promote_type(S, T)
+
+    J = [1:M;]
+    z = P[]                          # support points
+    f = P[]                          # function values at support points
+    C = P[]
+    w = P[]
+
+    errvec = P[]
+    R = F .- mean(F)
+    @inbounds for m = 1:mmax
+        j = argmax(abs.(F .- R))               # select next support point
+        push!(z, Z[j])
+        push!(f, F[j])
+        deleteat!(J, findfirst(isequal(j), J))   # update index vector
+
+        # next column of Cauchy matrix
+        C = isempty(C) ? reshape((1 ./ (Z .- Z[j])), (M,1)) : [C (1 ./ (Z .- Z[j]))]
+
+        Sf = diagm(f)                         # right scaling matrix
+        A = SF * C - C * Sf                   # Loewner matrix
+        G = svd(A[J, :])
+
+        # A[J, :] might not have full rank, so svd can return fewer than m columns
+        #w = G.V[:, m]                        
+        w = G.V[:, end]                       # weight vector = min sing vector
+
+        N = C * (w .* f)                      # numerator 
+        D = C * w
+        R .= F
+        R[J] .= N[J] ./ D[J]                  # rational approximation
+
+        err = norm(F - R, Inf)
+        verbose && println("Iteration: ", m, "  err: ", err)
+        errvec = [errvec; err]                # max error at sample points
+        err <= reltol && break                # stop if converged
+    end
+    r = AAAapprox(z, f, w, errvec)
+
+    # remove Frois. doublets if desired.  We do this in place
+    if clean
+        pol, res, zer = prz(r)            # poles, residues, and zeros
+        ii = findall(abs.(res) .< 1e-13)  # find negligible residues
+        length(ii) != 0 && cleanup!(r, pol, res, zer, Z, F)
+    end
+    return r
+end
+
+
+function prz(r::AAAapprox)
+    z, f, w = r.x, r.f, r.w        
+    m = length(w)
+    B = diagm(ones(m+1))
+    B[1, 1] = 0.0
+    E = [0.0  transpose(w); ones(m) diagm(z)]
+    pol, _ = eigen(E, B)
+    pol = pol[isfinite.(pol)] 
+    dz = 1e-5 * exp.(2im*pi*[1:4;]/4)
+
+    # residues
+    res = r(pol .+ transpose(dz)) * dz ./ 4 
+
+    E = [0 transpose(w .* f); ones(m) diagm(z)]
+    zer, _ = eigen(E, B)
+    zer = zer[isfinite.(zer)]
+    pol, res, zer
+end
+
+
+function reval(zz, z, f, w)
+    # evaluate r at zz
+    zv = size(zz) == () ? [zz] : vec(zz)  
+    CC = 1.0 ./ (zv .- transpose(z))           # Cauchy matrix
+    r = (CC * (w .* f)) ./ (CC * w)            # AAA approx as vector
+    r[isinf.(zv)] .= sum(f .* w) ./ sum(w)
+
+    ii = findall(isnan.(r))                    # find values NaN = Inf/Inf if any
+    @inbounds for j in ii
+        if !isnan(zv[j]) && ((v = findfirst(isequal(zv[j]), z)) !== nothing)
+            r[j] = f[v]  # force interpolation there
+        end
+    end
+    r = size(zz) == () ? r[1] : reshape(r, size(zz))             # the AAA approximation
+end
+
+
+# Only calculate the updated z, f, and w
+function cleanup!(r, pol, res, zer, Z, F)
+    z, f, w = r.x, r.f, r.w
+    m = length(z)
+    M = length(Z)
+    ii = findall(abs.(res) .< 1e-13)  # find negligible residues
+    ni = length(ii)
+    ni == 0 && return
+    println("$ni Froissart doublets. Number of residues = ", length(res))
+
+    @inbounds for j = 1:ni
+        azp = abs.(z .- pol[ii[j]] )
+        jj = findall(isequal(minimum(azp)), azp)
+        deleteat!(z, jj)    # remove nearest support points
+        deleteat!(f, jj)
+    end    
+
+    @inbounds for j = 1:length(z)
+        jj = findall(isequal(z[j]), Z)
+        deleteat!(F, jj)
+        deleteat!(Z, jj)
+    end
+    m = m - length(ii)
+    SF = spdiagm(M-m, M-m, 0 => F)
+    Sf = diagm(f)
+    C = 1 ./ (Z .- transpose(z))
+    A = SF*C - C*Sf
+    G = svd(A)
+    w[:] .= G.V[:, m]
+    println("cleanup: ", size(z), "  ", size(f), "  ", size(w))
+    return nothing
+end