Changing defaults and calling sequences

.gitignore

@@ -7,3 +7,4 @@ test/Manifest.toml
 docs/node_modules
 *.cov
 benchmark/Manifest.toml
+.idea

CLAUDE.md

@@ -55,10 +55,10 @@ ext/
 
 ### Concrete rational function types
 
-**`Barycentric{T,S}`** (barycentric.jl:13) — default method; alias `AAA = Barycentric`  
+**`Barycentric{T,S}`** (barycentric.jl:13) — alias `AAA = Barycentric`  
 Fields: `nodes`, `values`, `weights`, `w_times_f`
 
-**`Thiele{T,S}`** (thiele.jl) — continued fraction representation; alias `TCF = Thiele`  
+**`Thiele{T,S}`** (thiele.jl) — default method; continued fraction representation; alias `TCF = Thiele`  
 Fields: `nodes`, `values`, `weights`  
 Two evaluation strategies: classic and numerically stable `onediv`.
 
@@ -73,14 +73,15 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 - **`ContinuumApproximation{T,S,R}`** (approximation.jl:37) — continuous domain wrapper; fields: `original`, `domain`, `fun`, `allowed`, `path`, `history`
 - **`DiscreteApproximation{T,S,R}`** (approximation.jl:75) — discrete point set wrapper; fields: `data`, `domain`, `fun`, `test_index`, `allowed`, `history`
 - **`IterationRecord{R,S,T}`** (approximation.jl:5) — convergence history entry; fields: `interpolant`, `error`, `poles`
+- **`ConvergenceStatus`** (approximation.jl:45) — why an iteration stopped; fields: `reason`, `best`, `iterations`, `error`. `reason` is one of `:converged`, `:stagnated`, `:max_degree`, `:node_failure`, `:nan_weight`, `:refinement`, `:exhausted`, `:rewound`
 
 ---
 
 ## Public API
 
 ### Approximation construction
-- `approximate(f, domain; method, max_iter, tol, allowed, refinement, stagnation)` — main entry point
-- `approximate(f, domain, poles)` — least-squares with prescribed poles
+- `approximate(f, domain, method=Thiele(); max_iter, tol, allowed, refinement, stagnation)` — main entry point; `method` is an instance (`Thiele()` default, `Barycentric()`) passed as the last positional argument
+- `approximate(f, domain, poles)` — least-squares with prescribed poles (selectable via `PartialFractions()` as the last positional argument)
 - `aaa(y, z; kwargs...)` — legacy discrete AAA (deprecated)
 
 ### Rational function queries
@@ -98,6 +99,8 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 - `get_function(r)`, `domain(r)` — extract components
 - `rewind(r, index)` — revert to earlier iteration
 - `get_history(r)` — convergence history
+- `status(r)` — `ConvergenceStatus` for the run, or `nothing` if none was recorded
+- `isconverged(r)` — whether the iteration reached `tol`, as opposed to stagnating or exhausting `max_degree`
 - `test_points(r)` — test point locations
 
 ### Optimization
@@ -157,10 +160,12 @@ Fields: `polynomial::ArnoldiPolynomial`, `poles`, `residues`
 ## Design Decisions
 
 1. **Two-parameter type system**: `T` for float precision, `S` for value type — supports generic arithmetic.
-2. **Barycentric as default**: most efficient/stable; aliased `AAA` for historical compatibility.
+2. **Thiele as default**: continued-fraction method used when no selector is given; `Barycentric` (aliased `AAA`) remains available and selectable.
 3. **Continuum vs. Discrete split**: `ContinuumApproximation` and `DiscreteApproximation` reflect fundamentally different strategies.
 4. **Adaptive path discretization**: `DiscretizedPath` stores multiple refinement levels in matrix form.
 5. **`allowed` parameter**: generic function to filter pole locations; enables multiply-connected domains.
 6. **Convergence history**: optional recording enables `rewind()` and convergence plots.
+   `quitting_check` returns a `(reason, best)` tuple rather than an overloaded integer, and
+   `best_acceptable` is callable on its own so failure paths need not fake a `max_iter`.
 7. **Extension architecture**: plotting and autodiff are optional — no hard dependencies.
 8. **Precompilation workload**: uses `@compile_workload` for fast time-to-first-approximation.

Project.toml

@@ -9,6 +9,7 @@ ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
 ComplexValues = "41a84b80-6cf2-11e9-379d-9df124847946"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 GenericSchur = "c145ed77-6b09-5dd9-b285-bf645a82121e"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -37,6 +38,7 @@ ComplexValues = "0.3"
 ForwardDiff = "1"
 GenericLinearAlgebra = "0.3, 0.4"
 GenericSchur = "0.5"
+IntervalSets = "0.7"
 Logging = "1"
 Makie = "0.24"
 Plots = "1"
@@ -52,10 +54,11 @@ julia = "1"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ComplexRegions = "c64915e2-6c82-11e9-38e9-1f159a780463"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ReTest = "e0db7c4e-2690-44b9-bad6-7687da720f89"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 docs = ["CairoMakie"]
-test = ["Test", "ReTest", "LinearAlgebra", "ComplexRegions", "DoubleFloats"]
+test = ["Test", "ReTest", "LinearAlgebra", "ComplexRegions", "DoubleFloats", "IntervalSets"]

benchmark/benchmarks.jl

@@ -4,15 +4,27 @@ using Logging
 # silence convergence warnings so benchmark output stays clean
 global_logger(SimpleLogger(stderr, Logging.Error))
 
+# AirspeedVelocity runs this one script against both the base revision and the PR head,
+# so it has to work with either calling convention: the old API selected the interpolant
+# type with a `method = Barycentric` keyword, the current one takes an instance
+# (`Barycentric()`) as the last positional argument. Pick the form once, at load time, so
+# no branch survives into the benchmarked call. Note the deprecation shim is deliberately
+# *not* used here: it warns on every call, which would skew the timings it appears in.
+if hasmethod(Barycentric, Tuple{})
+    approx(f, domain, M; kw...) = approximate(f, domain, M(); kw...)
+else
+    approx(f, domain, M; kw...) = approximate(f, domain; method = M, kw...)
+end
+
 const SUITE = BenchmarkGroup()
 
 # --- construction cost, by method and domain ---
 SUITE["approximate"] = BenchmarkGroup()
 for (name, method) in (("aaa", Barycentric), ("thiele", Thiele))
     g = SUITE["approximate"][name] = BenchmarkGroup()
-    g["exp_interval"] = @benchmarkable approximate(exp, $unit_interval; method = $method, allowed = true)
-    g["tanh_steep"]   = @benchmarkable approximate(x -> tanh(50x), $unit_interval; method = $method, allowed = true)
-    g["abs_circle"]   = @benchmarkable approximate(z -> abs(z - 1.0001im), $unit_circle; method = $method, allowed = true)
+    g["exp_interval"] = @benchmarkable approx(exp, $unit_interval, $method; allowed = true)
+    g["tanh_steep"]   = @benchmarkable approx(x -> tanh(50x), $unit_interval, $method; allowed = true)
+    g["abs_circle"]   = @benchmarkable approx(z -> abs(z - 1.0001im), $unit_circle, $method; allowed = true)
 end
 
 # --- construction cost on a discrete point set ---
@@ -23,26 +35,26 @@ end
 SUITE["approximate_discrete"] = BenchmarkGroup()
 for (name, method) in (("aaa", Barycentric), ("thiele", Thiele))
     g = SUITE["approximate_discrete"][name] = BenchmarkGroup()
-    g["tanh_steep"]  = @benchmarkable approximate(x -> tanh(100x), $DISCRETE_PTS; method = $method, allowed = true)
-    g["abs_shift"]   = @benchmarkable approximate(x -> abs(x + 0.5 + 0.01im), $DISCRETE_PTS; method = $method, allowed = true)
-    g["sin_recip"]   = @benchmarkable approximate(x -> sin(1 / (1.05 - x)), $DISCRETE_PTS; method = $method, allowed = true)
+    g["tanh_steep"]  = @benchmarkable approx(x -> tanh(100x), $DISCRETE_PTS, $method; allowed = true)
+    g["abs_shift"]   = @benchmarkable approx(x -> abs(x + 0.5 + 0.01im), $DISCRETE_PTS, $method; allowed = true)
+    g["sin_recip"]   = @benchmarkable approx(x -> sin(1 / (1.05 - x)), $DISCRETE_PTS, $method; allowed = true)
 end
 
 # --- evaluation cost on a fixed approximant ---
 SUITE["evaluate"] = BenchmarkGroup()
-let r = approximate(x -> tanh(50x), unit_interval, method = Barycentric, allowed = true),
+let r = approx(x -> tanh(50x), unit_interval, Barycentric; allowed = true),
     z = collect(range(-1, 1, 1000))
 
     SUITE["evaluate"]["bary_vector"] = @benchmarkable $r.($z)
 end
-let r = approximate(x -> tanh(50x), unit_interval; method = Thiele, allowed = true),
+let r = approx(x -> tanh(50x), unit_interval, Thiele; allowed = true),
     z = collect(range(-1, 1, 1000))
 
     SUITE["evaluate"]["thiele_vector"] = @benchmarkable $r.($z)
 end
 
 # --- pole solve ---
 SUITE["poles"] = BenchmarkGroup()
-let r = approximate(x -> 1 / sqrt(x^2 + 0.01), unit_interval; method = Thiele)
+let r = approx(x -> 1 / sqrt(x^2 + 0.01), unit_interval, Thiele)
     SUITE["poles"]["thiele"] = @benchmarkable poles($r)
 end

docs/Project.toml

@@ -15,3 +15,4 @@ RationalFunctionApproximation = {path = ".."}
 [compat]
 Documenter = "1"
 DomainColoring = "2"
+DocumenterCitations = "~1.4"

docs/src/algorithms.md

@@ -26,8 +26,8 @@ The AAA algorithm is the best-known and most widely used method for rational app
 The `convergenceplot` function shows the errors of the approximants found during the AAA iteration.
 
 ```@example convergence
-f = x -> cos(exp(3x))
-r = approximate(f, unit_interval)
+f(x) = cos(exp(3x))
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
@@ -60,56 +60,56 @@ That's why the reported error is not very large. It's worth keeping in mind that
 It is possible for the iteration to stagnate with bad poles if the original function has a singularity very close to the domain.
 
 ```@example convergence
-f = x -> tanh(500*(x - 1//4))
-r = approximate(f, unit_interval)
+f(x) = tanh(500*(x - 1//4))
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
 This effect is thought to be mainly due to roundoff and conditioning of the problem. In this case, if we use more accurate floating-point arithmetic, we can see that the AAA convergence continues steadily past the previous plateau. In the following, we apply `Double64` arithmetic, having used exact rational numbers already in the definition of `f`:
 
 ```@example convergence
 using DoubleFloats, ComplexRegions
-r = approximate(f, Segment{Double64}(-1, 1))
+r = approximate(f, Double64(-1)..1, AAA())
 convergenceplot(r)
 ```
 
 In the extreme case of a function with a singularity on the domain, the convergence can be substantially affected:
 
 ```@example convergence
-f = x -> abs(x - 1/8)
-r = approximate(f, unit_interval)
+f(x) = abs(x - 1/8)
+r = approximate(f, -1..1, AAA())
 convergenceplot(r)
 ```
 
 In such a case, we might get improvement by increasing the number of allowed consecutive failures via the `stagnation` keyword argument:
 
 ```@example convergence
-r = approximate(f, unit_interval, stagnation=50)
+r = approximate(f, -1..1, AAA(); stagnation=50)
 convergenceplot(r)
 ```
 
 However, AAA is an $O(n^4)$ algorithm, so venturing into higher degrees can become costly.
 
 ## Thiele continued fractions (TCF)
 
-The TCF algorithm [SalazarCelisNumericalContinued2024](@cite) is much newer than AAA and less thoroughly battle-tested, even though it's based on a continued fraction representation of rational functions that is over a century old. Like AAA, it uses iterative greedy node selection, and the effects of that ordering look good in experiments so far but are poorly understood theoretically. In TCF's favor are its $O(n^3)$ complexity requirement and an algorithmic simplicity that requires nothing more than basic arithmetic.
+The greedy TCF algorithm [SalazarCelisNumericalContinued2024](@cite) is much newer than AAA despite being based on a continued fraction representation of rational functions that is over a century old. Like AAA, it uses iterative greedy node selection, and the effects of that ordering look very good in experiments. In TCF's favor are its $O(n^3)$ complexity requirement and an algorithmic simplicity that requires nothing more than basic arithmetic.
 
-To try greedy TCF, use `method=Thiele` or `method=TCF` as an argument to `approximate`.
+As of version 0.4 of the package, TCF is the default method used by `approximation`. You can also select it manually by passing `Thiele()` (or its alias `TCF()`) as the third positional argument.
 
 ```@example convergence
-f = x -> cos(41x - 5) * exp(-10x^2)
-r = approximate(f, unit_interval; method=TCF)
+f(x) = cos(41x - 5) * exp(-10x^2)
+r = approximate(f)    # -1..1 and TCF() by default
 convergenceplot(r)
 ```
 
-The $x$-axis of the convergence plot shows the degree of the denominator polynomial. Because the Thiele method alternates between interpolants of type $(n, n)$ and $(n+1, n)$, there are two dots in the plot for each degree. The dots corresponding to approximations of the diagonal and superdiagonal rational type are connected by lines; sometimes, they could be viewed as separate convergence curves.
+The $x$-axis of the convergence plot shows the degree of the denominator polynomial. Because the Thiele method alternates between interpolants of type $(n, n)$ and $(n+1, n)$, there are two dots in the plot for each degree, and line segments connect the dots sharing the same rational type, giving two convergence curves.
 
-Because TCF uses only addition, multiplication, and division, it is easy to use in extended precision arithmetic. Here, we use `allowed=true` to disable checking for poles, because doing so requires solving an eigenvalue problem that is far more expensive than the iteration itself.
+Because TCF uses only addition, multiplication, and division, it is easy to use in extended precision arithmetic. 
 
 ```@example convergence
-f = x -> atan(1e5*(x - 1//2))
-domain = Segment{BigFloat}(-1, 1)
-@elapsed r = approximate(f, domain; method=TCF, max_iter=400, allowed=true, stagnation=40)
+f(x) = atan(1e5*(x - 1//2))
+domain = big(-1)..1    # use BigFloats
+@elapsed r = approximate(f, domain; max_degree=150, stagnation=40)
 ```
 
 ```@example convergence
@@ -129,9 +129,9 @@ When posed on a discrete set of test points, this is a linear least-squares prob
 There is no iteration on the degree of the polynomial or rational parts of the approximant. In the continuum variant, though, the discretization of the boundary of the domain is refined iteratively until either the max-norm error is below a specified threshold or has stopped improving.
 
 ```@example convergence
-f = x -> tanh(x)
+f(x) = tanh(x)
 ζ = 1im * π * [-1/2, 1/2, -3/2, 3/2]
-r = approximate(f, Segment(-2, 2), ζ)
+r = approximate(f, -2..2, ζ)
 ```
 
 ```@example convergence
@@ -142,7 +142,7 @@ println("Max error: $(max_err(r))")
 To get greater accuracy, we can increase the degree of the polynomial part.
 
 ```@example convergence
-r = approximate(f, Segment(-2, 2), ζ; degree=20)
+r = approximate(f, -2..2, ζ; degree=20)
 max_err(r)
 ```
 
@@ -152,16 +152,16 @@ Note that the residues, which are all equal to 1 for the exact function, may not
 Pair.(residues(r)...)
 ```
 
-Suppose now we approximate $|x|$ using AAA. We can extract the poles of the result.
+Suppose now we approximate $|x|$ using TCF. We can extract the poles of the result, filtering out those that lie on the real axis.
 
 ```@example convergence
-r = approximate(abs, unit_interval, tol=1e-9)
-ζ = poles(r)
+r = approximate(abs; tol=1e-9, stagnation=30)
+ζ = filter(z -> abs(imag(z)) > 1e-8, poles(r))
 ```
 
-To what extent might these poles be suitable for a different function that has the same singularity?
+These poles might be suitable for a different function that has the same singularity:
 
 ```@example convergence
-s = approximate(x -> exp(abs(x)), unit_interval, ζ; degree=20)
+s = approximate(x -> exp(abs(x)), -1..1, ζ; degree=20)
 max_err(r)
 ```

docs/src/discrete.md

@@ -5,9 +5,9 @@ For many functions, discretization of the domain is straightforward. But if the
 The `approximate` function can take a vector of sample points as a domain. The given function is then evaluated only at those points, and the rational approximation is a fully discrete process that uses only the given data.
 
 ```@example mode
-using RationalFunctionApproximation, ComplexRegions
+using RationalFunctionApproximation
 x = -1:0.01:1
-f = x -> tanh(5 * (x - 0.2))
+f(x) = tanh(5 * (x - 0.2))
 r = approximate(f, x)
 ```
 You can alternatively provide just the discrete function values yourself. The domain is always given second:
@@ -20,6 +20,7 @@ r = approximate(y, x)
 As long as there are no singularities as close to the domain as the sample points are to one another, a basic discretization works well.
 
 ```@example mode
+using ComplexRegions    # to get dist() function
 I = unit_interval
 println("nearest pole is $(minimum(dist(z, I) for z in poles(r))) away")
 _, err = check(r);
@@ -29,7 +30,7 @@ println("max error on the given domain: ", maximum(abs, err))
 But if the distance to a singularity is comparable to the sample spacing, the quality of the approximation may suffer. Even worse, the method may not be aware that it has failed.
 
 ```@example mode
-f = x -> tanh(400 * (x - 0.2))
+f(x) = tanh(400 * (x - 0.2))
 r = approximate(f, x)
 println("nearest pole is $(minimum(dist(z, I) for z in poles(r))) away")
 _, err = check(r);
@@ -41,7 +42,7 @@ println("max error on finer test points: ", err)
 In the continuous mode, the adaptive sampling of the domain attempts to ensure that the approximation is accurate everywhere.
 
 ```@example mode
-r = approximate(f, I; tol=1e-12)
+r = approximate(f; tol=1e-12)
 err = maximum(abs(f(x)- r(x)) for x in range(-1, 1, 3000))
 println("max error on finer test points: ", err)
 ```