Refactor Remez exchange into structured Julia and reduce allocations

benchmark/README.md

@@ -0,0 +1,60 @@
+# Remez rewrite validation
+
+Run from a checkout containing the original revision:
+
+```sh
+julia --startup-file=no benchmark/remez.jl
+julia --startup-file=no benchmark/remez_compatibility.jl
+julia --project -e 'using Pkg; Pkg.test(; julia_args=["--check-bounds=yes"])'
+```
+
+The first two scripts use only standard libraries. Their optional first argument
+selects a baseline git revision; the default is
+`7c798756cca39251da14bb5d38146feb2cdd1717`. They load the original source and the
+working-tree source in separate modules in the same Julia process. They neither
+change the checkout nor install benchmark dependencies.
+
+The timing script checks exact coefficient equality for each case, warms up both
+implementations, and reports the median of 31 samples of five calls. Sample order
+alternates between old/new and new/old. Allocated bytes include the returned
+coefficient vector and are measured after compilation. Compilation time and
+process resident memory are not measured. Timing results depend on the CPU and
+Julia version; they are not assertions in the package tests.
+
+The differential script compares coefficients bit for bit, warning messages, and
+exception types/messages on 600 seeded random designs. It varies tap count,
+symmetry, number of bands, weights, grid density, partial/full frequency coverage,
+and iteration limits. Both scripts require the baseline revision in git history;
+the ordinary package tests do not.
+
+## Measured results
+
+Julia 1.12.7, `znver3`, one Julia thread. Microseconds per call:
+
+| Design | Original | Rewrite | Original bytes | Rewrite bytes |
+| --- | ---: | ---: | ---: | ---: |
+| 15-tap low-pass | 13.76 | 11.45 | 5,152 | 3,776 |
+| 151-tap low-pass | 863.71 | 620.49 | 45,656 | 34,328 |
+| 152-tap weighted low-pass | 1,297.99 | 922.68 | 45,592 | 34,328 |
+| 51-tap high-pass | 59.44 | 42.47 | 16,080 | 11,976 |
+| 180-tap band-pass | 2,430.06 | 1,771.71 | 53,064 | 39,848 |
+| 20-tap Hilbert | 10.71 | 9.20 | 6,304 | 4,736 |
+| 21-tap Hilbert | 9.35 | 8.53 | 6,304 | 4,736 |
+| 200-tap differentiator | 954.59 | 715.02 | 59,808 | 45,272 |
+| 201-tap differentiator | 1,001.03 | 770.01 | 55,712 | 42,200 |
+| 63-tap partial-band design | 78.45 | 56.03 | 15,432 | 11,400 |
+| 63-tap low-pass, density 64 | 270.58 | 186.17 | 65,512 | 49,032 |
+| 511-tap low-pass | 5,852.15 | 4,957.51 | 151,888 | 114,768 |
+| 201-tap inverse sinc | 1,240.34 | 929.87 | 59,024 | 44,808 |
+| Three-argument weighted low-pass | 1,311.01 | 934.69 | 45,944 | 34,680 |
+| Three-argument Hilbert, 20 taps | 11.42 | 9.76 | 6,672 | 5,104 |
+| Three-argument Hilbert, 21 taps | 9.75 | 8.51 | 6,672 | 5,104 |
+| Three-argument differentiator, 200 taps | 964.16 | 729.10 | 60,400 | 45,816 |
+| Three-argument differentiator, 201 taps | 1,041.57 | 782.92 | 56,304 | 42,744 |
+
+These cases take approximately 9–31% less time and allocate 24–27% fewer bytes.
+The rewrite fills the frequency grid directly and reuses interpolation workspace
+for the partial-band basis conversion. Its barycentric evaluator groups four
+independent divisions while preserving accumulation order with fused operations.
+The exchange step keeps the original search order, endpoint margin, and stopping
+rules, including the returned approximation when `maxiter` is reached.

benchmark/remez.jl

@@ -0,0 +1,90 @@
+# julia --startup-file=no benchmark/remez.jl [baseline git revision]
+# No non-stdlib dependencies; both implementations run in the same process.
+using Printf
+
+module BaselineRemez
+    revision = isempty(ARGS) ? "7c798756cca39251da14bb5d38146feb2cdd1717" : ARGS[1]
+    source = read(`git -C $(@__DIR__) show $revision:src/Filters/remez_fir.jl`, String)
+    include_string(@__MODULE__, source, "baseline_remez.jl")
+end
+
+module CurrentRemez
+    include("../src/Filters/remez_fir.jl")
+end
+
+const REMEZ_CASES = [
+    ("small LP", 15, [(0., .2) => 1., (.3, .5) => 0.], (;)),
+    ("LP odd", 151, [(0., .475) => 1., (.5, 1.) => 0.], (; Hz=2.)),
+    ("LP even weighted", 152, [(0., .475) => (1., 1.), (.5, 1.) => (0., 2.)], (; Hz=2.)),
+    ("HP", 51, [(0., .75) => 0., (.8, 1.) => 1.], (; Hz=2.)),
+    ("BP", 180, [(0., .375) => 0., (.4, .5) => 1., (.525, 1.) => 0.], (; Hz=2., maxiter=30)),
+    ("Hilbert even", 20, [(.1, .95) => 1.], (; Hz=2., neg=true)),
+    ("Hilbert odd", 21, [(.1, .95) => 1.], (; Hz=2., neg=true)),
+    ("diff even", 200, [(.01, .99) => (f -> f / 2, f -> 1 / f)], (; Hz=2., neg=true)),
+    ("diff odd", 201, [(.05, .95) => (f -> f / 2, f -> 1 / f)], (; Hz=2., neg=true)),
+    ("partial band", 63, [(.05, .2) => 1., (.25, .45) => 0.], (;)),
+    ("dense grid", 63, [(0., .2) => 1., (.23, .5) => 0.], (; grid_density=64)),
+    ("long LP", 511, [(0., .2) => 1., (.205, .5) => 0.], (; maxiter=100)),
+    ("inverse sinc", 201, [(0., 2880.) => (f -> inv(sinc(f / 4800)), 1.),
+                          (10000., 153600.) => (0., 100.)], (; Hz=307200.)),
+]
+
+function elapsed_batch(f, args, kwargs, repetitions)
+    start = time_ns()
+    for _ in 1:repetitions
+        f(args...; kwargs...)
+    end
+    return (time_ns() - start) / repetitions
+end
+
+function compare_case(name, args, kwargs; newkwargs=kwargs, samples=31, repetitions=5)
+    old, new = BaselineRemez.remez, CurrentRemez.remez
+    reference = old(args...; kwargs...)
+    @assert isequal(reference, new(args...; newkwargs...)) name
+    # Compile both the filters and the timing harness before collecting samples.
+    elapsed_batch(old, args, kwargs, 2)
+    elapsed_batch(new, args, newkwargs, 2)
+    oldbytes = @allocated old(args...; kwargs...)
+    newbytes = @allocated new(args...; newkwargs...)
+    oldtimes, newtimes = Float64[], Float64[]
+    GC.gc()
+    for i in 1:samples
+        # Alternate order to reduce systematic thermal / frequency bias.
+        if isodd(i)
+            push!(oldtimes, elapsed_batch(old, args, kwargs, repetitions))
+            push!(newtimes, elapsed_batch(new, args, newkwargs, repetitions))
+        else
+            push!(newtimes, elapsed_batch(new, args, newkwargs, repetitions))
+            push!(oldtimes, elapsed_batch(old, args, kwargs, repetitions))
+        end
+    end
+    oldtime = sort!(oldtimes)[cld(samples, 2)] / 1000
+    newtime = sort!(newtimes)[cld(samples, 2)] / 1000
+    @printf("%-24s %10.2f %10.2f %7.3f %10d %10d\n",
+            name, oldtime, newtime, newtime / oldtime, oldbytes, newbytes)
+    return (; oldtime, newtime, oldbytes, newbytes)
+end
+
+function main()
+    println("Julia ", VERSION, " / ", Sys.CPU_NAME, " / threads=", Threads.nthreads())
+    println("Times: median microseconds per call; bytes: allocations per call, after warmup.")
+    @printf("%-24s %10s %10s %7s %10s %10s\n",
+            "case", "old μs", "new μs", "new/old", "old bytes", "new bytes")
+    for (name, taps, bands, kwargs) in REMEZ_CASES
+        compare_case(name, (taps, bands), kwargs)
+    end
+    compare_case("three-argument weighted", (152, [0., .475, .5, 1.], [1., 0.]),
+                 (; weight=[1., 2.], Hz=2.))
+    for n in (20, 21)
+        compare_case("three-arg Hilbert $n", (n, [.1, .95], [1.]),
+                     (; Hz=2., filter_type=BaselineRemez.filter_type_hilbert);
+                     newkwargs=(; Hz=2., filter_type=CurrentRemez.filter_type_hilbert))
+    end
+    for (n, edges) in ((200, [.01, .99]), (201, [.05, .95]))
+        compare_case("three-arg diff $n", (n, edges, [1.]),
+                     (; Hz=2., filter_type=BaselineRemez.filter_type_differentiator);
+                     newkwargs=(; Hz=2., filter_type=CurrentRemez.filter_type_differentiator))
+    end
+end
+
+abspath(PROGRAM_FILE) == (@__FILE__) && main()

benchmark/remez_compatibility.jl

@@ -0,0 +1,36 @@
+# julia --startup-file=no benchmark/remez_compatibility.jl [baseline git revision]
+# Keep this differential check outside the package tests: it requires git history.
+using Test, Random, Logging
+include("remez.jl")
+
+function outcome(f, taps, bands, kwargs)
+    logger = Test.TestLogger()
+    result = with_logger(logger) do
+        try
+            f(taps, bands; kwargs...)
+        catch e
+            (typeof(e), sprint(showerror, e))
+        end
+    end
+    return result, [(entry.level, entry.message) for entry in logger.logs]
+end
+
+@testset "Remez differential compatibility" begin
+    rng = MersenneTwister(9138)
+    for trial in 1:600
+        taps = rand(rng, 6:201)
+        nbands = rand(rng, 1:4)
+        edges = sort!(rand(rng, 2nbands) ./ 2)
+        if isodd(trial)
+            edges[1], edges[end] = 0., .5
+        end
+        bands = [(edges[2i-1], edges[2i]) => (rand(rng), 10.0^rand(rng, -1:1))
+                 for i in 1:nbands]
+        kwargs = (; neg=rand(rng, Bool), grid_density=rand(rng, (4, 16, 32)),
+                   maxiter=rand(rng, (1, 2, 10, 25, 40)))
+        @testset "case $trial: taps=$taps, $kwargs" begin
+            @test isequal(outcome(BaselineRemez.remez, taps, bands, kwargs),
+                          outcome(CurrentRemez.remez, taps, bands, kwargs))
+        end
+    end
+end