From 3956cea3bb3b5b28e97e81d743705a0d4453ccdc Mon Sep 17 00:00:00 2001 From: "Luis M. B. Varona" Date: Thu, 29 Jan 2026 03:29:32 -0400 Subject: [PATCH] Switch back from @repl to jldoctest blocks in docstrings Since Documenter.jl stopped executing @repl blocks in the generated documentation, we switch back from @repl to jldoctest blocks in the docstrings so that users can still see the expected output in the REPL and in the source code. --- CHANGELOG.md | 1 + .../Exact/solvers/brute_force_search.jl | 25 +- .../Exact/solvers/caprara_salazar_gonzalez.jl | 34 +- .../Exact/solvers/del_corso_manzini.jl | 144 ++++++-- .../Exact/solvers/saxe_gurari_sudborough.jl | 70 +++- .../Heuristic/solvers/cuthill_mckee.jl | 338 ++++++++++++++---- .../solvers/gibbs_poole_stockmeyer.jl | 169 +++++++-- src/Minimization/core.jl | 162 +++++++-- src/Recognition/core.jl | 95 ++++- .../deciders/brute_force_search.jl | 37 +- .../deciders/caprara_salazar_gonzalez.jl | 37 +- src/Recognition/deciders/del_corso_manzini.jl | 95 +++-- .../deciders/saxe_gurari_sudborough.jl | 37 +- src/core.jl | 127 +++++-- src/utils.jl | 233 +++++++++--- 15 files changed, 1254 insertions(+), 350 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cd57c89..74aa624 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), ### Changed +- Switched back from `@repl` blocks to `jldoctest` blocks in the docstrings, since *Documenter.jl* stopped executing `@repl` blocks in the generated documentation (#217). - Switched from recomputing Del Corso–Manzini placement deadlines from scratch each step to maintaining them incrementally as nodes are placed, incorporating the last missing optimization from the paper (#216). - Changed the default recognition decider from Del Corso–Manzini to Caprara–Salazar-González, which is now considerably faster after the performance fixes (#214). - Updated the `_blb_connected` helper (used in `bandwidth_lower_bound`) to avoid unnecessary allocations, now requiring only `O(n)` auxiliary space instead of `O(n^2)` (#213). diff --git a/src/Minimization/Exact/solvers/brute_force_search.jl b/src/Minimization/Exact/solvers/brute_force_search.jl index 8fb5ccd..26f8ca4 100644 --- a/src/Minimization/Exact/solvers/brute_force_search.jl +++ b/src/Minimization/Exact/solvers/brute_force_search.jl @@ -33,13 +33,24 @@ is not *exactly* ``Θ(n²)``, although it is close.) The algorithm always iterates over all possible permutations, so it is infeasible to go above ``9×9`` or ``10×10`` without incurring multiple-hour runtimes. Nevertheless, we see that it is quite effective for, say, ``8×8``: -```@repl -using Random, SparseArrays -Random.seed!(628318); -(n, p) = (8, 0.2); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -minimize_bandwidth(A, Minimization.BruteForceSearch()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(628318); + +julia> (n, p) = (8, 0.2); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 3 + * Original Bandwidth: 6 + * Matrix Size: 8×8 ``` # Notes diff --git a/src/Minimization/Exact/solvers/caprara_salazar_gonzalez.jl b/src/Minimization/Exact/solvers/caprara_salazar_gonzalez.jl index a58cc84..aa5ee84 100644 --- a/src/Minimization/Exact/solvers/caprara_salazar_gonzalez.jl +++ b/src/Minimization/Exact/solvers/caprara_salazar_gonzalez.jl @@ -52,14 +52,32 @@ exponential growth in time complexity with respect to ``n``. # Examples We verify the optimality of the ordering found by Caprara–Salazar-González for a random ``8×8`` matrix via a brute-force search over all possible permutations up to reversal: -```@repl -using Random, SparseArrays -Random.seed!(5883); -(n, p) = (8, 0.25); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) -res_csg = minimize_bandwidth(A, Minimization.CapraraSalazarGonzalez()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(5883); + +julia> (n, p) = (8, 0.25); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 4 + * Original Bandwidth: 7 + * Matrix Size: 8×8 + +julia> res_csg = minimize_bandwidth(A, Minimization.CapraraSalazarGonzalez()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Caprara–Salazar-González + * Approach: exact + * Minimum Bandwidth: 4 + * Original Bandwidth: 7 + * Matrix Size: 8×8 ``` # Notes diff --git a/src/Minimization/Exact/solvers/del_corso_manzini.jl b/src/Minimization/Exact/solvers/del_corso_manzini.jl index a24110a..53ff419 100644 --- a/src/Minimization/Exact/solvers/del_corso_manzini.jl +++ b/src/Minimization/Exact/solvers/del_corso_manzini.jl @@ -54,14 +54,32 @@ Based on experimental results, the algorithm is feasible for ``n×n`` matrices u # Examples We verify the optimality of the ordering found by Del Corso–Manzini for a random ``9×9`` matrix via a brute-force search over all possible permutations up to reversal: -```@repl -using Random, SparseArrays -Random.seed!(0117); -(n, p) = (9, 0.5); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) -res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(0117); + +julia> (n, p) = (9, 0.5); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 5 + * Original Bandwidth: 8 + * Matrix Size: 9×9 + +julia> res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Del Corso–Manzini + * Approach: exact + * Minimum Bandwidth: 5 + * Original Bandwidth: 8 + * Matrix Size: 9×9 ``` We now generate (and shuffle) a random ``40×40`` matrix with minimum bandwidth ``10`` using @@ -71,16 +89,32 @@ cases, `random_banded_matrix(n, k)` *does* generate matrices with minimum bandwi Nevertheless, this example demonstrates that Del Corso–Manzini at the very least finds a quite good ordering, even though exact optimality—which *is* guaranteed by the original paper [DM99]—is not explicitly verified.) -```@repl -using Random -Random.seed!(0201); -(n, k) = (40, 10); -A = MatrixBandwidth.random_banded_matrix(n, k); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -minimize_bandwidth(A_shuffled, Minimization.DelCorsoManzini()) +```jldoctest +julia> using Random + +julia> Random.seed!(0201); + +julia> (n, k) = (40, 10); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +10 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +36 + +julia> minimize_bandwidth(A_shuffled, Minimization.DelCorsoManzini()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Del Corso–Manzini + * Approach: exact + * Minimum Bandwidth: 10 + * Original Bandwidth: 36 + * Matrix Size: 40×40 ``` # Notes @@ -184,14 +218,32 @@ for a random ``9×9`` matrix via a brute-force search over all possible permutat reversal. The depth parameter is not explicitly set; instead, some near-optimal value is automatically computed upon the first [`MatrixBandwidth.Minimization.minimize_bandwidth`](@ref) function call. -```@repl -using Random, SparseArrays -Random.seed!(548836); -(n, p) = (9, 0.2); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) -res_dcm_ps = minimize_bandwidth(A, Minimization.DelCorsoManziniWithPS()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(548836); + +julia> (n, p) = (9, 0.2); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 3 + * Original Bandwidth: 8 + * Matrix Size: 9×9 + +julia> res_dcm_ps = minimize_bandwidth(A, Minimization.DelCorsoManziniWithPS()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Del Corso–Manzini with perimeter search + * Approach: exact + * Minimum Bandwidth: 3 + * Original Bandwidth: 8 + * Matrix Size: 9×9 ``` We now generate (and shuffle) a random ``30×30`` matrix with minimum bandwidth ``8`` using @@ -201,18 +253,34 @@ finds a bandwidth-``8`` ordering, which is (we claim) optimal up to symmetric pe `< k`. Nevertheless, this example demonstrates that Del Corso–Manzini at the very least finds a quite good ordering, even though exact optimality—which *is* guaranteed by the original paper [DM99]—is not explicitly verified.) In this case, we set the depth parameter -to ``4`` beforehand instead of relying on [`Recognition.dcm_ps_optimal_depth`](@ref). - -```@repl -using Random -Random.seed!(78779); -(n, k, depth) = (30, 8, 4); -A = MatrixBandwidth.random_banded_matrix(n, k); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -minimize_bandwidth(A_shuffled, Minimization.DelCorsoManziniWithPS(depth)) +to ``3`` beforehand instead of relying on [`Recognition.dcm_ps_optimal_depth`](@ref). + +```jldoctest +julia> using Random + +julia> Random.seed!(78779); + +julia> (n, k, depth) = (30, 8, 3); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +8 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +25 + +julia> minimize_bandwidth(A_shuffled, Minimization.DelCorsoManziniWithPS(depth)) +Results of Bandwidth Minimization Algorithm + * Algorithm: Del Corso–Manzini with perimeter search + * Approach: exact + * Minimum Bandwidth: 8 + * Original Bandwidth: 25 + * Matrix Size: 30×30 ``` # Notes diff --git a/src/Minimization/Exact/solvers/saxe_gurari_sudborough.jl b/src/Minimization/Exact/solvers/saxe_gurari_sudborough.jl index d902e1f..bfc1d9f 100644 --- a/src/Minimization/Exact/solvers/saxe_gurari_sudborough.jl +++ b/src/Minimization/Exact/solvers/saxe_gurari_sudborough.jl @@ -56,14 +56,32 @@ aggressive pruning strategies keep their effective search space very small in pr # Examples We verify the optimality of the ordering found by Saxe–Gurari–Sudborough for a random ``9×9`` matrix via a brute-force search over all possible permutations up to reversal: -```@repl -using Random, SparseArrays -Random.seed!(52452); -(n, p) = (9, 0.5); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) -res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(52452); + +julia> (n, p) = (9, 0.5); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> res_bf = minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 5 + * Original Bandwidth: 8 + * Matrix Size: 9×9 + +julia> res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Approach: exact + * Minimum Bandwidth: 5 + * Original Bandwidth: 8 + * Matrix Size: 9×9 ``` We now generate (and shuffle) a random ``25×25`` matrix with minimum bandwidth ``5`` using @@ -73,16 +91,32 @@ cases, `random_banded_matrix(n, k)` generates matrices with minimum bandwidth `< appears to be one such case. Although we do not explicitly verify exact optimality—which *is* guaranteed by the original paper [GS84]—here via brute-force search, this example demonstrates that Saxe–Gurari–Sudborough at the very least finds a quite good ordering.) -```@repl -using Random -Random.seed!(937497); -(n, k, p) = (25, 5, 0.25); -A = MatrixBandwidth.random_banded_matrix(n, k; p=p); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -minimize_bandwidth(A_shuffled, Minimization.SaxeGurariSudborough()) +```jldoctest +julia> using Random + +julia> Random.seed!(937497); + +julia> (n, k, p) = (25, 5, 0.25); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k; p=p); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +5 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +19 + +julia> minimize_bandwidth(A_shuffled, Minimization.SaxeGurariSudborough()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Approach: exact + * Minimum Bandwidth: 4 + * Original Bandwidth: 19 + * Matrix Size: 25×25 ``` # Notes diff --git a/src/Minimization/Heuristic/solvers/cuthill_mckee.jl b/src/Minimization/Heuristic/solvers/cuthill_mckee.jl index 6cf1c5f..b4e958f 100644 --- a/src/Minimization/Heuristic/solvers/cuthill_mckee.jl +++ b/src/Minimization/Heuristic/solvers/cuthill_mckee.jl @@ -52,44 +52,143 @@ it hard to verify whether Cuthill–McKee finds a truly optimal ordering or simp near-optimal one. Nevertheless, the results are still very good in practice. Cuthill–McKee finds a good ordering for a ``30×30`` matrix: -```@repl -using Random -Random.seed!(13); -(n, k) = (30, 5); -A = MatrixBandwidth.random_banded_matrix(n, k); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee()) +```jldoctest +julia> using Random + +julia> Random.seed!(13); + +julia> (n, k) = (30, 5); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +5 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +25 + +julia> minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Cuthill–McKee + * Approach: heuristic + * Minimum Bandwidth: 8 + * Original Bandwidth: 25 + * Matrix Size: 30×30 ``` -Cuthill–McKee finds a good ordering for a structurally symmetric ``183×183`` matrix with +Cuthill–McKee finds a good ordering for a structurally symmetric ``276×276`` matrix with multiple (separate) connected components: -```@repl -using Random, SparseArrays -Random.seed!(37452); -(max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 7); -components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); - -for i in 1:num_ccs # Some components may themselves be disconnected - cc_size = rand(1:max_cc_size); - cc_band = rand(0:min(max_band, cc_size - 1)); - components[i] = sparse( - MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) - ); -end - -A = blockdiag(components...); # `A` has least 7 connected components -perm = randperm(sum(map(cc -> size(cc, 1), components))); -A_shuffled = A[perm, perm]; -res = minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee()); -A # The original matrix -A_shuffled # A far-from-optimal ordering of `A` -A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -res # Even better than the original bandwidth (which was, clearly, not yet optimal) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(37452); + +julia> (max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 7); + +julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); + +julia> for i in 1:num_ccs # Some components may themselves be disconnected + cc_size = rand(1:max_cc_size); + cc_band = rand(0:min(max_band, cc_size - 1)); + components[i] = sparse( + MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) + ); + end + +julia> A = blockdiag(components...); # `A` has least 7 connected components + +julia> perm = randperm(sum(map(cc -> size(cc, 1), components))); + +julia> A_shuffled = A[perm, perm]; + +julia> res = minimize_bandwidth(A_shuffled, Minimization.CuthillMcKee()); + +julia> A # The original matrix +276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries: +⎡⢾⡷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠘⢻⣲⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠘⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠈⠿⡧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠉⢯⡷⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠚⣤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢻⣶⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠯⡧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⣤⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠛⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣢⡀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡢⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢴⣷⡀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢿⡷⎦ + +julia> A_shuffled # A far-from-optimal ordering of `A` +276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries: +⎡⠁⢄⠀⢀⠀⠀⠀⢀⠠⠀⠀⠐⠀⠀⠀⠐⢀⡐⠀⠀⠀⢀⠀⠀⠀⠀⠐⠀⢠⠀⠀⠀⡄⠀⠀⠐⠀⠀⠂⠄⎤ +⎢⠀⢀⠱⠂⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⢨⠀⠀⠀⠀⠀⡀⠁⠠⠀⠘⠀⠀⠡⢀⠈⠀⠀⠀⠀⠀⠀⠄⠀⠁⠁⎥ +⎢⠀⠀⠀⠀⠑⢀⠀⠂⠀⠀⠀⠀⢐⠀⠀⠠⠈⠠⠀⠀⠀⠐⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⢢⢀⢀⠀⎥ +⎢⠀⢀⡀⠀⠠⠀⠁⠄⠀⠠⠀⠄⠀⠀⠀⠄⠀⠀⠀⠀⢀⠀⠀⢀⠀⠑⠀⠀⠐⠠⠀⠀⠠⠨⠂⠀⠀⠀⠀⠀⎥ +⎢⠀⠂⠀⠀⠀⠀⠀⡀⠱⢆⡀⠂⠀⠀⠀⠀⠀⠀⢀⢊⠀⠐⠐⠈⠀⠈⠀⢀⠄⠀⡀⠀⢁⢀⠠⠀⠃⠀⠊⠀⎥ +⎢⢀⠀⠀⠀⠀⠀⠀⠄⠠⠈⠑⠀⢀⠐⠀⠌⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⢀⠉⢀⠀⠠⠈⠀⠀⣁⠁⎥ +⎢⠀⠀⠀⠀⠐⠐⠀⠀⠀⠀⢀⠐⠁⠄⠈⠀⢌⠀⠆⠠⢀⠀⠄⠐⠰⠀⠀⠀⠁⠰⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⢀⠀⠂⠒⠀⡀⠀⠄⠀⠀⡀⠄⠂⠀⠐⢄⠁⢀⠀⠀⠀⡀⠀⠀⠀⠀⡠⠀⠀⠀⠀⠀⠀⠀⠀⢈⠀⠀⠀⠁⎥ +⎢⢀⠰⠀⠀⠂⡀⠀⠀⠀⠀⠀⠈⠂⠑⠁⢀⠐⠄⠄⠂⠂⠜⠄⠀⠀⠀⡄⠀⠀⢀⠀⠠⠀⢀⠄⠀⢀⠀⠂⡂⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⡠⢐⠀⠀⠈⡁⠀⠀⠠⠁⠀⢀⠀⠀⠀⠀⡀⠀⠀⢀⠀⠈⠃⠀⠸⠈⠠⠀⠀⠀⢄⠂⎥ +⎢⠀⢀⠄⠈⢀⠀⠀⠐⢀⠀⠀⠀⠀⠐⠀⠠⣈⠄⠀⠀⠐⢀⠀⡀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠊⠀⠠⠀⠐⠀⠀⎥ +⎢⠀⠀⠀⠂⢀⠀⠀⢀⡐⠀⠀⠀⢀⠁⠀⠀⠀⠁⠀⠀⠀⠠⠄⣥⠉⠈⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⡀⠀⠀⠀⎥ +⎢⠀⠀⠒⠀⠀⠀⢄⠀⡀⠀⠀⠀⠐⠂⠀⠀⠀⠀⠀⠈⠀⠀⡃⠀⠁⢀⠀⠀⢀⡀⢈⠈⠀⠀⠀⠂⠀⠠⠂⠂⎥ +⎢⠐⠀⠄⡀⠀⠀⠀⠀⠀⢀⠀⠈⠀⠀⠀⠊⠀⠉⠀⢀⠀⠀⠀⠀⠀⠀⠑⠀⠀⠀⠀⠀⢀⠀⠈⠀⠛⠃⢄⠀⎥ +⎢⠀⠒⡀⠐⠀⠀⠐⡀⠀⠁⠀⠀⢁⡀⠀⠀⠀⢀⡀⠀⠀⠀⠀⠀⠀⠰⠀⠀⠀⢄⠀⠰⠀⠠⠠⢀⠀⠀⢂⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⡄⠐⠂⠀⠀⠀⠀⡀⠉⠀⠐⠀⠀⠈⡂⠐⠀⠀⢀⡀⠀⣠⠀⠄⠠⠀⠀⡀⠀⠀⎥ +⎢⠀⠉⠀⠀⠀⠀⡀⡂⠁⢐⠀⠐⠀⠀⠀⠀⠀⢀⡒⠂⡠⠀⠀⠀⠀⠀⠀⠐⠀⡀⠀⠄⠑⠄⠀⠀⠀⠀⠀⠀⎥ +⎢⢀⠀⠀⠀⠀⡀⠈⠀⠀⠂⡀⠂⠀⠀⡀⢀⠀⠁⠀⠂⠀⡀⠀⠀⠠⠀⠂⠀⠀⢂⠀⠂⠀⠀⠁⢀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠁⠈⢒⠀⠀⠉⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⢀⠀⠀⠈⠀⡀⠿⠀⠀⠀⠀⠠⠀⠀⠀⠀⠱⠆⠀⠀⎥ +⎣⠈⠄⠅⠀⠀⠐⠀⠀⠊⠀⠅⠘⠀⠀⠄⠀⠨⠠⠠⠑⠀⠀⠀⠀⠨⠀⠀⠑⠈⠐⠀⠀⠀⠀⠀⠀⠀⠀⠔⢅⎦ + +julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` +276×276 SparseMatrixCSC{Float64, Int64} with 464 stored entries: +⎡⠱⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠉⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠘⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⡦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠺⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠚⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⣤⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⢄⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢄⠀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⎦ + +julia> bandwidth(A) +7 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +266 + +julia> res # Even better than the original bandwidth (which was, clearly, not yet optimal) +Results of Bandwidth Minimization Algorithm + * Algorithm: Cuthill–McKee + * Approach: heuristic + * Minimum Bandwidth: 5 + * Original Bandwidth: 266 + * Matrix Size: 276×276 ``` # Notes @@ -171,44 +270,143 @@ it hard to verify whether reverse Cuthill–McKee finds a truly optimal ordering near-optimal one. Nevertheless, the results are still very good in practice. Reverse Cuthill–McKee finds a good ordering for a ``35×35`` matrix: -```@repl -using Random -Random.seed!(87); -(n, k) = (35, 3); -A = MatrixBandwidth.random_banded_matrix(n, k); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee()) +```jldoctest +julia> using Random + +julia> Random.seed!(87); + +julia> (n, k) = (35, 3); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +3 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +30 + +julia> minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Reverse Cuthill–McKee + * Approach: heuristic + * Minimum Bandwidth: 3 + * Original Bandwidth: 30 + * Matrix Size: 35×35 ``` -Reverse Cuthill–McKee finds a good ordering for a ``235×235`` matrix with multiple +Reverse Cuthill–McKee finds a good ordering for a ``233×233`` matrix with multiple (separate) connected components: -```@repl -using Random, SparseArrays -Random.seed!(5747); -(max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 8); -components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); - -for i in 1:num_ccs # Some components may themselves be disconnected - cc_size = rand(1:max_cc_size); - cc_band = rand(0:min(max_band, cc_size - 1)); - components[i] = sparse( - MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) - ); -end - -A = blockdiag(components...); # `A` has least 8 connected components -perm = randperm(sum(map(cc -> size(cc, 1), components))); -A_shuffled = A[perm, perm]; -res = minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee()); -A # The original matrix -A_shuffled # A far-from-optimal ordering of `A` -A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -res # Gets very close to the original bandwidth +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(5747); + +julia> (max_cc_size, max_band, p, num_ccs) = (60, 9, 0.2, 8); + +julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); + +julia> for i in 1:num_ccs # Some components may themselves be disconnected + cc_size = rand(1:max_cc_size); + cc_band = rand(0:min(max_band, cc_size - 1)); + components[i] = sparse( + MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) + ); + end + +julia> A = blockdiag(components...); # `A` has least 8 connected components + +julia> perm = randperm(sum(map(cc -> size(cc, 1), components))); + +julia> A_shuffled = A[perm, perm]; + +julia> res = minimize_bandwidth(A_shuffled, Minimization.ReverseCuthillMcKee()); + +julia> A # The original matrix +233×233 SparseMatrixCSC{Float64, Int64} with 571 stored entries: +⎡⢾⣳⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠘⢿⡷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠘⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠘⢴⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠙⠻⢂⣀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⣮⣿⣲⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠺⢿⡷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠰⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢡⣶⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠊⠀⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠛⡄⡭⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣻⣾⡄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠩⢿⣷⣆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣯⡿⡆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠩⠪⣢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢯⣳⡄⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢿⣷⣀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠚⡠⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⣠⠀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⎦ + +julia> A_shuffled # A far-from-optimal ordering of `A` +233×233 SparseMatrixCSC{Float64, Int64} with 571 stored entries: +⎡⠑⠀⠀⠀⠀⠀⡄⠄⠀⠐⠈⠘⠀⠀⠐⠀⠁⠀⠀⠐⠈⢀⡀⠁⠄⠀⡁⠀⢀⠀⠐⢃⢠⠀⡁⠀⡀⠀⠀⠀⎤ +⎢⠀⠀⠁⢄⠀⠐⠀⠀⠈⠀⡐⠀⠈⡀⠀⠡⡀⠀⢀⠀⢀⠔⡂⠀⠀⠂⡀⢂⠅⠀⠀⠀⠀⢀⠀⠀⠄⠀⡀⠀⎥ +⎢⠀⠀⢀⠀⠀⠀⢠⠀⠀⡀⠀⡀⠄⠀⠄⠀⠠⡀⠠⠄⠀⢠⠈⠀⠐⠀⠀⢒⠀⠀⠁⠀⡀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠍⠀⠀⠀⠒⠐⢄⠀⠀⢀⠠⠄⠀⡠⠀⠀⢂⠂⠀⠩⠘⠀⢠⠠⠀⠄⠐⠀⠀⢀⡀⠀⠀⢀⠀⠀⠀⠀⠂⎥ +⎢⢀⠀⠂⠀⠀⠠⠀⠀⠑⢀⠀⠀⠠⠀⠄⠀⢀⠀⠀⠈⠀⢈⢐⠀⠀⠠⠀⠠⠀⠀⣐⠀⠀⡀⡀⠀⠀⢀⠈⠣⎥ +⎢⣂⠀⠐⠈⠀⠠⠀⡐⠀⠀⢁⢔⠠⠀⡁⠀⠀⠀⠐⠂⠀⠀⠀⠄⠄⠈⠀⠠⢀⠈⠂⠀⢀⠀⠀⠀⢀⠀⠈⠊⎥ +⎢⠀⠀⠂⠠⠀⠁⠀⠁⠀⠂⠀⠂⠀⢀⠂⠀⠀⠐⠁⠀⠠⠀⡀⠡⠀⠁⠀⠀⠀⠀⠀⠀⢂⢂⠐⠀⠄⠀⠀⠀⎥ +⎢⠐⠀⠄⡀⠀⠁⠀⠊⠀⠁⠁⠈⠈⠀⠁⠀⠀⠈⠁⠀⠈⠁⠄⠀⠀⠄⠠⠁⡀⠀⠨⢁⠁⡀⠈⠤⠀⠀⠄⠢⎥ +⎢⠁⠀⠀⠈⠀⠢⠠⢀⠀⠐⠀⠀⢀⠀⡀⠀⠊⠀⠠⠀⠀⠰⠀⢀⠀⠂⠀⠐⠀⠁⠀⠈⠀⠀⠂⠀⠄⠀⠀⠀⎥ +⎢⢀⠀⠀⠐⠀⠆⠈⠀⡀⠀⠰⠀⠁⠀⠁⠀⠀⠂⡐⢈⠠⢀⠀⠁⠀⠐⠀⠐⠠⠀⠌⠄⠂⢀⡀⠠⠀⢀⢄⠐⎥ +⎢⠂⢀⢀⠔⠀⣀⣃⠂⡀⢀⠀⠀⠀⠂⠆⠀⢀⡀⠀⢂⠑⠀⡀⠂⠀⠀⠄⠀⢀⠀⠀⠀⠀⠀⣂⠀⡂⠀⠀⠀⎥ +⎢⠄⠈⠈⠈⠂⠀⠀⣀⠐⠐⠀⠄⠄⡈⠀⠁⠀⢀⠄⠀⠠⠈⠄⢅⡀⢀⠀⠢⠀⠄⠀⠀⠀⠂⠀⢀⠀⢀⠅⠀⎥ +⎢⠀⠁⠠⠀⠐⠀⠀⠂⠀⡀⡀⠁⠄⠀⠀⠄⠠⠀⢀⠀⠀⠀⠀⢈⠀⠄⡀⢀⠐⢀⠀⠖⠀⠀⠀⠀⠀⠀⢐⠀⎥ +⎢⠁⠈⠠⢈⢠⢀⢀⠁⠀⡀⠀⡀⠀⠀⠄⠂⢀⠀⢀⠀⠀⠁⠠⡀⠀⢈⢁⢔⠀⠈⠈⠂⡂⠀⠀⡢⠀⠀⠈⠈⎥ +⎢⠀⠐⠁⠁⠀⠀⠀⠀⠀⠀⡀⠐⠀⠀⠀⠈⠄⠀⠀⠂⠀⠐⠀⠄⠐⢀⡀⠀⠀⠀⠂⠀⠀⠰⡂⠄⠀⠠⢀⠀⎥ +⎢⠴⢀⠀⠀⠁⠀⠀⠰⠐⠘⠈⠀⠀⠀⠆⢂⡀⠀⠂⠅⠀⠀⠀⠀⢠⠄⠢⠀⠈⠀⠁⠀⠀⠀⠀⠊⠀⠈⠀⠀⎥ +⎢⠀⠒⠀⢀⠀⠈⠀⠀⠀⠠⠀⠐⠨⢐⠁⠠⠀⠀⠈⢀⠀⠀⠠⠀⠀⠀⠈⠈⢀⡀⠀⠀⠐⠄⠠⡀⠀⡀⠠⠈⎥ +⎢⠁⠈⠀⠀⠀⠀⠀⠐⠀⠈⠀⠀⠐⠀⠂⡄⠈⠀⠀⡈⠈⠘⠀⢀⠀⠀⠠⡠⠈⠌⡠⠀⠀⠢⠁⠀⠂⠀⠈⠀⎥ +⎢⠀⠈⠀⠁⠀⠀⠀⠀⠀⢀⠀⠐⠀⠁⠀⠀⠀⠁⠀⢀⠈⠈⠀⢀⠀⠀⠀⠀⠀⡀⡀⠀⠀⠠⠈⠀⠊⢀⠀⠀⎥ +⎣⠀⠀⠀⠈⠀⠀⠠⠀⠦⡀⡢⠀⠀⠀⠠⡁⠀⠀⢀⠑⠀⠀⠁⠁⠐⠐⡂⠀⠀⠐⠀⠀⡀⠂⠂⠀⠀⠀⠑⠄⎦ + +julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` +233×233 SparseMatrixCSC{Float64, Int64} with 571 stored entries: +⎡⠁⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠐⡠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠡⣢⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⡶⣠⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢺⣾⡻⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢻⣲⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢻⣲⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠿⣧⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠿⣣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢮⣷⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⢾⡷⢤⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣵⣿⣦⡀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣤⣿⡄⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢯⣷⡄⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⎦ + +julia> bandwidth(A) +9 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +202 + +julia> res # Gets very close to the original bandwidth +Results of Bandwidth Minimization Algorithm + * Algorithm: Reverse Cuthill–McKee + * Approach: heuristic + * Minimum Bandwidth: 11 + * Original Bandwidth: 202 + * Matrix Size: 233×233 ``` # Notes diff --git a/src/Minimization/Heuristic/solvers/gibbs_poole_stockmeyer.jl b/src/Minimization/Heuristic/solvers/gibbs_poole_stockmeyer.jl index 7df31a7..9429968 100644 --- a/src/Minimization/Heuristic/solvers/gibbs_poole_stockmeyer.jl +++ b/src/Minimization/Heuristic/solvers/gibbs_poole_stockmeyer.jl @@ -66,44 +66,143 @@ it hard to verify whether Gibbs–Poole–Stockmeyer finds a truly optimal order near-optimal one. Nevertheless, the results are still very good in practice. Gibbs–Poole–Stockmeyer finds a good ordering for a ``40×40`` matrix: -```@repl -using Random -Random.seed!(561); -(n, k) = (40, 7); -A = MatrixBandwidth.random_banded_matrix(n, k); -perm = randperm(n); -A_shuffled = A[perm, perm]; -bandwidth(A) -bandwidth(A_shuffled) -minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer()) -``` +```jldoctest +julia> using Random -Gibbs–Poole–Stockmeyer finds a good ordering for a ``748×748`` matrix with multiple -(separate) connected components: +julia> Random.seed!(561); + +julia> (n, k) = (40, 7); + +julia> A = MatrixBandwidth.random_banded_matrix(n, k); + +julia> perm = randperm(n); + +julia> A_shuffled = A[perm, perm]; + +julia> bandwidth(A) +7 + +julia> bandwidth(A_shuffled) +37 + +julia> minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Gibbs–Poole–Stockmeyer + * Approach: heuristic + * Minimum Bandwidth: 8 + * Original Bandwidth: 37 + * Matrix Size: 40×40 ``` -using Random, SparseArrays -Random.seed!(271828); -(max_cc_size, max_band, p, num_ccs) = (120, 13, 0.3, 11); -components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); - -for i in 1:num_ccs # Some components may themselves be disconnected - cc_size = rand(1:max_cc_size); - cc_band = rand(0:min(max_band, cc_size - 1)); - components[i] = sparse( - MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) - ); -end -A = blockdiag(components...); # `A` has least 8 connected components -perm = randperm(sum(map(cc -> size(cc, 1), components))); -A_shuffled = A[perm, perm]; -res = minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer()); -A # The original matrix -A_shuffled # A far-from-optimal ordering of `A` -A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` -bandwidth(A) -bandwidth(A_shuffled) # Much larger after shuffling -res # Gets very close to the original bandwidth +Gibbs–Poole–Stockmeyer finds a good ordering for a ``733×733`` matrix with multiple +(separate) connected components: +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(271828); + +julia> (max_cc_size, max_band, p, num_ccs) = (120, 13, 0.3, 11); + +julia> components = Vector{SparseMatrixCSC{Float64, Int64}}(undef, num_ccs); + +julia> for i in 1:num_ccs # Some components may themselves be disconnected + cc_size = rand(1:max_cc_size); + cc_band = rand(0:min(max_band, cc_size - 1)); + components[i] = sparse( + MatrixBandwidth.random_banded_matrix(cc_size, cc_band; p=p) + ); + end + +julia> A = blockdiag(components...); # `A` has least 8 connected components + +julia> perm = randperm(sum(map(cc -> size(cc, 1), components))); + +julia> A_shuffled = A[perm, perm]; + +julia> res = minimize_bandwidth(A_shuffled, Minimization.GibbsPooleStockmeyer()); + +julia> A # The original matrix +733×733 SparseMatrixCSC{Float64, Int64} with 2864 stored entries: +⎡⠿⣧⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠘⠻⣦⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠙⢻⣶⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⢆⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⣦⡀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦ + +julia> A_shuffled # A far-from-optimal ordering of `A` +733×733 SparseMatrixCSC{Float64, Int64} with 2864 stored entries: +⎡⠑⢄⠐⣆⡇⠂⡒⠆⠄⠅⠒⠜⡥⠠⢂⡒⠢⢈⢑⠐⠰⢤⠒⠀⢰⢠⠂⢢⡀⠪⠂⡀⢂⠔⡗⡔⠂⣴⠄⢀⎤ +⎢⠰⢤⠱⠆⡐⠃⠠⡀⢱⡊⣜⢖⡈⡐⢀⠀⢎⣈⡹⢠⣭⡜⡌⢄⠰⠂⠌⢧⢉⣂⢁⡥⠈⢘⡄⢕⣙⠬⡐⠓⎥ +⎢⠩⠉⠴⠈⠟⢅⢊⠴⠔⠎⠘⠿⢼⡈⡁⠚⠘⠲⠽⠉⡯⠐⡫⢛⡈⠀⠂⠥⠠⠘⢂⣅⡬⡢⢫⡀⢈⡌⡔⠀⎥ +⎢⠸⠌⠀⠢⢊⡔⠑⣤⠀⠃⠬⠈⠤⢠⠁⢐⡦⢒⢸⠄⠵⢧⢸⠉⠤⠂⠖⠀⣉⢭⠔⢦⠐⠖⢯⠠⠆⡈⠀⠆⎥ +⎢⠄⠅⡱⠲⡰⠅⠤⠀⡕⢍⠀⠀⣏⠘⡄⠰⢉⢡⢒⠊⠭⡧⠾⢸⣣⣔⡈⡄⢜⡯⠄⠁⡐⠤⠆⢎⢳⡀⠀⣀⎥ +⎢⣘⠄⢲⢝⣶⡄⡂⠃⠀⠀⢴⣷⡑⠪⡠⢆⢰⣨⢻⠀⣰⣬⠋⠔⣛⡈⠪⠄⢼⢵⠀⠅⢁⣤⠀⣥⣅⣳⣰⠄⎥ +⎢⠁⡋⢂⠨⡒⠳⠀⣃⣋⠙⡱⡈⡛⣬⠖⡄⢊⠸⢩⢽⡙⢏⠊⢔⠩⠋⢘⣀⠈⡕⡅⣠⠈⡋⠨⡓⣺⣫⠓⠀⎥ +⎢⢨⠰⠀⠐⣡⠈⢁⢀⢀⡉⠠⢎⠘⠥⠑⢄⡁⠈⢁⠀⠒⠺⠬⠋⠠⠒⠄⡀⠊⢠⠡⠠⠒⢌⢂⠀⢪⢫⠡⠂⎥ +⎢⡈⢂⡊⢱⢲⡀⢨⢋⠇⣐⡐⣲⣊⡐⡁⠈⠰⣦⠸⣢⠑⢛⡔⠁⢎⡽⡅⢠⠠⡖⢆⣂⢠⠀⢓⠤⡜⡮⠺⠔⎥ +⎢⢑⠐⠓⣊⡗⠃⠒⠖⡸⠐⠛⠒⣇⣖⠁⠐⠲⣢⢑⢔⡪⣈⡋⢀⣙⣔⠂⣒⠻⡜⠁⡇⢅⠕⠂⢰⣅⣐⢈⡋⎥ +⎢⠐⣆⣃⠿⢋⠋⠵⣇⠧⡧⡐⣾⡷⢌⣸⡀⣵⢀⡊⢪⠻⣦⠪⡌⢹⠁⠂⡔⣐⡨⠇⡆⠪⣋⡶⡅⡘⡡⠂⠑⎥ +⎢⠘⠀⠂⢍⣯⢊⡖⠒⣚⣃⢋⠄⢊⢄⡦⠃⠔⠉⠋⢈⡊⠦⠫⢆⡈⠁⠀⣈⠒⢘⠍⢂⡙⡙⠤⡠⡝⠌⠄⢀⎥ +⎢⠐⣒⠰⠂⠂⠈⠠⠃⢉⢾⡛⠸⡧⠂⢠⠂⣎⡵⢓⢼⠗⠒⠆⠈⠱⣦⠆⣼⠠⢐⠀⠂⢸⢠⡃⡫⢀⠳⠣⣀⎥ +⎢⠨⣀⠦⣅⠌⡄⠘⠁⠂⠬⠊⠆⠒⢰⠀⠡⠁⣉⢨⢠⢈⠤⡀⢠⣈⣥⠑⣤⣚⠄⡄⡖⠐⠩⣬⠄⠆⠀⠁⡄⎥ +⎢⡠⡈⠣⢰⣀⠂⡇⣜⡶⡵⢖⣗⢆⠤⠊⣀⢠⠦⣛⠦⡐⡸⣘⢀⢀⢂⠚⠜⠑⣤⠚⡈⢄⠈⡤⢲⠐⣔⣉⢄⎥ +⎢⠈⠠⠅⡴⠌⢴⠰⣅⠄⠁⠄⠄⠁⣩⠁⡂⠨⢱⠥⠤⠩⠥⠣⢁⠠⠀⢠⠭⡚⠠⠑⢄⠼⠷⠉⠗⣎⠅⣐⠀⎥ +⎢⢈⠔⣂⢀⠢⡫⢰⠄⠐⡌⠁⣴⡦⠠⡘⢄⠀⠒⢅⠕⡮⢢⣗⠨⠒⣒⡔⡀⡀⠑⢶⡇⡑⢌⢀⢳⡫⠀⠲⠐⎥ +⎢⢙⠭⢄⢍⠋⠲⠋⡓⡨⢅⠄⣤⢦⠢⠈⠐⠙⡔⢈⣀⠜⠯⠀⡣⡭⡨⠂⠟⢠⣋⢧⠄⢤⣐⢏⢕⡭⣍⢀⠄⎥ +⎢⢈⣤⡓⡜⡂⠴⡈⠡⠙⠲⢥⣹⡾⣺⡮⣒⡲⡭⢁⢹⠖⡨⡓⠍⢤⡐⠈⠁⢐⢤⠎⠝⠋⠊⡇⢯⠛⢄⢐⡆⎥ +⎣⠀⢁⢴⠈⠐⠉⠠⠄⠀⢠⠐⠞⠙⠀⠡⠂⢚⠆⡦⠰⢌⠀⠀⢁⠉⢢⠁⠤⠃⢜⠐⠘⢘⠂⠀⠔⠰⠴⠛⣤⎦ + +julia> A_shuffled[res.ordering, res.ordering] # A near-optimal reordering of `A_shuffled` +733×733 SparseMatrixCSC{Float64, Int64} with 2864 stored entries: +⎡⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠛⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠑⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢿⣷⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠿⣧⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢻⣶⣀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⢿⣷⡀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡄⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠉⢿⣷⣄⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢻⣶⣄⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢿⣷⎦ + +julia> bandwidth(A) +12 + +julia> bandwidth(A_shuffled) # Much larger after shuffling +703 + +julia> res # Gets very close to the original bandwidth +Results of Bandwidth Minimization Algorithm + * Algorithm: Gibbs–Poole–Stockmeyer + * Approach: heuristic + * Minimum Bandwidth: 15 + * Original Bandwidth: 703 + * Matrix Size: 733×733 ``` # Notes diff --git a/src/Minimization/core.jl b/src/Minimization/core.jl index 22dc7b4..1faa66a 100644 --- a/src/Minimization/core.jl +++ b/src/Minimization/core.jl @@ -40,41 +40,151 @@ than heuristic methods, but even more slowly). Certainly, exact solvers will always produce the same optimal bandwidth (but likely different orderings): -```@repl -using Random, SparseArrays -Random.seed!(38); -(n, p) = (20, 0.05); -A = sprand(n, n, p); -A = A .+ A' # Ensure structural symmetry -res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini()) -A[res_dcm.ordering, res_dcm.ordering] -res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough()) -A[res_sgs.ordering, res_sgs.ordering] +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(1011); + +julia> (n, p) = (20, 0.05); + +julia> A = sprand(n, n, p); + +julia> A = A .+ A' # Ensure structural symmetry +20×20 SparseMatrixCSC{Float64, Int64} with 33 stored entries: +⎡⠀⠀⠄⠂⠨⠄⢂⡐⠈⠁⎤ +⎢⠠⠁⠀⠀⠐⠀⠀⠀⠀⠀⎥ +⎢⠂⠆⠐⠀⠄⡡⠀⠁⠀⠀⎥ +⎢⢈⠰⠀⠀⠄⠀⠀⡠⠀⠀⎥ +⎣⠆⠀⠀⠀⠀⠀⠀⠀⠀⠄⎦ + +julia> res_dcm = minimize_bandwidth(A, Minimization.DelCorsoManzini()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Del Corso–Manzini + * Approach: exact + * Minimum Bandwidth: 2 + * Original Bandwidth: 18 + * Matrix Size: 20×20 + +julia> A[res_dcm.ordering, res_dcm.ordering] +20×20 SparseMatrixCSC{Float64, Int64} with 33 stored entries: +⎡⠄⡡⢄⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠑⢠⠒⢄⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠑⢠⠒⣀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠘⢠⠒⣀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠘⠁⠀⎦ + +julia> res_sgs = minimize_bandwidth(A, Minimization.SaxeGurariSudborough()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Approach: exact + * Minimum Bandwidth: 2 + * Original Bandwidth: 18 + * Matrix Size: 20×20 + +julia> A[res_sgs.ordering, res_sgs.ordering] +20×20 SparseMatrixCSC{Float64, Int64} with 33 stored entries: +⎡⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠐⢊⠐⢄⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠑⢊⠐⡄⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠉⢪⠒⣀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠘⢠⠖⎦ ``` However, the answers of (meta)heuristic solvers may differ from each other: -```@repl -using Random, SparseArrays -Random.seed!(405); -(n, p) = (400, 0.002); -A = sprand(n, n, p); -A = A .+ A' # Ensure structural symmetry; -minimize_bandwidth(A, Minimization.GibbsPooleStockmeyer()) -minimize_bandwidth(A, Minimization.ReverseCuthillMcKee()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(736362); + +julia> (n, p) = (400, 0.002); + +julia> A = sprand(n, n, p); + +julia> A = A .+ A' # Ensure structural symmetry; + +julia> minimize_bandwidth(A, Minimization.GibbsPooleStockmeyer()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Gibbs–Poole–Stockmeyer + * Approach: heuristic + * Minimum Bandwidth: 30 + * Original Bandwidth: 388 + * Matrix Size: 400×400 + +julia> minimize_bandwidth(A, Minimization.ReverseCuthillMcKee()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Reverse Cuthill–McKee + * Approach: heuristic + * Minimum Bandwidth: 31 + * Original Bandwidth: 388 + * Matrix Size: 400×400 ``` If no solver is specified, then the heuristic Gibbs–Poole–Stockmeyer algorithm is used by default: -```@repl -using Random, SparseArrays -Random.seed!(80); -(n, p) = (500, 0.001); -A = sprand(n, n, p); -A = A .+ A' # Ensure structural symmetry -res = minimize_bandwidth(A) -A[res.ordering, res.ordering] +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(2327); + +julia> (n, p) = (500, 0.001); + +julia> A = sprand(n, n, p); + +julia> A = A .+ A' # Ensure structural symmetry +500×500 SparseMatrixCSC{Float64, Int64} with 484 stored entries: +⎡⠀⠀⠀⠄⠀⠀⠀⠂⠀⠁⠁⠀⠀⠀⠈⠂⠔⠀⠀⠛⢠⡠⠀⠀⠀⠀⠀⠀⡀⠀⠀⠐⠁⠐⠠⠀⠁⠀⢀⢀⎤ +⎢⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠁⠑⠀⠀⠀⠀⠠⠁⠐⠀⠀⠀⠂⠀⠠⡠⠀⠀⠀⢐⠀⡀⣀⠀⠀⠀⠨⠀⠂⠈⎥ +⎢⠀⠀⠀⠀⠁⠀⡀⠠⢀⢂⠐⡀⠐⠐⠈⠀⠁⠀⡠⠀⢀⠀⡈⠀⠀⠀⠀⠀⠀⠨⡀⠀⠁⠐⠂⠀⠀⠀⠐⠀⎥ +⎢⠠⠀⠀⠀⠀⡈⠠⠂⠀⡀⠀⡂⠢⡀⠀⡂⠀⠀⠄⠀⠀⠀⠀⠀⠂⠀⠠⠀⢄⠐⠂⢀⠀⠁⠀⠀⢐⠀⠈⠀⎥ +⎢⠄⠀⠀⠀⠠⢐⠀⠠⠀⠀⠀⡄⠀⠀⠀⠀⠀⠠⠄⢀⠀⢂⠀⠈⠁⠀⠀⠐⠠⠀⠀⠀⠀⠀⠐⠀⠁⠀⠀⠀⎥ +⎢⠁⠀⢅⠀⠐⠠⠠⠠⠀⠤⠀⠀⠂⠀⠀⠂⠀⠀⠈⠂⠀⠀⠄⠂⠀⠀⠀⠠⡠⠀⡈⠀⠀⠀⠀⠆⢀⠂⠀⠀⎥ +⎢⠀⠀⠀⠀⢐⠀⠈⠢⠀⠀⠈⠀⠀⠀⠀⠀⠐⠀⢀⠂⠀⠀⠀⠂⠀⢀⠀⡔⠀⠀⠐⢘⠀⠠⡀⠁⠀⡀⠂⠀⎥ +⎢⠢⠀⠀⠀⠂⠀⠠⠠⠀⠀⠠⠀⠀⠀⢠⠒⠀⠀⠀⠀⠀⠀⠀⠀⠠⠀⠠⢀⠀⠀⠀⠈⠰⠁⠀⠈⠐⠀⠀⠁⎥ +⎢⠐⠁⠄⠂⠁⠀⠀⠀⠀⡀⠀⠀⠐⠀⠀⠀⠀⠀⠂⠄⠀⠀⠀⠀⠀⠀⠀⠁⠀⠐⠴⠢⠈⠀⠀⠁⠀⡀⠀⠐⎥ +⎢⣤⠀⠐⠀⠀⠊⠀⠁⠀⢁⠢⠀⠠⠐⠀⠀⠈⠄⠀⠀⠀⠀⠀⠀⠀⠀⡀⠐⠀⠀⢀⠐⠐⢀⠁⠄⠀⠀⠂⡂⎥ +⎢⠀⡲⠀⠀⠀⠐⠀⠀⠠⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠂⠀⠂⠀⠁⡄⠀⠂⠀⢀⡀⠘⠀⠀⠀⠠⎥ +⎢⠀⠀⠈⠀⠂⠈⠀⠀⡀⠀⠠⠁⠠⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠉⠀⠀⠐⠐⠀⠀⠈⠢⠀⠀⠐⠀⠀⠀⎥ +⎢⠀⠀⠀⡢⠀⠀⠈⠀⠁⠀⠀⠀⠀⢀⠀⠂⠀⠀⠀⠀⠈⠀⡄⠀⠊⠀⠀⠐⠀⡐⠀⠀⠀⠀⠆⡠⠌⠤⢰⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠂⢀⠀⠀⡀⢀⠤⠀⢂⠄⠀⢀⠈⠈⠀⠀⠀⢀⠀⠀⡠⠁⢈⠀⠀⠀⠠⠐⠂⠂⡀⠀⠀⎥ +⎢⠀⠈⢀⢀⡀⡀⢀⠑⠀⠂⠀⠊⠀⠀⠀⠀⢀⠀⠀⠀⠁⠤⢐⠀⢀⠠⡁⢀⠀⠀⠁⢀⠀⠀⠀⠄⠀⡀⠐⠀⎥ +⎢⢀⠀⠀⠠⠀⠈⠈⢀⠀⠀⠂⠈⣐⢀⡀⠀⠰⡃⢀⠐⠠⠀⠀⠀⠀⠀⠀⠀⠁⢀⠄⠁⠀⠠⢀⠀⠀⡀⡉⠀⎥ +⎢⢁⠀⠀⠘⢁⠀⠄⠀⠀⠀⠀⠀⠀⡀⠔⠂⠂⠀⠐⢀⠀⢀⠢⡀⠀⠀⠀⡀⠀⠀⠀⡀⡀⠈⠀⠁⠂⠀⡀⡀⎥ +⎢⠀⠂⠀⠀⠈⠀⠀⠀⠐⠀⠠⠄⠄⠈⡀⠀⠄⠀⠁⠄⣀⠈⠀⠀⠈⡡⠰⠀⠀⠄⠀⠐⠄⠀⠀⠀⡘⡀⠀⠀⎥ +⎢⠁⠀⠂⠂⠀⠀⠐⠐⠁⠀⠠⠐⠀⠠⠐⠀⠀⠠⠀⠀⠀⠀⠐⠀⠂⡅⠈⠠⠀⠠⠀⠠⠈⠀⠒⠨⣀⠘⠀⠀⎥ +⎣⠀⢐⡈⠀⠐⠀⠂⠀⠀⠀⠀⠀⠈⠀⠄⠀⢀⠀⠨⠠⠀⡀⠀⠀⠐⠒⠀⠀⠐⠀⠃⠈⠀⠨⠀⠀⠀⠀⠊⠀⎦ + +julia> res = minimize_bandwidth(A) +Results of Bandwidth Minimization Algorithm + * Algorithm: Gibbs–Poole–Stockmeyer + * Approach: heuristic + * Minimum Bandwidth: 6 + * Original Bandwidth: 481 + * Matrix Size: 500×500 + +julia> A[res.ordering, res.ordering] +500×500 SparseMatrixCSC{Float64, Int64} with 484 stored entries: +⎡⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤ +⎢⠀⠀⠁⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠁⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠋⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠰⣢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⢄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⣦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠻⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠱⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⢆⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠫⣦⡀⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⡀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣦⎦ ``` # Notes diff --git a/src/Recognition/core.jl b/src/Recognition/core.jl index dd7db0a..a8f04b9 100644 --- a/src/Recognition/core.jl +++ b/src/Recognition/core.jl @@ -41,28 +41,87 @@ Multiple algorithms to decide whether a given matrix has bandwidth at most `k` a available. Naturally, they will always agree, but the final orderings produced (in the case of an affirmative) may differ: -```@repl -using Random, SparseArrays -Random.seed!(52); -(n, p) = (8, 0.2); -A = sprand(Bool, n, n, p); -A = A .|| A' # Ensure structural symmetry -k = 3; -res_csg = has_bandwidth_k_ordering(A, k, Recognition.CapraraSalazarGonzalez()) -A[res_csg.ordering, res_csg.ordering] -res_sgs = has_bandwidth_k_ordering(A, k, Recognition.SaxeGurariSudborough()) -A[res_sgs.ordering, res_sgs.ordering] +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(52); + +julia> (n, p) = (8, 0.2); + +julia> A = sprand(Bool, n, n, p); + +julia> A = A .|| A' # Ensure structural symmetry +8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries: + 1 ⋅ ⋅ 1 ⋅ 1 1 ⋅ + ⋅ ⋅ ⋅ ⋅ 1 1 ⋅ ⋅ + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 + 1 ⋅ ⋅ ⋅ ⋅ ⋅ 1 ⋅ + ⋅ 1 ⋅ ⋅ ⋅ ⋅ 1 1 + 1 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + 1 ⋅ ⋅ 1 1 ⋅ ⋅ 1 + ⋅ ⋅ 1 ⋅ 1 ⋅ 1 1 + +julia> k = 3; + +julia> res_csg = has_bandwidth_k_ordering(A, k, Recognition.CapraraSalazarGonzalez()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Caprara–Salazar-González + * Bandwidth Threshold k: 3 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 6 + * Matrix Size: 8×8 + +julia> A[res_csg.ordering, res_csg.ordering] +8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries: + ⋅ ⋅ 1 1 ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ 1 ⋅ 1 ⋅ ⋅ ⋅ + 1 1 1 1 ⋅ ⋅ ⋅ ⋅ + 1 ⋅ 1 ⋅ ⋅ 1 1 ⋅ + ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ + ⋅ ⋅ ⋅ 1 1 ⋅ 1 ⋅ + ⋅ ⋅ ⋅ 1 ⋅ 1 1 1 + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 ⋅ + +julia> res_sgs = has_bandwidth_k_ordering(A, k, Recognition.SaxeGurariSudborough()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Bandwidth Threshold k: 3 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 6 + * Matrix Size: 8×8 + +julia> A[res_sgs.ordering, res_sgs.ordering] +8×8 SparseMatrixCSC{Bool, Int64} with 22 stored entries: + ⋅ 1 ⋅ 1 ⋅ ⋅ ⋅ ⋅ + 1 ⋅ 1 ⋅ 1 ⋅ ⋅ ⋅ + ⋅ 1 1 ⋅ 1 1 ⋅ ⋅ + 1 ⋅ ⋅ ⋅ ⋅ ⋅ 1 ⋅ + ⋅ 1 1 ⋅ ⋅ ⋅ 1 1 + ⋅ ⋅ 1 ⋅ ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ 1 1 ⋅ 1 1 + ⋅ ⋅ ⋅ ⋅ 1 ⋅ 1 ⋅ ``` If no decider is specified, then the Caprara–Salazar-González algorithm is used by default: -```@repl -using Random, SparseArrays -Random.seed!(174); -(n, p, k) = (20, 0.1, 4); -A = sprand(n, n, p); -A = A .+ A' # Ensure structural symmetry; -has_bandwidth_k_ordering(A, k) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(174); + +julia> (n, p, k) = (20, 0.1, 4); + +julia> A = sprand(n, n, p); + +julia> A = A .+ A' # Ensure structural symmetry; + +julia> has_bandwidth_k_ordering(A, k) +Results of Bandwidth Recognition Algorithm + * Algorithm: Caprara–Salazar-González + * Bandwidth Threshold k: 4 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 15 + * Matrix Size: 20×20 ``` # Notes diff --git a/src/Recognition/deciders/brute_force_search.jl b/src/Recognition/deciders/brute_force_search.jl index f2220d0..9276d87 100644 --- a/src/Recognition/deciders/brute_force_search.jl +++ b/src/Recognition/deciders/brute_force_search.jl @@ -35,15 +35,34 @@ minimum) possible permutations—in these cases, it is infeasible to go above `` ``10×10`` without incurring multiple-hour runtimes. (Even when ``k`` is considerably larger than the true minimum, it is unlikely that a bandwidth-``k`` ordering will be found in a reasonable time frame.) Nevertheless, we see that it is quite effective for, say, ``8×8``: -```@repl -using Random, SparseArrays -Random.seed!(314159); -(n, p) = (8, 0.5); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -(k_false, k_true) = (3, 5); -has_bandwidth_k_ordering(A, k_false, Recognition.BruteForceSearch()) -has_bandwidth_k_ordering(A, k_true, Recognition.BruteForceSearch()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(314159); + +julia> (n, p) = (8, 0.5); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> (k_false, k_true) = (3, 5); + +julia> has_bandwidth_k_ordering(A, k_false, Recognition.BruteForceSearch()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Brute-force search + * Bandwidth Threshold k: 3 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 6 + * Matrix Size: 8×8 + +julia> has_bandwidth_k_ordering(A, k_true, Recognition.BruteForceSearch()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Brute-force search + * Bandwidth Threshold k: 5 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 6 + * Matrix Size: 8×8 ``` # Notes diff --git a/src/Recognition/deciders/caprara_salazar_gonzalez.jl b/src/Recognition/deciders/caprara_salazar_gonzalez.jl index 6b3cdfe..153e380 100644 --- a/src/Recognition/deciders/caprara_salazar_gonzalez.jl +++ b/src/Recognition/deciders/caprara_salazar_gonzalez.jl @@ -39,15 +39,34 @@ efficient pruning techniques and compatibility checks result in approximately ex growth in time complexity with respect to ``n``. # Examples -```@repl -using Random, SparseArrays -Random.seed!(17); -(n, p) = (10, 0.17); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -(k_false, k_true) = (3, 6); -has_bandwidth_k_ordering(A, k_false, Recognition.CapraraSalazarGonzalez()) -has_bandwidth_k_ordering(A, k_true, Recognition.CapraraSalazarGonzalez()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(17); + +julia> (n, p) = (10, 0.17); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> (k_false, k_true) = (3, 6); + +julia> has_bandwidth_k_ordering(A, k_false, Recognition.CapraraSalazarGonzalez()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Caprara–Salazar-González + * Bandwidth Threshold k: 3 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 9 + * Matrix Size: 10×10 + +julia> has_bandwidth_k_ordering(A, k_true, Recognition.CapraraSalazarGonzalez()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Caprara–Salazar-González + * Bandwidth Threshold k: 6 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 9 + * Matrix Size: 10×10 ``` # Notes diff --git a/src/Recognition/deciders/del_corso_manzini.jl b/src/Recognition/deciders/del_corso_manzini.jl index 403e44f..8d71419 100644 --- a/src/Recognition/deciders/del_corso_manzini.jl +++ b/src/Recognition/deciders/del_corso_manzini.jl @@ -37,15 +37,34 @@ Based on experimental results, the algorithm is feasible for ``n×n`` matrices u ``n ≈ 100`` or so. # Examples -```@repl -using Random, SparseArrays -Random.seed!(7878); -(n, p) = (40, 0.1); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -(k_false, k_true) = (13, 26); -has_bandwidth_k_ordering(A, k_false, Recognition.DelCorsoManzini()) -has_bandwidth_k_ordering(A, k_true, Recognition.DelCorsoManzini()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(7878); + +julia> (n, p) = (40, 0.1); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> (k_false, k_true) = (13, 26); + +julia> has_bandwidth_k_ordering(A, k_false, Recognition.DelCorsoManzini()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Del Corso–Manzini + * Bandwidth Threshold k: 13 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 34 + * Matrix Size: 40×40 + +julia> has_bandwidth_k_ordering(A, k_true, Recognition.DelCorsoManzini()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Del Corso–Manzini + * Bandwidth Threshold k: 26 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 34 + * Matrix Size: 40×40 ``` # Notes @@ -136,30 +155,54 @@ Based on experimental results, the algorithm is feasible for ``n×n`` matrices u # Examples Here, Del Corso–Manzini with perimeter search ascertains that A random ``30×30`` matrix has -a minimum bandwidth greater than ``9``. The depth parameter is not explicitly set; instead, +a minimum bandwidth greater than ``6``. The depth parameter is not explicitly set; instead, some near-optimal value is automatically computed upon the first [`has_bandwidth_k_ordering`](@ref) function call. -```@repl -using Random, SparseArrays -Random.seed!(5847); -(n, p) = (30, 0.05); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -k = 6; -has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(5847); + +julia> (n, p) = (30, 0.05); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> k = 6; + +julia> has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Del Corso–Manzini with perimeter search + * Bandwidth Threshold k: 6 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 27 + * Matrix Size: 30×30 ``` Now, Del Corso–Manzini with perimeter search recognizes that a random ``35×35`` matrix has a minimum bandwidth at most ``8``. In this case, we explicitly set the depth parameter to ``4`` beforehand instead of relying on [`Recognition.dcm_ps_optimal_depth`](@ref). -```@repl -using Random, SparseArrays -Random.seed!(23552); -(n, p, depth) = (35, 0.02, 4); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -k = 8; -has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS(depth)) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(23552); + +julia> (n, p, depth) = (35, 0.02, 4); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> k = 8; + +julia> has_bandwidth_k_ordering(A, k, Recognition.DelCorsoManziniWithPS(depth)) +Results of Bandwidth Recognition Algorithm + * Algorithm: Del Corso–Manzini with perimeter search + * Bandwidth Threshold k: 8 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 32 + * Matrix Size: 35×35 ``` # Notes diff --git a/src/Recognition/deciders/saxe_gurari_sudborough.jl b/src/Recognition/deciders/saxe_gurari_sudborough.jl index 48f5035..b398c70 100644 --- a/src/Recognition/deciders/saxe_gurari_sudborough.jl +++ b/src/Recognition/deciders/saxe_gurari_sudborough.jl @@ -43,15 +43,34 @@ for larger ``k``, given that their aggressive pruning strategies keep their effe space very small in practice. # Examples -```@repl -using Random, SparseArrays -Random.seed!(274); -(n, p) = (20, 0.08); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -(k_false, k_true) = (3, 5); -has_bandwidth_k_ordering(A, k_false, Recognition.SaxeGurariSudborough()) -has_bandwidth_k_ordering(A, k_true, Recognition.SaxeGurariSudborough()) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(274); + +julia> (n, p) = (20, 0.08); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> (k_false, k_true) = (3, 5); + +julia> has_bandwidth_k_ordering(A, k_false, Recognition.SaxeGurariSudborough()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Bandwidth Threshold k: 3 + * Has Bandwidth ≤ k Ordering: false + * Original Bandwidth: 12 + * Matrix Size: 20×20 + +julia> has_bandwidth_k_ordering(A, k_true, Recognition.SaxeGurariSudborough()) +Results of Bandwidth Recognition Algorithm + * Algorithm: Saxe–Gurari–Sudborough + * Bandwidth Threshold k: 5 + * Has Bandwidth ≤ k Ordering: true + * Original Bandwidth: 12 + * Matrix Size: 20×20 ``` # Notes diff --git a/src/core.jl b/src/core.jl index 94db9a3..b52ce3d 100644 --- a/src/core.jl +++ b/src/core.jl @@ -27,15 +27,42 @@ Given an ``n×n`` input matrix, this relatively simple algorithm runs in ``O(n² # Examples `bandwidth` correctly identifies the bandwidth of a pentadiagonal matrix as ``2`` and does not attempt to find a minimizing permutation upon shuffling of its rows and columns: -```@repl -using Random -Random.seed!(242622); -(n, k) = (8, 2); -perm = randperm(n); -A = (!iszero).(MatrixBandwidth.random_banded_matrix(n, k)) -bandwidth(A) -A_shuffled = A[perm, perm] -bandwidth(A_shuffled) +```jldoctest +julia> using Random + +julia> Random.seed!(242622); + +julia> (n, k) = (8, 2); + +julia> perm = randperm(n); + +julia> A = (!iszero).(MatrixBandwidth.random_banded_matrix(n, k)) +8×8 BitMatrix: + 1 0 0 0 0 0 0 0 + 0 1 0 1 0 0 0 0 + 0 0 0 1 1 0 0 0 + 0 1 1 1 0 1 0 0 + 0 0 1 0 0 0 0 0 + 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + +julia> bandwidth(A) +2 + +julia> A_shuffled = A[perm, perm] +8×8 BitMatrix: + 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 + 0 0 0 1 0 0 0 0 + 0 0 1 0 0 0 1 0 + 0 0 0 0 1 0 0 0 + 0 0 0 0 0 1 1 0 + 1 0 0 1 0 1 1 0 + 0 0 0 0 0 0 0 0 + +julia> bandwidth(A_shuffled) +6 ``` # Notes @@ -99,23 +126,50 @@ Given an ``n×n`` input matrix, this relatively simple algorithm runs in ``O(n² # Examples `profile` computes the column profile of a matrix by default: -```@repl -using Random, SparseArrays -Random.seed!(2287); -(n, p) = (25, 0.05); -A = sprand(n, n, p) -profile(A) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(2287); + +julia> (n, p) = (25, 0.05); + +julia> A = sprand(n, n, p) +25×25 SparseMatrixCSC{Float64, Int64} with 29 stored entries: +⎡⠀⠀⠀⠀⠀⠀⠐⠀⠒⠀⡀⠀⠀⎤ +⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠐⠂⎥ +⎢⠀⠀⠀⢀⠌⠀⠀⠀⢀⠈⠀⠀⠀⎥ +⎢⠠⢄⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥ +⎢⠀⠀⠀⠂⡀⠀⠀⠌⠀⠈⠀⠀⠀⎥ +⎢⠀⠀⠀⠀⠀⠀⠢⡀⢄⡈⠀⠀⠀⎥ +⎣⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦ + +julia> profile(A) 211 ``` The dimension (`:row` or `:col`) can also be explicitly specified: -```@repl -using Random, SparseArrays -Random.seed!(3647); -(n, p) = (25, 0.05); -A = sprand(n, n, p) -profile(A, dim=:row) -profile(A, dim=:col) +```jldoctest +julia> using Random, SparseArrays + +julia> Random.seed!(3647); + +julia> (n, p) = (25, 0.05); + +julia> A = sprand(n, n, p) +25×25 SparseMatrixCSC{Float64, Int64} with 31 stored entries: +⎡⠄⠀⠀⠀⠀⠀⠀⠘⠀⠀⠀⠁⠀⎤ +⎢⠄⢀⠀⠀⠁⠀⠀⠀⠀⢀⠀⠀⠀⎥ +⎢⠀⢀⡂⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⎥ +⎢⠀⠀⠀⠀⠀⡀⠂⠀⠀⠀⠀⠀⠀⎥ +⎢⠁⠀⠁⠀⠁⠀⠀⠁⠄⢀⠈⠀⠀⎥ +⎢⠂⠐⠐⠐⠠⠀⠄⠀⠀⠀⠠⣀⠀⎥ +⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎦ + +julia> profile(A, dim=:row) +147 + +julia> profile(A, dim=:col) +175 ``` # References @@ -190,14 +244,27 @@ in ``O(n²)`` time). cases but not universally so.) # Examples -```@repl -using Random, SparseArrays, Combinatorics -Random.seed!(21); -(n, p) = (9, 0.4); -A = sprand(n, n, p); -A = A + A' # Ensure structural symmetry; -minimize_bandwidth(A, Minimization.BruteForceSearch()) -bandwidth_lower_bound(A) # Always less than or equal to the true minimum bandwidth +```jldoctest +julia> using Random, SparseArrays, Combinatorics + +julia> Random.seed!(21); + +julia> (n, p) = (9, 0.4); + +julia> A = sprand(n, n, p); + +julia> A = A + A' # Ensure structural symmetry; + +julia> minimize_bandwidth(A, Minimization.BruteForceSearch()) +Results of Bandwidth Minimization Algorithm + * Algorithm: Brute-force search + * Approach: exact + * Minimum Bandwidth: 5 + * Original Bandwidth: 8 + * Matrix Size: 9×9 + +julia> bandwidth_lower_bound(A) # Always less than or equal to the true minimum bandwidth +4 ``` # Notes diff --git a/src/utils.jl b/src/utils.jl index fefb04e..d3cbf5f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -34,24 +34,57 @@ permutation inducing a minimum bandwidth less than `k`, especially for small val # Examples Generate a ``6×6`` matrix with bandwidth ``1`` and the maximum number of nonzero entries: -```@repl -using Random -A = MatrixBandwidth.random_banded_matrix(6, 1; p=1, rng=MersenneTwister(1228)) -bandwidth(A) +```jldoctest +julia> using Random + +julia> A = MatrixBandwidth.random_banded_matrix(6, 1; p=1, rng=MersenneTwister(1228)) +6×6 Matrix{Float64}: + 0.310239 0.346413 0.0 0.0 0.0 0.0 + 0.509981 0.917073 0.390771 0.0 0.0 0.0 + 0.0 0.760045 0.808396 0.0195686 0.0 0.0 + 0.0 0.0 0.222338 0.853164 0.806888 0.0 + 0.0 0.0 0.0 0.421603 0.132165 0.805813 + 0.0 0.0 0.0 0.0 0.305339 0.0799183 + +julia> bandwidth(A) +1 ``` Generate a ``7×7`` matrix with bandwidth ``3`` and band density ``0.3``: -```@repl -using Random -A = MatrixBandwidth.random_banded_matrix(7, 3; p=0.3, rng=MersenneTwister(0402)) -bandwidth(A) +```jldoctest +julia> using Random + +julia> A = MatrixBandwidth.random_banded_matrix(7, 3; p=0.3, rng=MersenneTwister(0402)) +7×7 Matrix{Float64}: + 0.0 0.132699 0.0 0.0 0.0 0.0 0.0 + 0.869352 0.0 0.324319 0.926496 0.0 0.0 0.0 + 0.0 0.891878 0.0 0.658102 0.0 0.0 0.0 + 0.0 0.88859 0.399559 0.0 0.0 0.284285 0.703377 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.489594 0.0 0.0 0.393573 + 0.0 0.0 0.0 0.412412 0.0 0.47063 0.0 + +julia> bandwidth(A) +3 ``` Generate an ``8×8`` diagonal (bandwidth ``0``) matrix with default band density (``0.5``): -```@repl -using Random -A = MatrixBandwidth.random_banded_matrix(8, 0; rng=MersenneTwister(0102)) -bandwidth(A) +```jldoctest +julia> using Random + +julia> A = MatrixBandwidth.random_banded_matrix(8, 0; rng=MersenneTwister(0102)) +8×8 Matrix{Float64}: + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0762399 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.373113 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.726309 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.41974 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.293132 + +julia> bandwidth(A) +0 ``` # Notes @@ -127,11 +160,29 @@ Find the indices of all connected components of the graph whose adjacency matrix indices belonging to a connected component. # Examples -```@repl -using Graphs -g = complement(complete_multipartite_graph([3, 4, 2])) -A = Bool.(adjacency_matrix(g)) -MatrixBandwidth.connected_components(A) +```jldoctest +julia> using Graphs + +julia> g = complement(complete_multipartite_graph([3, 4, 2])) +{9, 10} undirected simple Int64 graph + +julia> A = Bool.(adjacency_matrix(g)) +9×9 SparseArrays.SparseMatrixCSC{Bool, Int64} with 20 stored entries: + ⋅ 1 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + 1 ⋅ 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + 1 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ ⋅ 1 1 1 ⋅ ⋅ + ⋅ ⋅ ⋅ 1 ⋅ 1 1 ⋅ ⋅ + ⋅ ⋅ ⋅ 1 1 ⋅ 1 ⋅ ⋅ + ⋅ ⋅ ⋅ 1 1 1 ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1 ⋅ + +julia> MatrixBandwidth.connected_components(A) +3-element Vector{Vector{Int64}}: + [1, 2, 3] + [4, 5, 6, 7] + [8, 9] ``` """ function connected_components(A::AbstractMatrix{Bool}) @@ -189,19 +240,67 @@ that each level in the triple-nested loop iterates over ``O(n)`` entries. # Examples Floyd–Warshall finds the shortest distances between all pairs of nodes in a connected graph: -```@repl -using Graphs -g = ladder_graph(5) -A = Bool.(adjacency_matrix(g)) -Int.(MatrixBandwidth.floyd_warshall_shortest_paths(A)) +```jldoctest +julia> using Graphs + +julia> g = ladder_graph(5) +{10, 13} undirected simple Int64 graph + +julia> A = Bool.(adjacency_matrix(g)) +10×10 SparseArrays.SparseMatrixCSC{Bool, Int64} with 26 stored entries: + ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ ⋅ + 1 ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ + ⋅ 1 ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ ⋅ + ⋅ ⋅ 1 ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ + ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ ⋅ ⋅ 1 + 1 ⋅ ⋅ ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ + ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ 1 ⋅ ⋅ + ⋅ ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ 1 ⋅ + ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ 1 + ⋅ ⋅ ⋅ ⋅ 1 ⋅ ⋅ ⋅ 1 ⋅ + +julia> Int.(MatrixBandwidth.floyd_warshall_shortest_paths(A)) +10×10 Matrix{Int64}: + 0 1 2 3 4 1 2 3 4 5 + 1 0 1 2 3 2 1 2 3 4 + 2 1 0 1 2 3 2 1 2 3 + 3 2 1 0 1 4 3 2 1 2 + 4 3 2 1 0 5 4 3 2 1 + 1 2 3 4 5 0 1 2 3 4 + 2 1 2 3 4 1 0 1 2 3 + 3 2 1 2 3 2 1 0 1 2 + 4 3 2 1 2 3 2 1 0 1 + 5 4 3 2 1 4 3 2 1 0 ``` Floyd–Warshall assigns `Inf` to pairs of nodes in different connected components: -```@repl -using Graphs -g = complement(wheel_graph(8)) -A = Bool.(adjacency_matrix(g)) -MatrixBandwidth.floyd_warshall_shortest_paths(A) +```jldoctest +julia> using Graphs + +julia> g = complement(wheel_graph(8)) +{8, 14} undirected simple Int64 graph + +julia> A = Bool.(adjacency_matrix(g)) +8×8 SparseArrays.SparseMatrixCSC{Bool, Int64} with 28 stored entries: + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ + ⋅ ⋅ ⋅ 1 1 1 1 ⋅ + ⋅ ⋅ ⋅ ⋅ 1 1 1 1 + ⋅ 1 ⋅ ⋅ ⋅ 1 1 1 + ⋅ 1 1 ⋅ ⋅ ⋅ 1 1 + ⋅ 1 1 1 ⋅ ⋅ ⋅ 1 + ⋅ 1 1 1 1 ⋅ ⋅ ⋅ + ⋅ ⋅ 1 1 1 1 ⋅ ⋅ + +julia> MatrixBandwidth.floyd_warshall_shortest_paths(A) +8×8 Matrix{Float64}: + 0.0 Inf Inf Inf Inf Inf Inf Inf + Inf 0.0 2.0 1.0 1.0 1.0 1.0 2.0 + Inf 2.0 0.0 2.0 1.0 1.0 1.0 1.0 + Inf 1.0 2.0 0.0 2.0 1.0 1.0 1.0 + Inf 1.0 1.0 2.0 0.0 2.0 1.0 1.0 + Inf 1.0 1.0 1.0 2.0 0.0 2.0 1.0 + Inf 1.0 1.0 1.0 1.0 2.0 0.0 2.0 + Inf 2.0 1.0 1.0 1.0 1.0 2.0 0.0 ``` """ function floyd_warshall_shortest_paths(A::AbstractMatrix{Bool}) @@ -239,11 +338,25 @@ Check whether `A[i, j]` is nonzero if and only if `A[j, i]` is nonzero for all ` - `::Bool`: whether `A` is structurally symmetric. # Examples -```@repl -A = [4 0 9 -2; 0 0 1 0; 3 -1 5 0; 4 0 0 3] -MatrixBandwidth.is_structurally_symmetric(A) -B = [1.12 2.36 0.00; 5.99 0.0 0.0; 0.0 3.1 -7.49] -MatrixBandwidth.is_structurally_symmetric(B) +```jldoctest +julia> A = [4 0 9 -2; 0 0 1 0; 3 -1 5 0; 4 0 0 3] +4×4 Matrix{Int64}: + 4 0 9 -2 + 0 0 1 0 + 3 -1 5 0 + 4 0 0 3 + +julia> MatrixBandwidth.is_structurally_symmetric(A) +true + +julia> B = [1.12 2.36 0.00; 5.99 0.0 0.0; 0.0 3.1 -7.49] +3×3 Matrix{Float64}: + 1.12 2.36 0.0 + 5.99 0.0 0.0 + 0.0 3.1 -7.49 + +julia> MatrixBandwidth.is_structurally_symmetric(B) +false ``` # Notes @@ -269,9 +382,20 @@ Convert `A` to a boolean matrix and set all its diagonal entries to `false`. `false`. # Examples -```@repl -A = [0 2 0 7; 0 -8 0 3; -1 9 0 0; 0 0 0 5] -MatrixBandwidth.offdiag_nz_support(A) +```jldoctest +julia> A = [0 2 0 7; 0 -8 0 3; -1 9 0 0; 0 0 0 5] +4×4 Matrix{Int64}: + 0 2 0 7 + 0 -8 0 3 + -1 9 0 0 + 0 0 0 5 + +julia> MatrixBandwidth.offdiag_nz_support(A) +4×4 BitMatrix: + 0 1 0 1 + 0 0 0 1 + 1 1 0 0 + 0 0 0 0 ``` # Notes @@ -298,18 +422,33 @@ Identify the highest supertype of `subtype` that is also a subtype of `abstractt - `::Type`: the direct subtype of `abstracttype` that is a supertype of `subtype`. # Examples -```@repl -abstract type Parent end -abstract type Child1 <: Parent end -abstract type Grandchild1 <: Child1 end -struct Grandchild2 <: Child1 end -abstract type Child2 <: Parent end -struct Child3 <: Parent end -MatrixBandwidth.find_direct_subtype(Parent, Child1) -MatrixBandwidth.find_direct_subtype(Parent, Grandchild1) -MatrixBandwidth.find_direct_subtype(Parent, Grandchild2) -MatrixBandwidth.find_direct_subtype(Parent, Child2) -MatrixBandwidth.find_direct_subtype(Parent, Child3) +```jldoctest +julia> abstract type Parent end + +julia> abstract type Child1 <: Parent end + +julia> abstract type Grandchild1 <: Child1 end + +julia> struct Grandchild2 <: Child1 end + +julia> abstract type Child2 <: Parent end + +julia> struct Child3 <: Parent end + +julia> MatrixBandwidth.find_direct_subtype(Parent, Child1) +Child1 + +julia> MatrixBandwidth.find_direct_subtype(Parent, Grandchild1) +Child1 + +julia> MatrixBandwidth.find_direct_subtype(Parent, Grandchild2) +Child1 + +julia> MatrixBandwidth.find_direct_subtype(Parent, Child2) +Child2 + +julia> MatrixBandwidth.find_direct_subtype(Parent, Child3) +Child3 ``` """ function find_direct_subtype(abstracttype::Type, subtype::Type)